text:
"""
$(SIGNATURES)

Kinetic flux vector splitting (KFVS) flux
"""
function flux_kfvs!(
    face::Interface1F,
    ctrL::T,
    ctrR::T,
    gas::AbstractProperty,
    vs::AbstractVelocitySpace1D,
    p,
    dt = 1.0,
) where {T<:ControlVolume1F}

    dxL, dxR = p[1:2]

    flux_kfvs!(
        face.fw,
        face.ff,
        ctrL.f + ctrL.sf * dxL,
        ctrR.f - ctrR.sf * dxR,
        vs.u,
        vs.weights,
        dt,
        ctrL.sf,
        ctrR.sf,
    )

    return nothing

end

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

    dxL, dxR = p[1:2]

    flux_kfvs!(
        face.fw,
        face.fh,
        face.fb,
        ctrL.h .+ ctrL.sb .* dxL,
        ctrL.b .+ ctrL.sb .* dxL,
        ctrR.h .- ctrR.sb .* dxR,
        ctrR.b .- ctrR.sb .* dxR,
        vs.u,
        vs.weights,
        dt,
        ctrL.sh,
        ctrL.sb,
        ctrR.sh,
        ctrR.sb,
    )

    return nothing

end

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

    dxL, dxR, len, n, dirc = p[1:5]
    sfL = @view ctrL.sf[:, :, dirc]
    sfR = @view ctrR.sf[:, :, dirc]

    flux_kfvs!(
        face.fw,
        face.ff,
        ctrL.f + sfL * dxL,
        ctrR.f - sfR * dxR,
        vs.u .* n[1] .+ vs.v .* n[2],
        vs.v .* n[1] .- vs.u .* n[2],
        vs.weights,
        dt,
        len,
        sfL,
        sfR,
    )

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

    return nothing

end

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

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

    flux_kfvs!(
        face.fw,
        face.fh,
        face.fb,
        ctrL.h .+ shL .* dxL,
        ctrL.b .+ sbL .* dxL,
        ctrR.h .- shR .* dxR,
        ctrR.b .- sbR .* dxR,
        vs.u .* n[1] .+ vs.v .* n[2],
        vs.v .* n[1] .- vs.u .* n[2],
        vs.weights,
        dt,
        len,
        shL,
        sbL,
        shR,
        sbR,
    )

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

    return nothing

end

function flux_kfvs!(
    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_kfvs!(
            face.fw,
            face.ff,
            ctrL.f + ctrL.sf * dxL,
            ctrR.f - ctrR.sf * dxR,
            vs.u,
            vs.v,
            vs.w,
            vs.weights,
            dt,
            1.0,
            ctrL.sf,
            ctrR.sf,
        )
    else
        len, n, dirc = p[3:5]
        sfL = @view ctrL.sf[:, :, :, dirc]
        sfR = @view ctrR.sf[:, :, :, dirc]

        flux_kfvs!(
            face.fw,
            face.ff,
            ctrL.f + sfL * dxL,
            ctrR.f - sfR * dxR,
            vs.u .* n[1] .+ vs.v .* n[2],
            vs.v .* n[1] .- vs.u .* n[2],
            vs.w,
            vs.weights,
            dt,
            len,
            sfL,
            sfR,
        )
        face.fw .= global_frame(face.fw, n[1], n[2])
    end

    return nothing

end

# ------------------------------------------------------------
# Low-level backends
# ------------------------------------------------------------

"""
$(SIGNATURES)

Kinetic flux vector splitting (KFVS) flux
"""
function flux_kfvs(
    fL::Y,
    fR::Y,
    u::AV{<:FN},
    dt,
    sfL = zero(fL)::Y,
    sfR = zero(fR)::Y,
) where {Y<:AV{<:FN}}

    ff = similar(fL)
    flux_kfvs!(ff, fL, fR, u, dt, sfL, sfR)

    return ff

end

"""
$(SIGNATURES)

Kinetic flux vector splitting (KFVS) flux

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

1F1V for pure DOM
"""
function flux_kfvs!(
    ff::AV{<:FN},
    fL::Y,
    fR::Y,
    u::AV{<:FN},
    dt,
    sfL = zero(fL)::Y,
    sfR = zero(fR)::Y,
) where {Y<:AV{<:FN}} # 1F1V flux for pure DOM

    # upwind reconstruction
    δ = heaviside.(u)

    f = @. fL * δ + fR * (1.0 - δ)
    sf = @. sfL * δ + sfR * (1.0 - δ)

    # calculate flux
    @. ff = dt * u * f - 0.5 * dt^2 * u^2 * sf

    return nothing

end

"""
$(SIGNATURES)

1F1V
"""
function flux_kfvs!(
    fw::AV{<:FN},
    ff::AV{<:FN},
    fL::Z,
    fR::Z,
    u::A,
    ω::A,
    dt,
    sfL = zero(fL)::Z,
    sfR = zero(fR)::Z,
) where {Z<:AV{<:FN},A<:AV{<:FN}}

    # upwind reconstruction
    δ = heaviside.(u)

    f = @. fL * δ + fR * (1.0 - δ)
    sf = @. sfL * δ + sfR * (1.0 - δ)

    # calculate fluxes
    @. ff = dt * u * f - 0.5 * dt^2 * u^2 * sf
    flux_conserve!(fw, ff, u, ω)
    #fw[1] = dt * sum(ω .* u .* f) - 0.5 * dt^2 * sum(ω .* u .^ 2 .* sf)
    #fw[2] = dt * sum(ω .* u .^ 2 .* f) - 0.5 * dt^2 * sum(ω .* u .^ 3 .* sf)
    #fw[3] = dt * 0.5 * sum(ω .* u .^ 3 .* f) - 0.5 * dt^2 * 0.5 * sum(ω .* u .^ 4 .* sf)

    return nothing

end

"""
$(SIGNATURES)

Mixture
"""
function flux_kfvs!(
    fw::AM{<:FN},
    ff::AM{<:FN},
    fL::Z,
    fR::Z,
    u::A,
    ω::A,
    dt,
    sfL = zero(fL)::Z,
    sfR = zero(fR)::Z,
) where {Z<:AM{<:FN},A<:AM{<:FN}}

    for j in axes(fw, 2)
        _fw = @view fw[:, j]
        _ff = @view ff[:, j]
        _fL = @view fL[:, j]
        _fR = @view fR[:, j]
        _u = @view u[:, j]
        _ω = @view ω[:, j]
        _sfL = @view sfL[:, j]
        _sfR = @view sfR[:, j]

        flux_kfvs!(_fw, _ff, _fL, _fR, _u, _ω, dt, _sfL, _sfR)
    end

    return nothing

end

"""
$(SIGNATURES)

2F1V
"""
function flux_kfvs!(
    fw::AV{<:FN},
    fh::Y,
    fb::Y,
    hL::Z,
    bL::Z,
    hR::Z,
    bR::Z,
    u::A,
    ω::A,
    dt,
    shL = zero(hL)::Z,
    sbL = zero(bL)::Z,
    shR = zero(hR)::Z,
    sbR = zero(bR)::Z,
) where {Y<:AV{<:FN},Z<:AV{<:FN},A<:AV{<:FN}}

    # upwind reconstruction
    δ = heaviside.(u)

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

    sh = @. shL * δ + shR * (1.0 - δ)
    sb = @. sbL * δ + sbR * (1.0 - δ)

    # calculate fluxes
    @. fh = dt * u * h - 0.5 * dt^2 * u^2 * sh
    @. fb = dt * u * b - 0.5 * dt^2 * u^2 * sb
    flux_conserve!(fw, fh, fb, u, ω)
    #=fw[1] = dt * sum(ω .* u .* h) - 0.5 * dt^2 * sum(ω .* u .^ 2 .* sh)
    fw[2] = dt * sum(ω .* u .^ 2 .* h) - 0.5 * dt^2 * sum(ω .* u .^ 3 .* sh)
    fw[3] =
        dt * 0.5 * (sum(ω .* u .^ 3 .* h) + sum(ω .* u .* b)) -
        0.5 * dt^2 * 0.5 * (sum(ω .* u .^ 4 .* sh) + sum(ω .* u .^ 2 .* sb))=#

    return nothing

end

"""
$(SIGNATURES)

Mixture
"""
function flux_kfvs!(
    fw::AM{<:FN},
    fh::Y,
    fb::Y,
    hL::Z,
    bL::Z,
    hR::Z,
    bR::Z,
    u::A,
    ω::A,
    dt,
    shL = zero(hL)::Z,
    sbL = zero(bL)::Z,
    shR = zero(hR)::Z,
    sbR = zero(bR)::Z,
) where {Y<:AM{<:FN},Z<:AM{<:FN},A<:AM{<:FN}}

    for j in axes(fw, 2)
        _fw = @view fw[:, j]
        _fh = @view fh[:, j]
        _fb = @view fb[:, j]
        _hL = @view hL[:, j]
        _bL = @view bL[:, j]
        _hR = @view hR[:, j]
        _bR = @view bR[:, j]
        _u = @view u[:, j]
        _ω = @view ω[:, j]
        _shL = @view shL[:, j]
        _sbL = @view sbL[:, j]
        _shR = @view shR[:, j]
        _sbR = @view sbR[:, j]

        flux_kfvs!(_fw, _fh, _fb, _hL, _bL, _hR, _bR, _u, _ω, dt, _shL, _sbL, _shR, _sbR)
    end

    return nothing

end

"""
$(SIGNATURES)

3F1V @ Rykov
"""
function flux_kfvs!(
    fw::AV{<:FN},
    fh::Y,
    fb::Y,
    fr::Y,
    hL::Z,
    bL::Z,
    rL::Z,
    hR::Z,
    bR::Z,
    rR::Z,
    u::A,
    ω::A,
    dt,
    shL = zero(hL)::Z,
    sbL = zero(bL)::Z,
    srL = zero(rL)::Z,
    shR = zero(hR)::Z,
    sbR = zero(bR)::Z,
    srR = zero(rR)::Z,
) where {Y<:AV{<:FN},Z<:AV{<:FN},A<:AV{<:FN}}

    # upwind reconstruction
    δ = heaviside.(u)

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

    sh = @. shL * δ + shR * (1.0 - δ)
    sb = @. sbL * δ + sbR * (1.0 - δ)
    sr = @. srL * δ + srR * (1.0 - δ)

    # macro fluxes
    fw[1] = dt * sum(ω .* u .* h) - 0.5 * dt^2 * sum(ω .* u .^ 2 .* sh)
    fw[2] = dt * sum(ω .* u .^ 2 .* h) - 0.5 * dt^2 * sum(ω .* u .^ 3 .* sh)
    fw[3] =
        dt * 0.5 * (sum(ω .* u .^ 3 .* h) + sum(ω .* u .* b)) -
        0.5 * dt^2 * 0.5 * (sum(ω .* u .^ 4 .* sh) + sum(ω .* u .^ 2 .* sb))
    fw[4] = dt * 0.5 * sum(ω .* u .* r) - 0.5 * dt^2 * sum(ω .* u .^ 2 .* sr)

    # micro fluxes
    @. fh = dt * u * h - 0.5 * dt^2 * u^2 * sh
    @. fb = dt * u * b - 0.5 * dt^2 * u^2 * sb
    @. fr = dt * u * r - 0.5 * dt^2 * u^2 * sr

    return nothing

end

"""
$(SIGNATURES)

1F3V
"""
function flux_kfvs!(
    fw::AV{<:FN},
    ff::AA{<:FN,3},
    fL::Z,
    fR::Z,
    u::A,
    v::A,
    w::A,
    ω::A,
    dt::Real,
    len = 1.0::Real,
    sfL = zero(fL)::Z,
    sfR = zero(fR)::Z,
) where {Z<:AA{<:FN,3},A<:AA{<:FN,3}}

    # upwind reconstruction
    δ = heaviside.(u)

    f = @. fL * δ + fR * (1.0 - δ)
    sf = @. sfL * δ + sfR * (1.0 - δ)

    # calculate fluxes
    @. ff = (dt * u * f - 0.5 * dt^2 * u^2 * sf) * len
    flux_conserve!(fw, ff, u, v, w, ω)
    #=fw[1] = dt * sum(ω .* u .* f) - 0.5 * dt^2 * sum(ω .* u .^ 2 .* sf)
    fw[2] = dt * sum(ω .* u .^ 2 .* f) - 0.5 * dt^2 * sum(ω .* u .^ 3 .* sf)
    fw[3] = dt * sum(ω .* u .* v .* f) - 0.5 * dt^2 * sum(ω .* u .^ 2 .* v .* sf)
    fw[4] = dt * sum(ω .* u .* w .* f) - 0.5 * dt^2 * sum(ω .* u .^ 2 .* w .* sf)
    fw[5] =
        dt * 0.5 * sum(ω .* u .* (u .^ 2 .+ v .^ 2 .+ w .^ 2) .* f) -
        0.5 * dt^2 * 0.5 * sum(ω .* u .^ 2 .* (u .^ 2 .+ v .^ 2 .+ w .^ 2) .* sf)
    @. fw *= len=#

    return nothing

end

"""
$(SIGNATURES)

4F1V
"""
function flux_kfvs!(
    fw::AV{<:FN},
    fh0::Y,
    fh1::Y,
    fh2::Y,
    fh3::Y,
    h0L::Z,
    h1L::Z,
    h2L::Z,
    h3L::Z,
    h0R::Z,
    h1R::Z,
    h2R::Z,
    h3R::Z,
    u::A,
    ω::A,
    dt,
    sh0L = zero(h0L)::Z,
    sh1L = zero(h1L)::Z,
    sh2L = zero(h2L)::Z,
    sh3L = zero(h3L)::Z,
    sh0R = zero(h0R)::Z,
    sh1R = zero(h1R)::Z,
    sh2R = zero(h2R)::Z,
    sh3R = zero(h3R)::Z,
) where {Y<:AV{<:FN},Z<:AV{<:FN},A<:AV{<:FN}}

    # upwind reconstruction
    δ = heaviside.(u)

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

    sh0 = @. sh0L * δ + sh0R * (1.0 - δ)
    sh1 = @. sh1L * δ + sh1R * (1.0 - δ)
    sh2 = @. sh2L * δ + sh2R * (1.0 - δ)
    sh3 = @. sh3L * δ + sh3R * (1.0 - δ)

    # calculate fluxes
    fw[1] = dt * sum(ω .* u .* h0) - 0.5 * dt^2 * sum(ω .* u .^ 2 .* sh0)
    fw[2] = dt * sum(ω .* u .^ 2 .* h0) - 0.5 * dt^2 * sum(ω .* u .^ 3 .* sh0)
    fw[3] = dt * sum(ω .* u .* h1) - 0.5 * dt^2 * sum(ω .* u .^ 2 .* sh1)
    fw[4] = dt * sum(ω .* u .* h2) - 0.5 * dt^2 * sum(ω .* u .^ 2 .* sh2)
    fw[5] =
        dt * 0.5 * (sum(ω .* u .^ 3 .* h0) + sum(ω .* u .* h3)) -
        0.5 * dt^2 * 0.5 * (sum(ω .* u .^ 4 .* sh0) + sum(ω .* u .^ 2 .* sh3))

    @. fh0 = dt * u * h0 - 0.5 * dt^2 * u^2 * sh0
    @. fh1 = dt * u * h1 - 0.5 * dt^2 * u^2 * sh1
    @. fh2 = dt * u * h2 - 0.5 * dt^2 * u^2 * sh2
    @. fh3 = dt * u * h3 - 0.5 * dt^2 * u^2 * sh3

    return nothing

end

"""
$(SIGNATURES)

Mixture
"""
function flux_kfvs!(
    fw::AM{<:FN},
    fh0::Y,
    fh1::Y,
    fh2::Y,
    fh3::Y,
    h0L::Z,
    h1L::Z,
    h2L::Z,
    h3L::Z,
    h0R::Z,
    h1R::Z,
    h2R::Z,
    h3R::Z,
    u::A,
    ω::A,
    dt,
    sh0L = zero(h0L)::Z,
    sh1L = zero(h1L)::Z,
    sh2L = zero(h2L)::Z,
    sh3L = zero(h3L)::Z,
    sh0R = zero(h0R)::Z,
    sh1R = zero(h1R)::Z,
    sh2R = zero(h2R)::Z,
    sh3R = zero(h3R)::Z,
) where {Y<:AM{<:FN},Z<:AM{<:FN},A<:AM{<:FN}}

    for j in axes(fw, 2)
        _fw = @view fw[:, j]
        _fh0 = @view fh0[:, j]
        _fh1 = @view fh1[:, j]
        _fh2 = @view fh2[:, j]
        _fh3 = @view fh3[:, j]

        flux_kfvs!(
            _fw,
            _fh0,
            _fh1,
            _fh2,
            _fh3,
            h0L[:, j],
            h1L[:, j],
            h2L[:, j],
            h3L[:, j],
            h0R[:, j],
            h1R[:, j],
            h2R[:, j],
            h3R[:, j],
            u[:, j],
            ω[:, j],
            dt,
            sh0L[:, j],
            sh1L[:, j],
            sh2L[:, j],
            sh3L[:, j],
            sh0R[:, j],
            sh1R[:, j],
            sh2R[:, j],
            sh3R[:, j],
        )
    end

    return nothing

end

"""
$(SIGNATURES)

1F2V
"""
function flux_kfvs!(
    fw::AV{<:FN},
    ff::Union{AV{<:FN},AM{<:FN}},
    fL::Z,
    fR::Z,
    u::A,
    v::A,
    ω::A,
    dt,
    len,
    sfL = zero(fL)::Z,
    sfR = zero(fR)::Z,
) where {Z<:Union{AV{<:FN},AM{<:FN}},A<:Union{AV{<:FN},AM{<:FN}}}

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

    f = @. fL * δ + fR * (1.0 - δ)
    sf = @. sfL * δ + sfR * (1.0 - δ)

    #--- calculate fluxes ---#
    @. ff = (dt * u * f - 0.5 * dt^2 * u^2 * sf) * len
    flux_conserve!(fw, ff, u, v, ω)
    #=fw[1] = dt * sum(ω .* u .* f) - 0.5 * dt^2 * sum(ω .* u .^ 2 .* sf)
    fw[2] = dt * sum(ω .* u .^ 2 .* f) - 0.5 * dt^2 * sum(ω .* u .^ 3 .* sf)
    fw[3] = dt * sum(ω .* v .* u .* f) - 0.5 * dt^2 * sum(ω .* v .* u .^ 2 .* sf)
    fw[4] =
        dt * 0.5 * sum(ω .* u .* (u .^ 2 .+ v .^ 2) .* f) -
        0.5 * dt^2 * 0.5 * sum(ω .* u .^ 2 .* (u .^ 2 .+ v .^ 2) .* sf)
    fw .*= len=#

    return nothing

end

"""
$(SIGNATURES)

2F2V
"""
function flux_kfvs!(
    fw::AV{<:FN},
    fh::Y,
    fb::Y,
    hL::Z,
    bL::Z,
    hR::Z,
    bR::Z,
    u::A,
    v::A,
    ω::A,
    dt,
    len,
    shL = zero(hL)::Z,
    sbL = zero(bL)::Z,
    shR = zero(hR)::Z,
    sbR = zero(bR)::Z,
) where {
    Y<:Union{AV{<:FN},AM{<:FN}},
    Z<:Union{AV{<:FN},AM{<:FN}},
    A<:Union{AV{<:FN},AM{<:FN}},
}

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

    h = @. hL * δ + hR * (1.0 - δ)
    b = @. bL * δ + bR * (1.0 - δ)
    sh = @. shL * δ + shR * (1.0 - δ)
    sb = @. sbL * δ + sbR * (1.0 - δ)

    #--- calculate fluxes ---#
    @. fh = (dt * u * h - 0.5 * dt^2 * u^2 * sh) * len
    @. fb = (dt * u * b - 0.5 * dt^2 * u^2 * sb) * len
    flux_conserve!(fw, fh, fb, u, v, ω)
    #=fw[1] = dt * sum(ω .* u .* h) - 0.5 * dt^2 * sum(ω .* u .^ 2 .* sh)
    fw[2] = dt * sum(ω .* u .^ 2 .* h) - 0.5 * dt^2 * sum(ω .* u .^ 3 .* sh)
    fw[3] = dt * sum(ω .* v .* u .* h) - 0.5 * dt^2 * sum(ω .* v .* u .^ 2 .* sh)
    fw[4] =
        dt * 0.5 * (sum(ω .* u .* (u .^ 2 .+ v .^ 2) .* h) + sum(ω .* u .* b)) -
        0.5 *
        dt^2 *
        0.5 *
        (sum(ω .* u .^ 2 .* (u .^ 2 .+ v .^ 2) .* sh) + sum(ω .* u .^ 2 .* sb))
    fw .*= len=#

    return nothing

end

"""
$(SIGNATURES)

3F2V
"""
function flux_kfvs!(
    fw::AV{<:FN},
    fh0::Y,
    fh1::Y,
    fh2::Y,
    h0L::Z,
    h1L::Z,
    h2L::Z,
    h0R::Z,
    h1R::Z,
    h2R::Z,
    u::A,
    v::A,
    ω::A,
    dt,
    len,
    sh0L = zero(h0L)::Z,
    sh1L = zero(h1L)::Z,
    sh2L = zero(h2L)::Z,
    sh0R = zero(h0R)::Z,
    sh1R = zero(h1R)::Z,
    sh2R = zero(h2R)::Z,
) where {Y<:AM{<:FN},Z<:AM{<:FN},A<:AM{<:FN}}

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

    h0 = @. h0L * δ + h0R * (1.0 - δ)
    h1 = @. h1L * δ + h1R * (1.0 - δ)
    h2 = @. h2L * δ + h2R * (1.0 - δ)
    sh0 = @. sh0L * δ + sh0R * (1.0 - δ)
    sh1 = @. sh1L * δ + sh1R * (1.0 - δ)
    sh2 = @. sh2L * δ + sh2R * (1.0 - δ)

    #--- calculate fluxes ---#
    fw[1] = dt * sum(ω .* u .* h0) - 0.5 * dt^2 * sum(ω .* u .^ 2 .* sh0)
    fw[2] = dt * sum(ω .* u .^ 2 .* h0) - 0.5 * dt^2 * sum(ω .* u .^ 3 .* sh0)
    fw[3] = dt * sum(ω .* u .* v .* h0) - 0.5 * dt^2 * sum(ω .* u .^ 2 .* v .* sh0)
    fw[4] = dt * sum(ω .* u .* h1) - 0.5 * dt^2 * sum(ω .* u .^ 2 .* sh1)
    fw[5] =
        dt * 0.5 * (sum(ω .* u .* (u .^ 2 .+ v .^ 2) .* h0) + sum(ω .* u .* h2)) -
        0.5 *
        dt^2 *
        0.5 *
        (sum(ω .* u .^ 2 .* (u .^ 2 .+ v .^ 2) .* sh0) + sum(ω .* u .^ 2 .* sh2))

    fw .*= len
    @. fh0 = (dt * u * h0 - 0.5 * dt^2 * u^2 * sh0) * len
    @. fh1 = (dt * u * h1 - 0.5 * dt^2 * u^2 * sh1) * len
    @. fh2 = (dt * u * h2 - 0.5 * dt^2 * u^2 * sh2) * len

    return nothing

end

"""
$(SIGNATURES)

Mixture
"""
function flux_kfvs!(
    fw::AM{<:FN},
    fh0::Y,
    fh1::Y,
    fh2::Y,
    h0L::Z,
    h1L::Z,
    h2L::Z,
    h0R::Z,
    h1R::Z,
    h2R::Z,
    u::A,
    v::A,
    ω::A,
    dt,
    len,
    sh0L = zero(h0L)::Z,
    sh1L = zero(h1L)::Z,
    sh2L = zero(h2L)::Z,
    sh0R = zero(h0R)::Z,
    sh1R = zero(h1R)::Z,
    sh2R = zero(h2R)::Z,
) where {Y<:AA{<:FN,3},Z<:AA{<:FN,3},A<:AA{<:FN,3}}

    #--- reconstruct initial distribution ---#
    δ = heaviside.(u)

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

    sh0 = @. sh0L * δ + sh0R * (1.0 - δ)
    sh1 = @. sh1L * δ + sh1R * (1.0 - δ)
    sh2 = @. sh2L * δ + sh2R * (1.0 - δ)

    for j = 1:2
        fw[1, j] =
            dt * sum(ω[:, :, j] .* u[:, :, j] .* h0[:, :, j]) -
            0.5 * dt^2 * sum(ω[:, :, j] .* u[:, :, j] .^ 2 .* sh0[:, :, j])
        fw[2, j] =
            dt * sum(ω[:, :, j] .* u[:, :, j] .^ 2 .* h0[:, :, j]) -
            0.5 * dt^2 * sum(ω[:, :, j] .* u[:, :, j] .^ 3 .* sh0[:, :, j])
        fw[3, j] =
            dt * sum(ω[:, :, j] .* v[:, :, j] .* u[:, :, j] .* h0[:, :, j]) -
            0.5 * dt^2 * sum(ω[:, :, j] .* u[:, :, j] .^ 2 .* v[:, :, j] .* sh0[:, :, j])
        fw[4, j] =
            dt * sum(ω[:, :, j] .* u[:, :, j] .* h1[:, :, j]) -
            0.5 * dt^2 * sum(ω[:, :, j] .* u[:, :, j] .^ 2 .* sh1[:, :, j])
        fw[5, j] =
            dt *
            0.5 *
            (
                sum(
                    ω[:, :, j] .* u[:, :, j] .* (u[:, :, j] .^ 2 .+ v[:, :, j] .^ 2) .*
                    h0[:, :, j],
                ) + sum(ω[:, :, j] .* u[:, :, j] .* h2[:, :, j])
            ) -
            0.5 *
            dt^2 *
            0.5 *
            (
                sum(
                    ω[:, :, j] .* u[:, :, j] .^ 2 .* (u[:, :, j] .^ 2 .+ v[:, :, j] .^ 2) .*
                    sh0[:, :, j],
                ) + sum(ω[:, :, j] .* u[:, :, j] .^ 2 .* sh2[:, :, j])
            )

        @. fh0[:, :, j] =
            dt * u[:, :, j] * h0[:, :, j] - 0.5 * dt^2 * u[:, :, j]^2 * sh0[:, :, j]
        @. fh1[:, :, j] =
            dt * u[:, :, j] * h1[:, :, j] - 0.5 * dt^2 * u[:, :, j]^2 * sh1[:, :, j]
        @. fh2[:, :, j] =
            dt * u[:, :, j] * h2[:, :, j] - 0.5 * dt^2 * u[:, :, j]^2 * sh2[:, :, j]
    end

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

    return nothing

end


text:
# This file is a part of MGVI.jl, licensed under the MIT License (MIT).

module ModelPolyfit

using Distributions
using Random
using ValueShapes

export model, true_params, starting_point

const _x1_grid = [Float64(i)/10 for i in 1:25]
const _x2_grid = [Float64(i)/10 + 0.1 for i in 1:15]
const _common_grid = sort(vcat(_x1_grid, _x2_grid))

function _mean(x_grid, p)
    p[1]*10 .+ p[2]*40 .* x_grid .+ p[3]*600 .* x_grid.^2 .+ p[4]*80 .* x_grid.^3
end

function model(p)
    dist1 = Product(Normal.(_mean(_x1_grid, p), p[5]^2*60))
    dist2 = Product(Normal.(_mean(_x2_grid, p), p[5]^2*60))
    NamedTupleDist(a=dist1,
                   b=dist2)
end

const true_params =  [
 -0.3
 -1.5
 0.2
 -0.5
 0.3]

const starting_point = [
  0.2
  0.5
  -0.1
  0.3
 -0.6
]

end # module

ext: jl


text:
workspace();
using BoxChopper;

box = [0 1; 0 1; 0 1]
r₀ = [0,0,0]
nout = [-1,-1,1]

nr₀ = nout⋅r₀
vcbits = BoxChopper.calc_vcbits(box, nout, nr₀)
nvc = count_ones(vcbits)  # number of vertices contained

needflip = false
if nvc ≥ 5
    needflip = true
    nout = -nout
    nr₀ = -nr₀
    vcbits = ~vcbits
    nvc = 8-nvc
end

@code_warntype BoxChopper.rvol_gensect(box, nout, nr₀, vcbits)


text:
# # [Parallel processing tutorial](@id tutorial-parallel)
#
# See also: [Overview of parallel processing in Transducers.jl](@ref
# overview-parallel)

using LiterateTest                                                     #src
using Test                                                             #src

# ## Quick examples
# ### Sequential processing
using Transducers

xs = randn(10_000_000)
foldl(+, Map(sin), xs)

# ### Thread-based parallelism

# Just replace `foldl` with `foldxt`, to make use of multiple cores:

foldxt(+, Map(sin), xs)

# (In my laptop (4 core machine) I start seeing some speedup around
# `length(xs) ≥ 100_000` for this transducer and reducing function.)

# ### Process-based parallelism

using Distributed
#md addprocs(4)

foldxd(+, Map(sin), xs)

# (Note: there is likely no speedup for light-weight computation and
# large input data like this, when using `foldxd`.)

# ### Parallel processing with iterator comprehensions

# You can also use parallel processing functions in Transducers.jl
# such as [`foldxt`](@ref), [`foldxd`](@ref), [`tcollect`](@ref),
# [`dcollect`](@ref), [`tcopy`](@ref) and [`dcopy`](@ref) with
# iterator comprehensions:

foldxt(+, (sin(x) for x in xs if abs(x) < 1); basesize = 500_000)

#-

if VERSION >= v"1.3"                                                   #src
    @test begin
        foldxt(+, (x * y for x in 1:3, y in 1:3))
    end == 36
end                                                                    #src

#-

tcollect(sin(x) for x in xs if abs(x) < 1)

#-

using StructArrays: StructVector
table = StructVector(a = [1, 2, 3], b = [5, 6, 7])

@test begin
    tcopy((A = row.a + 1, B = row.b - 1) for row in table if isodd(row.a))
end == StructVector(A = [2, 4], B = [4, 6])

# ## When can I use `foldxt` and `foldxd`?

# ### Requirement 1: Associative reducing step function
#
# Parallel reductions such as [`foldxt`](@ref) and [`foldxd`](@ref)
# requires _associative_ [reducing step function](@ref Glossary).
# Recall that
# [associativity](https://en.wikipedia.org/wiki/Associative_property)
# means that the _grouping_ of evaluations of binary operator does not
# matter:

op = +  # for example
a, b, c = 1, 2, 3  # for example

@assert op(op(a, b), c) == op(a, op(b, c))

# Given this property, computations like `a + b + c + d` can be done
# with different "groupings":
#
# ```
#   a + b + c + d
#
# = ((a + b) + c) + d                +
#                                   / \
#                                  +   d
#                                 / \
#                                +   c           foldl-like grouping
#                               / \
#                              a   b
#
# = (a + b) + (c + d)                +
#                                   / \
#                                  /   \
#                                 /     \        reduce-like grouping
#                                +       +
#                               / \     / \
#                              a   b   c   d
# ```
#
# Notice that, in the last grouping, computation of `c + d` does not
# have to wait for the result of `a + b`.  This is why we need
# associativity for parallel execution.
#
# !!! note
#     Do not confuse associativity with _commutativity_ `op(a, b) =
#     op(b, a)`.  For example, matrix multiplication `*(::Matrix,
#     ::Matrix)` is associative but _not_ commutative in general.
#     However, since `foldxt` only requires associativity, it is valid
#     to use `foldxt(*, xf, matrices)`.

# As reducing function `+` is associative, it can be used with
# `foldxt` (and `foldxd`):

@dedent y1 = begin
    foldxt(+, Map(identity), 1:10; init = 0, basesize = 1)
end

# and the result is the same as the sequential version:

@dedent y2 = begin
    foldl(+, Map(identity), 1:10; init = 0)
end
@test y1 == y2

# Note: `basesize` is for forcing `foldxt` to avoid falling back to
# `foldl` for small length container such as `1:10`.

# On the other hand, binary function `-` is not associative.  Thus,
# `foldxt` cannot be used instead of `foldl` (they produce different
# result):

@dedent y3 = begin
    foldxt(-, Map(identity), 1:10; init = 0, basesize = 1)
end
#-

@dedent y4 = begin
    foldl(+, Map(identity), 1:10; init = 0)
end
@test y3 != y4

# ### Requirement 2: stateless transducers

# Parallel reduction only work with stateless transducers
# [`Map`](@ref), [`Filter`](@ref), [`Cat`](@ref), etc. and you will
# get an error when using stateful transducers such as `Scan` with
# `foldxt` or `foldxd`:

@evaltest_throw "foldxt(+, Scan(+), 1:10; basesize = 1)" begin
    @test occursin(
        "Stateful transducer Scan(+) does not support `combine`",
        sprint(showerror, ans),
    )
end

# Stateful transducers cannot be used with `foldxt` because it is
# impossible to start processing input collection from the middle when
# the transducers need to know all previous elements (= stateful).
#
# [`ScanEmit`](@ref) is a stateful transducer but it is assumed that
# it is used in a context that outputs can be treated as stateless
# (see: [Splitting a string into words and counting them in
# parallel](@ref parallel-word-count)).

# ## [Example: parallel `collect`](@id tutorial-parallel-collect)
#
# !!! note
#     This section explains the _implementation ideas_ of parallel
#     `collect`.  Pre-defined functions such as [`tcopy`](@ref) and
#     [`dcopy`](@ref) should cover many use-cases.
#
# Suppose (pretend) there is a compute-heavy transducer:

xf_compute = opcompose(Filter(!ismissing), Map(x -> x^2))
nothing                                                              # hide

# Transducers.jl supports applying this to an input container and then
# collecting the results into another container.  It can be done
# sequentially ([`collect`](@ref), [`copy`](@ref), etc.) and in
# parallel using threads ([`tcollect`](@ref), [`tcopy`](@ref)) or
# using multiple processes ([`dcollect`](@ref), [`dcopy`](@ref)).  For
# example:

xs = [abs(x) > 1 ? missing : x for x in randn(10_000)]
y1 = collect(xf_compute, xs)
nothing                                                              # hide

# Doing this in parallel is as easy as using `tcollect` or `dcollect`.
# However, it is easy to do this manually, too:

using BangBang: append!!

singleton_vector(x) = [x]
y2 = foldxt(append!!, xs |> xf_compute |> Map(singleton_vector))
@assert y1 == y2

# This code illustrates the common pattern in parallel processing:
#
# 1. Put a result from the transducer in a "singleton solution".
#    Here, it is `[x]`.
#
# 2. Then "merge" the (singleton) solution into the exsiting one.
#    This is done by `append!!` in the above example.

# To illustrate how `foldxt(append!!, xs |> ... |>
# Map(singleton_vector))` works, let's create a reducing function that
# records the arguments and returned values of `append!!`:

chan = Channel(Inf)

function append_and_log!!(a, b)
    # As arguments and output may be mutated later, we use `copy` to
    # record the snapshots of their values at this moment:
    #+
    a0 = copy(a)
    b0 = copy(b)
    c = append!!(a, b)
    put!(chan, (a0, b0) => copy(c))
    return c
end
nothing                                                              # hide

# This function can be used instead of `append!!`.  Let's try simpler
# and shorter example.  This is equivalent to `collect(1:4)`:

foldxt(append_and_log!!, Map(singleton_vector), 1:4; basesize = 1, init = Union{}[])

# (See below for why we are using `init = Union{}[]` here.)
#
# Here is the list of arguments and returned value of `append!!` in
# this reduction:

records = Pair[]
while isready(chan)
    push!(records, take!(chan))
end
@testset "check records" begin
    if VERSION >= v"1.3"
        @test Set(records) == Set([
            (Union{}[], [1]) => [1],
            (Union{}[], [2]) => [2],
            (Union{}[], [3]) => [3],
            (Union{}[], [4]) => [4],
            ([1], [2]) => [1, 2],
            ([3], [4]) => [3, 4],
            ([1, 2], [3, 4]) => [1, 2, 3, 4],
        ])
    end
end
records

# This recorded inputs and outputs of `append!!` show that its "call
# tree" is:
#
# ```
#           [1,2,3,4] <------------- append!!([1,2], [3,4]) == [1,2,3,4]
#          /         \
#     [1,2]           [3,4] <------- append!!([3], [4]) == [3, 4]
#    /     \         /     \
#  [1]     [2]     [3]     [4] <---- append!!([], [4]) == [4]
#  / \     / \     / \     / \
# [] [1]  [] [2]  [] [3]  [] [4]
# ```
#
# Compare this to the example `a + b + c + d` above.

# ### Optimization and generic container handling
#
# Above usage of `foldxt` is not quite efficient as `singleton_vector`
# allocates small objects in the heap.  Thus, it makes sense to use
# immutable objects for the singleton solutions so that Julia compiler
# can eliminate allocation of the intermediate singleton solutions.
# Here, this can be done by simply using `SVector`instead of
# `singleton_vector`:

using StaticArrays: SVector

@evaltest "foldxt(append!!, Map(SVector), 1:4)" begin
    @test ans isa SVector
    @test ans == 1:4
end

# However, notice that the return value is a static vector.  This is
# not ideal when the input collection is large.  The output collection
# type can be specified by `init`.  We can simply use `init =
# Union{}[]` in this case:

@evaltest "foldxt(append!!, Map(SVector), 1:4; init = Union{}[])" begin
    @test ans isa Vector
    @test ans == 1:4
end

# Note that passing `Vector` to `init` of `foldxt` is usually a wrong
# choice as it would mean that the same object is simultaneously
# mutated by different threads.  However, since `Vector{Union{}}`
# cannot have any element (as there is no object of type `Union{}`),
# using `Union{}[]` for `init` is an exception and it is a good way to
# indicate that output vector should use the "smallest" `eltype`
# required.  That is to say, `append!!` widens the vector "just
# enough" to fit the resulting elements.
#
# For generic containers (e.g., various table types), use
# [`BangBang.Empty`](https://juliafolds.github.io/BangBang.jl/dev/#BangBang.NoBang.singletonof-Union{Tuple{T},%20Tuple{Type{T},Any}}%20where%20T)
# as the empty initial value.  This is useful for creating a table
# object such as
# [`DataFrame`](https://github.com/JuliaData/DataFrames.jl) as the
# result of parallel processing:

using BangBang: Empty
using DataFrames: DataFrame

ans = begin
    foldxt(append!!, Map(x -> SVector((a = x,))), 1:4; init = Empty(DataFrame))
end
@testset "foldxt(append!!, ...; init = Empty(DataFrame))" begin
    @test ans isa DataFrame
    @test ans.a == 1:4
end
Text(ans)                                                            # hide

# It is slightly more tricky to make this approach work with other
# table types such as
# [`StructArrays`](https://github.com/JuliaArrays/StructArrays.jl) and
# [`TypedTables`](https://github.com/JuliaData/TypedTables.jl).  Use
# [`tcopy`](@ref) or [`dcopy`](@ref) to work with generic containers.

#src It doesn't work:
#src
#src # Similar approach works with
#src # [`StructArrays.StructVector`](https://github.com/JuliaArrays/StructArrays.jl):
#src
#src using StructArrays: StructVector
#src
#src y6 =                                                                   #src
#src foldxt(append!!, Map(x -> SVector((a = x,))), 1:4; init = Empty(StructVector))
#src @test y6 is StructVector                                               #src
#src @test y6.a == 1:4                                                      #src
#src
#src # or [`TypedTables.Table`](https://github.com/JuliaData/TypedTables.jl):
#src
#src using TypedTables: Table
#src
#src y7 =                                                                   #src
#src foldxt(append!!, Map(x -> SVector((a = x,))), 1:4; init = Empty(Table))
#src @test y7 is Table                                                      #src
#src @test y7.a == 1:4                                                      #src
#src
#src # or any type that can be constructed from a
#src # [table](https://github.com/JuliaData/Tables.jl) and behaves like a
#src # vector (e.g., support `push` and `vcat`).

# ## Example: ad-hoc histogram

# Following example counts number of occurrence of each leading digit
# in a distribution of random numbers.  First, let's create "singleton
# solutions" using transducers:

xs = 1_000_000 * randn(10_000_000)
dicts1 = xs |> Map(abs) |> Filter(x -> x > 1) |> Map() do x
    y = digits(floor(Int, x))[end]
    Dict(y => 1)
end
nothing                                                              # hide

# The singleton solutions can be merged using `mergewith!(+, a, b)`.
# Conveniently, `mergewith!(+)` is the curried form `(args...) ->
# mergewith!(+, args...)`:

using Compat: mergewith!  # not required in Julia >= 1.5
rf! = mergewith!(+)
@test begin
    rf!(Dict(:a => 1, :b => 2), Dict(:b => 3, :c => 4))
end == Dict(:a => 1, :b => 5, :c => 4)

# This is the form of binary function appropriate for `foldl` and
# `foldxt`.
#
# Note that it is OK to use in-place function `mergewith!` here because
# the dictionary passed as `a` is created by `Dict(y => 1)` and not shared by
# anyone.  When there is no such guarantee, passing [`init =
# OnInit(Dict{Int,Int})`](@ref OnInit) is a good option.  Note that
# passing `init = Dict{Int,Int}()` to `foldxt` is not correct as
# multiple tasks would share and try to mutate the same dictionary
# this way.

# Let's try this with parallel `foldxt`:

counts1 = foldxt(mergewith!(+), dicts1)
nothing                                                              # hide

# Compare the result with `foldl`:
counts2 = foldl(mergewith!(+), dicts1)
@assert counts1 == counts2

# Hopefully the result is close to the [Benford's law -
# Wikipedia](https://en.wikipedia.org/wiki/Benford%27s_law):
let n = sum(values(counts1))
    sort!(keys(counts1) .=> values(counts1) ./ n)
end

# Since we are counting only nine elements, it is actually better to
# use fixed-size container such as a tuple in this case:

dicts2 = xs |> Map(abs) |> Filter(x -> x > 1) |> Map() do x
    y = digits(floor(Int, x))[end]
    ntuple(i -> i == y, 9)
end

counts3 = foldxt(dicts2; init=ntuple(_ -> 0, 9)) do a, b
    map(+, a, b)
end
@assert Dict(zip(1:9, counts3)) == counts1

# Note that, as tuples are immutable, it is valid to pass it as `init`
# of `foldxt`.

# ### MicroCollections.jl for efficient singleton solution
#
# When the appropriate "bins" are not known, `mergewith!(+)`-based
# strategy is more appropriate.  However, it is not ideal to allocate
# a small container like `Dict(y => 1)` in the heap for each
# iteration.
# [MicroCollections.jl](https://github.com/JuliaFolds/MicroCollections.jl)
# provides singleton (and empty) containers that are designed for this
# usecase.  The `SingletonDict` is "upcast" to the mutable `Dict` in
# the first invocation when merged with BangBang.jl functions:

using BangBang: mergewith!!
using MicroCollections: SingletonDict

acc1 = mergewith!!(+, SingletonDict(:a => 1), SingletonDict(:b => 1))
@test acc1 isa Dict
Text(acc1)                                                           # hide

# This dictionary is reused in the subsequent iterations:

acc2 = mergewith!!(+, acc1, SingletonDict(:b => 1))
Text(acc2)                                                           # hide
#-

acc3 = mergewith!!(+, acc2, SingletonDict(:c => 1))
Text(acc3)                                                           # hide

# The first result is reused across these iterations (within a single
# thread).

@assert acc1 === acc2 === acc3

# Finally, `Dict`s from different threads are merged using the same
# function `mergewith!!(+)`:

acc4 = Dict(:a => 5, :c => 3)  # from different thread
acc5 = mergewith!!(+, acc3, acc4)
Text(acc5)                                                           # hide

# Thus, `dicts1` can be optimized simply by replacing `Dict(y => 1)`
# with `SingletonDict(y => 1)`:

dicts3 = xs |> Map(abs) |> Filter(x -> x > 1) |> Map() do x
    y = digits(floor(Int, x))[end]
    SingletonDict(y => 1)
end

counts4 = foldxt(mergewith!!(+), dicts3)
@assert counts1 == counts4

# ## Example: early termination

# Find the _first_ element that is multiple of three:

@test begin
    foldxt(ReduceIf(x -> x % 3 == 0), 1:10; init = nothing, basesize = 1) do _, x
        ## # Uncomment for demo:
        ## x == 3 ? sleep(0.1) : @show x  # give other tasks a chance to finish first
        return x
    end
end == 3

# This snippet always returns `3`, even though the reductions for `c =
# 6` and `c = 9` may finish first.


text:
#
# Simple Pendulum Problem
#
using OrdinaryDiffEq, Plots

#Constants
const g = 9.81
L = 1.0

#Initial Conditions
u₀ = [0,π/2]
tspan = (0.0,6.3)

#Define the problem
function simplependulum(du,u,p,t)
    θ  = u[1]
    dθ = u[2]
    du[1] = dθ
    du[2] = -(g/L)*sin(θ)
end

#Pass to solvers
prob = ODEProblem(simplependulum,u₀, tspan)
sol = solve(prob,Tsit5())

#Plot
plot(sol,linewidth=2,title ="Simple Pendulum Problem", xaxis = "Time", yaxis = "Height", label = ["Theta","dTheta"])

max_stars_repo_name: sergev/vak-opensource

max_issues_repo_name: sergev/vak-opensource

max_forks_repo_name: sergev/vak-opensource


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

using Markdown
using InteractiveUtils

# ╔═╡ ff2837d2-4dd4-11eb-2442-db690b83f6ee
using DataStructures

# ╔═╡ cb44de38-4dc9-11eb-0027-eb8cdb77d454
md"""
# Collections

Data structures are containers that organize and group data together in different ways. Collection is a very general structure that store item that doesn't have any order.

- List is an order collection without a fixed length
- Set is an unorder collection with only unique elements

#### Two Relevant Models of Computation

A computational model describes how units of computations, memories, and communications are organized. This help to determine what operations an algorithm is allowed and hence measure the computational complexity (time and space) of the algorithm.

1. [Random Access Machine (RAM)](https://en.wikipedia.org/wiki/Random-access_machine)
RAM is an abstract machine modelled by a big array (in register) with $\Theta(1)$ time to load word from register from one to another, perform operation(+,−,∗,/,&,|,ˆ). Python list is modelled after RAM

2. [Pointer Machine](https://en.wikipedia.org/wiki/Pointer_machine)
dynamically allocated objects(namedtuple) with O(1) fields (contain word or pointer). Objects with O(1) attributes modelled after Pointer Machine (a.next)
"""

# ╔═╡ 6b81197a-4dca-11eb-1fca-93d5b121174c
md"""
## Array modelled after RAM

An array is a list with indexing. Depending on programming language array might have different properties such as having fixed length, containing only same type of element. But common to all array has locations each element allowing direct access to each element through index (location). They are stored in memory in one _contiguous block_ modelled after RAM.

- Indexable O(1) i.e constant time access: `getindex(l, idx)` and `setindex!(l, data, idx)`
- Insertion of element
  - Insertion at the end is easier, might run into problems when array has a set size
  - Insertion in the beginning and middle is difficult as we need to move everything behind the insertion position a step back. O(n) operation
- Deletion
  - Deletion at any point needs to move everything behind the insertion point a step ahead.

"""

# ╔═╡ b50142de-4dd4-11eb-3a00-a541ecda8cc0
fruitslist = ["🍎","🍌","🍍","🍐", "🍑","🍓"]

# ╔═╡ e889cb4e-4dd4-11eb-2584-05626659c3b3
# pop, append

# ╔═╡ f1164634-4dd4-11eb-0d20-773a92d68ead
md"""
## Linked List modelled after Pointer Machine

Adding and removing element from an array is slightly complicated, need to shift the remaining whole array.

Linked list (A list of sequentially linked nodes) allows efficient insertion of nodes to the front of the list

While Array stores a value and has an index. Linked list stores a value and a pointer to the next element. This can allow linked list to be stored in non-contiguous memory. This also makes it efficient to add and remove element.

![LinkedList](https://upload.wikimedia.org/wikipedia/commons/thumb/6/6d/Singly-linked-list.svg/408px-Singly-linked-list.svg.png)

Insertion takes O(1) time.
"""

# ╔═╡ 082bca6a-4dd5-11eb-2a6c-1b272a805e00
begin
	# imutable list

	l1 = nil()  # create a empty linked list
	l1 = cons(1, l1)
	l1 = cons(2, l1)
	l1 = cons(4, l1)
	l1 = cons(3, l1)
	l1 = cons(6, l1)
	l2 = map((x) -> x*2, l1)
end

# ╔═╡ 24f49228-4de0-11eb-207e-3982680d86c8
md"""
### Pascals Triangle

Triangular array of binomial coefficients where its nth row and kth column can be constructed from previous row as follows:

$ \binom{n}{k}  = \binom{n-1}{k-1} + \binom{n-1}{k} $

![Pascal's Triangle](https://upload.wikimedia.org/wikipedia/commons/0/0d/PascalTriangleAnimated2.gif)
"""

# ╔═╡ ebf4807a-4dda-11eb-2531-add87fbfb72a
function Pascaltriangle(n::Int)
	l = nil(Array{Int});  n == 0 && return l
	l = cons([1], l); n==1 && return l
	l = cons([1,1],l);
	for i = 3:n
		lnew = ones(length(first(l))+1)
		for j =2:i-1
			lnew[j] = first(l)[j-1] + first(l)[j]
		end
		l= cons(lnew, l)
	end
	return l
end

# ╔═╡ 628f92c8-4ddc-11eb-19be-4bcde1b5d680
Pascaltriangle(8)

# ╔═╡ 1380e2ee-4dd5-11eb-1c65-776ea288e5a2
md"""
### Doubly Linked List

Allows traversing in both direction.
![Doubly Linked List](https://upload.wikimedia.org/wikipedia/commons/thumb/5/5e/Doubly-linked-list.svg/610px-Doubly-linked-list.svg.png)

The `MutableLinkedList` type implements a doubly linked list with mutable nodes.

- Indexable O(n):
  - `getindex(l, idx[or range])` or `l[idx]`
  - `setindex!(l,data,idx)` or `l[idx]=data` is defined.
- Insertion:
  - in the begining O(1) : `pushfirst!(l, data)`
  - in the end O(1) : `push!(l, data) ` or `append!(l, elts...) ` or `append!(l1,l2)`
  - in the middle O(n):
- Removal:
  - from the begining O(1) : `popfirst!(l)`
  - from the end O(1) : `pop!(l) `
  - from the middle O(n): `delete!(l,idx[or range])`
"""

# ╔═╡ 31faaa82-4dd5-11eb-31f8-e149054b2990
begin
	l = MutableLinkedList{Int64}(collect(1:100)...)    # initialize

	first(l), last(l)
end

# ╔═╡ 4b74d0f0-4dd5-11eb-239e-1d3e2b4f769e
md"""
#### Reversing a linked list

"""

# ╔═╡ 577fb59c-4dd8-11eb-1800-237d644677ae
reverse(l)

# ╔═╡ 52c3b400-4dd5-11eb-0519-03a4b8f97579
md"""
#### Flatten a nested linked list
"""

# ╔═╡ 5a791804-4dd5-11eb-3f7a-013eddac9a12

# ╔═╡ 5ed65876-4dd5-11eb-389f-458a3913abe2
md"""
### Circular Linked List

Last node of the list rather than having a null reference points back first node of the list.

![Circular Linked List](https://upload.wikimedia.org/wikipedia/commons/thumb/d/df/Circularly-linked-list.svg/350px-Circularly-linked-list.svg.png)
"""

# ╔═╡ d45c7102-4de2-11eb-2d7e-9d4e5843c16f
md"""
### Multiply linked list

As extention of Doubly Linked List, where each node contains two or more link fields that allows to connect the same set of data records in different order.

"""

# ╔═╡ Cell order:
# ╟─cb44de38-4dc9-11eb-0027-eb8cdb77d454
# ╟─6b81197a-4dca-11eb-1fca-93d5b121174c
# ╠═b50142de-4dd4-11eb-3a00-a541ecda8cc0
# ╠═e889cb4e-4dd4-11eb-2584-05626659c3b3
# ╟─f1164634-4dd4-11eb-0d20-773a92d68ead
# ╠═ff2837d2-4dd4-11eb-2442-db690b83f6ee
# ╠═082bca6a-4dd5-11eb-2a6c-1b272a805e00
# ╟─24f49228-4de0-11eb-207e-3982680d86c8
# ╠═ebf4807a-4dda-11eb-2531-add87fbfb72a
# ╠═628f92c8-4ddc-11eb-19be-4bcde1b5d680
# ╟─1380e2ee-4dd5-11eb-1c65-776ea288e5a2
# ╠═31faaa82-4dd5-11eb-31f8-e149054b2990
# ╟─4b74d0f0-4dd5-11eb-239e-1d3e2b4f769e
# ╠═577fb59c-4dd8-11eb-1800-237d644677ae
# ╠═52c3b400-4dd5-11eb-0519-03a4b8f97579
# ╠═5a791804-4dd5-11eb-3f7a-013eddac9a12
# ╟─5ed65876-4dd5-11eb-389f-458a3913abe2
# ╟─d45c7102-4de2-11eb-2d7e-9d4e5843c16f


text:
using Genie
using Genie.Router
using HTTP
import Genie.Renderer.Json: json
using Genie.Renderer.Json, Genie.Requests
using Genie.Renderer.Html
using LibPQ, Tables
using DataFrames, PrettyTables
using Genie.Renderer
conn = LibPQ.Connection("dbname=danenfcgd5khab host=ec2-35-153-114-74.compute-1.amazonaws.com port=5432 user=hbwwyuvguzemdw password=514ffbe17667034ddf0db74ae2d5157c2caf0374c5ac127d0ce38a679345786e sslmode=require")

function launchServer(port)

    Genie.config.run_as_server = true
    Genie.config.server_host = "0.0.0.0"
    Genie.config.server_port = port

    println("port set to $(port)")

    route("/") do
        "Hi there! This is server 1"
    end

    route("/jsontest") do
        (:message => "Hi there!this is a json test") |> json
    end

    route("/echo", method = POST) do
        message = jsonpayload()
        (:echo => (message["message"] * " ")^message["repeat"]) |> json
    end

    route("/send") do
        response = HTTP.request(
            "POST",
            "https://julia-apiserver1.herokuapp.com/echo",
            [("Content-Type", "application/json")],
            """{"message":"hello", "repeat":3}""",
        )

        response.body |> String |> json
    end

    route("/create") do
        result = execute(conn, """
        CREATE TABLE customers (
    customer_id bpchar NOT NULL,
    company_name character varying(40) NOT NULL,
    contact_name character varying(30),
    contact_title character varying(30),
    address character varying(60),
    city character varying(15),
    region character varying(15),
    postal_code character varying(10),
    country character varying(15),
    phone character varying(24),
    fax character varying(24)
);
        """)
    end

    route("/delete") do
        result = execute(conn, """
        DROP TABLE IF EXISTS customers;
        """)
    end

    route("/adddata") do
        result = execute(conn, """
        INSERT INTO customers VALUES ('ALFKI', 'Alfreds Futterkiste', 'Maria Anders', 'Sales Representative', 'Obere Str. 57', 'Berlin', NULL, '12209', 'Germany', '030-0074321', '030-0076545');
        INSERT INTO customers VALUES ('ANATR', 'Ana Trujillo Emparedados y helados', 'Ana Trujillo', 'Owner', 'Avda. de la Constitución 2222', 'México D.F.', NULL, '05021', 'Mexico', '(5) 555-4729', '(5) 555-3745');
        INSERT INTO customers VALUES ('ANTON', 'Antonio Moreno Taquería', 'Antonio Moreno', 'Owner', 'Mataderos  2312', 'México D.F.', NULL, '05023', 'Mexico', '(5) 555-3932', NULL);
        INSERT INTO customers VALUES ('AROUT', 'Around the Horn', 'Thomas Hardy', 'Sales Representative', '120 Hanover Sq.', 'London', NULL, 'WA1 1DP', 'UK', '(171) 555-7788', '(171) 555-6750');
        INSERT INTO customers VALUES ('BERGS', 'Berglunds snabbköp', 'Christina Berglund', 'Order Administrator', 'Berguvsvägen  8', 'Luleå', NULL, 'S-958 22', 'Sweden', '0921-12 34 65', '0921-12 34 67');
        """)
    end

    holder = "SELECT * FROM COMPANY;"

    form = """
    <form action="/command" method="POST" enctype="multipart/form-data">
    <textarea id="name" style="resize: none;" name="name" rows="4" cols="50" />$holder</textarea>

    <input type="submit" value="Run Postgres SQL >>" />
    </form>
    """

    restore = """
    DROP TABLE IF EXISTS COMPANY;
    CREATE TABLE COMPANY(
    ID INT PRIMARY KEY     NOT NULL,
    NAME           TEXT    NOT NULL,
    AGE            INT     NOT NULL,
    ADDRESS        CHAR(50),
    SALARY         REAL,
    JOIN_DATE	  DATE
    );
    INSERT INTO COMPANY (ID,NAME,AGE,ADDRESS,SALARY,JOIN_DATE) VALUES (1, 'Paul', 32, 'California', 20000.00,'2001-07-13');
    INSERT INTO COMPANY (ID,NAME,AGE,ADDRESS,JOIN_DATE) VALUES (2, 'Allen', 25, 'Texas', '2007-12-13');
    INSERT INTO COMPANY (ID,NAME,AGE,ADDRESS,SALARY,JOIN_DATE) VALUES (3, 'Teddy', 23, 'Norway', 20000.00, DEFAULT );
    INSERT INTO COMPANY (ID,NAME,AGE,ADDRESS,SALARY,JOIN_DATE) VALUES (4, 'Mark', 25, 'Rich-Mond ', 65000.00, '2007-12-13' ), (5, 'David', 27, 'Texas', 85000.00, '2007-12-13');
    """

    restoreform = """
    <form action="/restore" method="POST" enctype="multipart/form-data">
    <input type="submit" value="Restore Database" />
    </form>
    """

    route("/restore", method = POST) do
        result = execute(conn, """ $restore""")
        redirect(:get_com)

    end

    route("/command",named=:get_com) do
        html("SQL STATEMENT:\n" * form * "\n Restore Database:\n" * restoreform)
    end

    route("/command", method = POST) do
        holder = "$(postpayload(:name, "Anon"))"
        form = """
        <form action="/command" method="POST" enctype="multipart/form-data">
        <textarea id="name" style="resize: none;" name="name" rows="4" cols="50" />$holder</textarea>

        <input type="submit" value="Run Postgres SQL >>" />
        </form>
        """
        result = DataFrame(execute(conn, """ $(postpayload(:name, "Anon")) """))
        data = pretty_table(String, result; backend = Val(:html))
        html("SQL STATEMENT:\n" * form * "\n " * "Output From Query:" * data * "\n Restore Database:\n" * restoreform)
    end

    Genie.AppServer.startup(async = false)
end
conn = LibPQ.Connection("dbname=danenfcgd5khab host=ec2-35-153-114-74.compute-1.amazonaws.com port=5432 user=hbwwyuvguzemdw password=514ffbe17667034ddf0db74ae2d5157c2caf0374c5ac127d0ce38a679345786e sslmode=require")

launchServer(parse(Int, ARGS[1]))
launchServer( 8001)

max_stars_repo_name: ashwani-rathee/Sigil

max_issues_repo_name: ashwani-rathee/Sigil

max_forks_repo_name: ashwani-rathee/Sigil


text:
mutable struct Quad <: AbstractTweenEquation
    name::String
    # A function that accepts a time argument and returns a float
    compute::Function # f(t::Float64) -> Float64
    type::Int64

    function Quad(type::Int64)
        o = new()

        if type == EASE_IN
            o.type = EASE_IN
            o.name = "Quad.easeIn"
            o.compute = _quad_ease_in
        elseif type == EASE_OUT
            o.type = EASE_OUT
            o.name = "Quad.easeOut"
            o.compute = _quad_ease_out
        else
            o.type = EASE_INOUT
            o.name = "Quad.easeInOut"
            o.compute = _quad_ease_inout
        end

        o
    end
end

function _quad_ease_in(t::Float64)
    t * t
end

function _quad_ease_out(t::Float64)
    -t * (t - 2.0)
end

function _quad_ease_inout(t::Float64)
    a = t * 2.0
    if a < 1.0
        return 0.5 * a * a
    end

    -0.5 * ((t - 1.0) * (t - 2.0) - 1.0)
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
module TestMetrics

using Phylo
using Test

@testset "Metrics" begin
    trees = open(parsenexus, Phylo.path("H1N1.trees"))
    tree1, tree2 = trees["TREE1"], trees["TREE2"]
    species = ["H1N1_A_PUERTORICO_8_1934", "H1N1_A_BRAZIL_11_1978"];
    mrca1 = mrca(tree1, species)

    # MRCA should be the root
    @test mrca1 == getnodename(tree1, getroot(tree1))
    nodes = getnodenames(tree1)
    desc1 = getdescendants(tree1, mrca1)

    # MRCA + descendants is whole tree
    @test all(sort(nodes) .== sort(desc1 ∪ [mrca1]))
    leaves = getleafnames(tree1)
    leaves1 = sort(leaves ∩ desc1)

    # All leaves are in descendants
    @test length(leaves1) == length(leaves)

    # All leaves are the same in both trees
    mrca2 = mrca(tree2, reverse(species))
    leaves2 = sort(leaves ∩ getdescendants(tree2, mrca2))
    @test all(leaves1 .== leaves2)

    # MRCA for Sendai strains has 8 descendants, including Fukuoka
    desendais =
        getdescendants(tree1,
                       mrca(tree1,
                                       filter(x -> occursin("SENDAI", x),
                                       leaves)))
    @test length(desendais) == 8
    @test length(filter(x -> occursin("FUKUOKA", x), desendais)) == 1
end

end


text:
module Renderables

import Term:
    split_lines,
    reshape_text,
    ltrim_str,
    fillin,
    join_lines,
    unescape_brackets_with_space,
    TERM_DEBUG_ON

import ..Style: get_style_codes, MarkupStyle, apply_style
import Term: highlight as highlighter
import ..Consoles: console_width
import ..Measures: Measure
import ..Segments: Segment

export AbstractRenderable, Renderable, RenderableText

# ------------------------------- abstract type ------------------------------ #

"""
    AbstractRenderable
"""
abstract type AbstractRenderable end

Measure(renderable::AbstractRenderable) = renderable.measure

info(r::AbstractRenderable)::String =
    "\e[38;5;117m$(typeof(r)) <: AbstractRenderable\e[0m \e[2m(w:$(r.measure.w), h:$(r.measure.h))\e[0m"

"""
    Base.string(r::AbstractRenderable)::String

Creates a string representation of a renderable
"""
Base.string(r::AbstractRenderable)::String = join([seg.text for seg in r.segments], "\n")

"""
    print(io::IO, renderable::AbstractRenderable)

Print a renderable to an IO
"""
function Base.print(io::IO, renderable::AbstractRenderable; highlight = true)
    ren = unescape_brackets_with_space(string(renderable))
    return println(io, ren)
end

"""
    show(io::IO, renderable::AbstractRenderable)

Show a renderable.
"""
Base.show(io::IO, renderable::AbstractRenderable) = print(io, info(renderable))

"""
    show(io::IO, mime::MIME"text/plain", renderable::AbstractRenderable)

Show a renderable and some information about its shape.
"""
function Base.show(io::IO, ::MIME"text/plain", renderable::AbstractRenderable)
    if TERM_DEBUG_ON[]
        println(io, string(renderable))
        println(io, info(renderable))
    else
        print(io, string(renderable))
    end
end

# ------------------------- generic renderable object ------------------------ #

"""
    Renderable

Generic `Renderable` object.
"""

mutable struct Renderable <: AbstractRenderable
    segments::Vector{Segment}
    measure::Measure
end

Base.convert(::Renderable, x) = Renderable(x)

"""
    Renderable(
        str::String; width::Union{Nothing,Int} = nothing
    )

Convenience method to construct a RenderableText
"""
Renderable(str::AbstractString; width::Union{Nothing,Int} = nothing) =
    RenderableText(str; width = width)

Renderable(ren::AbstractRenderable) = ren
Renderable() = Renderable(Vector{Segment}[], Measure(0, 0))

# ---------------------------------------------------------------------------- #
#                                TEXT RENDERABLE                               #
# ---------------------------------------------------------------------------- #

"""
    RenderableText

`Renderable` representing a text.

See also [`Renderable`](@ref), [`TextBox`](@ref)
"""

mutable struct RenderableText <: AbstractRenderable
    segments::Vector
    measure::Measure
    style::Union{Nothing,String}
end

"""
    RenderableText(text::String; width::Union{Nothing, Int, Symbol}=nothing)

Construct a `RenderableText` out of a string.

If a `width` is passed the text is resized to match the width.
"""
function RenderableText(
    text::AbstractString;
    style::Union{Nothing,String} = nothing,
    width::Union{Nothing,Int} = nothing,
)
    text = apply_style(text)

    # reshape text
    if !isnothing(width)
        width = min(console_width(stdout) - 1, width)
        text = reshape_text(text, width)
    end
    text = fillin(text)

    # create renderable
    segments = if isnothing(style)
        Segment.(
            isnothing(width) ? split_lines(text) :
            map(ln -> rpad(ln, width - textwidth(ln) + 1), split_lines(text))
        )
    else
        style_init, style_finish = get_style_codes(MarkupStyle(style))
        map(ln -> Segment(style_init * ln * style_finish), split_lines(text))
    end

    # @info "a" Measure(segments) segments[1]

    return RenderableText(segments, Measure(segments), style)
end

"""
    RenderableText(rt::RenderableText; width::Union{Nothing,Int} = nothing)

Construct a RenderableText by possibly re-shaping a RenderableText
"""
function RenderableText(
    rt::RenderableText;
    style::Union{Nothing,String} = nothing,
    width::Union{Nothing,Int} = nothing,
)
    return if rt.style == style && rt.measure.w == width
        rt
    else
        text = join_lines([seg.text for seg in rt.segments])
        RenderableText(text; style = style, width = width)
    end
end

# -------------------------------- union type -------------------------------- #
RenderablesUnion = Union{AbstractString,AbstractRenderable}

end


text:
#!/usr/bin/env julia

using Formatting: printfmt

const DATA_DIR = "../08_data"
const INPUT_FILE = "input.txt"
# const INPUT_FILE = "input_test.txt"

macro assert(bool_expr)
    message = string("Assertion: ", bool_expr, " failed")
    :($(esc(bool_expr)) || error($message))
end

@enum CommandType nop acc jmp

mutable struct Command
    cmd_type::CommandType
    cmd_val::Int32
    #--------------------
    function Command(line)
        (a, b) = split(strip(line), " ")
        return new(getproperty(Main, Symbol(a)), parse(Int32, b))
    end
end

mutable struct VM
    vm_code::Array{Command}
    vm_pc::Int32
    vm_acc::Int32
    vm_executed::Set{Int32}
    #--------------------
    function VM(code)
        return new(code, 1, 0, Set{Int32}())
    end
end

# ----------------------------------------

function is_jmp(x)
    return x[2].cmd_type == jmp
end

function fst(pair)
    return pair[1]
end

function vm_exec_cmd(vm::VM)
    push!(vm.vm_executed, vm.vm_pc)
    cmd = vm.vm_code[vm.vm_pc]
    ct = cmd.cmd_type
    if ct == nop
        # println("About to execute noop")
        vm.vm_pc += 1
    elseif ct == acc
        # println("About to execute acc ", cmd.cmd_val)
        vm.vm_acc += cmd.cmd_val
        vm.vm_pc += 1
    elseif ct == jmp
        # println("About to execute jmp ", cmd.cmd_val)
        vm.vm_pc += cmd.cmd_val
    else
        error("Internal error: Unrecognized command type")
    end
end

function vm_run(vm::VM)
    vm.vm_pc = 1
    while true
        if vm.vm_pc in vm.vm_executed
            return ("loop", vm.vm_acc)
        elseif vm.vm_pc == length(vm.vm_code) + 1
            return ("completed", vm.vm_acc)
        else
            vm_exec_cmd(vm)
        end
    end
end

function vm_run_until_loop(vm::VM)
    vm.vm_pc = 1
    while !(vm.vm_pc in vm.vm_executed)
        vm_exec_cmd(vm)
    end
    println("PC=", vm.vm_pc, " for 2nd time. acc=", vm.vm_acc)
end

# ----------------------------------------

function part1()
    println("========== PART 1 ==========")
    lines = readlines("$DATA_DIR/$INPUT_FILE")
    code = map(Command, lines)
    vm = VM(code)
    vm_run_until_loop(vm)
end

function part2()
    println("========== PART 2 ==========")
    lines = readlines("$DATA_DIR/$INPUT_FILE")
    println("# Number of lines=", length(lines))
    code = map(Command, lines)
    jmp_inds = [ecmd[1] for ecmd in enumerate(code) if ecmd[2].cmd_type == jmp]
    for cmd_ind in 1:1:length(code)
        if code[cmd_ind].cmd_type == jmp
            # println("Swapping cmd#$cmd_ind from jmp to nop")
            code_fixed = deepcopy(code)
            code_fixed[cmd_ind].cmd_type = nop
        elseif code[cmd_ind].cmd_type == nop
            # println("Swapping cmd#$cmd_ind from nop to jmp")
            code_fixed = deepcopy(code)
            code_fixed[cmd_ind].cmd_type = jmp
        else
            continue
        end

        vm = VM(code_fixed)
        (exit_code, exit_value) = vm_run(vm)
        if exit_code == "completed"
            println("# Swapped cmd #$(cmd_ind)")
            printfmt(
                "(exit_code, exit_value)=($exit_code, $exit_value)\n"
                , exit_code, exit_value
            )
            exit(0)
        end
    end
    println("Did not find cmd swap that causes successful program termination.\n")
end

# ----------------------------------------

function main()
    part1()
    part2()
end

main()

max_stars_repo_name: jaycoskey/AdventOfCode

max_issues_repo_name: jaycoskey/AdventOfCode

max_forks_repo_name: jaycoskey/AdventOfCode


text:
function mapsURL(path::String; type="terrain", token::String=ENV["GOOGLE_MAPS_API"],
                 size=800, scale=2, MarkersStart=(0.0, 0.0), MarkersEnd=(0.0, 0.0))
    mapURL = "https://maps.googleapis.com/maps/api/staticmap?maptype="
    mapType = type
    mapPath = "&path=enc:" * path * "&"
    mapKey = "key=$token" * "&"
    mapSize = "size=$size" * "x$size" * "&"
    mapScale = "scale=$scale"

    mapMarkersStart = ""
    if MarkersStart != (0.0, 0.0)
        mapMarkersStart = "&markers=color:yellow|label:S|$(MarkersStart[1]),$(MarkersStart[2])"
    end

    mapMarkersEnd = ""
    if MarkersEnd != (0.0, 0.0)
        mapMarkersEnd = "&markers=color:green|label:F|$(MarkersEnd[1]),$(MarkersEnd[2])"
    end

    return mapURL * mapType * mapPath * mapKey * mapSize * mapScale *
           mapMarkersStart * mapMarkersEnd
end

function getMapImage(URL::String; pathFig="/tmp/polyline.png")
    download(URL, pathFig)
end


text:
using Test
using Random
using Statistics
using BenchmarkTools
using PolynomialRoots

using QuarticSolver

#==
These tests simply compare PolynomialRoots using BigFloats to the QuarticSolver
using Float64 (and the Complex versions of both).
==#

# Generate random coefficients (PolynomialRoots requires coeffs in reverse order)
rng = MersenneTwister(235)

rcoeffs = rand(rng, Float64, 5, 1000)
rncoeffs = randn(rng, Float64, 5, 1000)
ccoeffs = rand(rng, ComplexF64, 5, 1000)
cncoeffs = randn(rng, ComplexF64, 5, 1000)
largeccoeffs = rand(rng, ComplexF64, 5, 1000)*100

revrcoeffs = reverse(rcoeffs, dims=1)
revrncoeffs = reverse(rncoeffs, dims=1)
revccoeffs = reverse(ccoeffs, dims=1)
revcncoeffs = reverse(cncoeffs, dims=1)
largerevccoeffs = reverse(largeccoeffs, dims=1)

@testset "real, positive coefficient samples" begin
    pr_rresults = [sort(roots(collect(big.(coeffs))),
                        by=x->(round(real(x), digits=6), round(imag(x), digits=6))) for coeffs in eachcol(revrcoeffs)]
    qs_rresults = [sort([QuarticSolver.solvequartic(coeffs...)...],
                        by=x->(round(real(x), digits=6), round(imag(x), digits=6))) for coeffs in eachcol(rcoeffs)]

    rmatches = isapprox.(pr_rresults, qs_rresults, atol=1e-8)
    rfraction = sum(rmatches)/length(rmatches)
    rdiff = pr_rresults - qs_rresults
    avgdiff = mean(rdiff)
    avgdiffreal = mean(real.(rdiff))
    avgdiffimag = mean(imag.(rdiff))
    mindiffreal = minimum(real.(rdiff))
    mindiffimag = minimum(imag.(rdiff))
    maxdiffreal = maximum(real.(rdiff))
    maxdiffimag = maximum(imag.(rdiff))
    maxabsdiffreal = maximum(maximum(abs.(x) for x in real.(rdiff)))
    maxabsdiffimag = maximum(maximum(abs.(x) for x in imag.(rdiff)))

    argmax(real.(rdiff))
    argmax(imag.(rdiff))
    argmin(real.(rdiff))
    argmin(imag.(rdiff))

    @test rfraction == 1
    @test all(avgdiffreal .< 1e-10)
    @test all(avgdiffimag .< 1e-10)
    @test maxabsdiffreal < 1e-8
    @test maxabsdiffimag < 1e-8
end

@testset "real, positive and negative coefficient samples" begin
    pr_rnresults = [sort(roots(collect(big.(coeffs))),
                         by=x->(round(real(x), digits=6), round(imag(x), digits=6))) for coeffs in eachcol(revrncoeffs)]
    qs_rnresults = [sort([QuarticSolver.solvequartic(coeffs...)...],
                         by=x->(round(real(x), digits=6), round(imag(x), digits=6))) for coeffs in eachcol(rncoeffs)]

    rnmatches = isapprox.(pr_rnresults, qs_rnresults, atol=1e-8)
    rnfraction = sum(rnmatches)/length(rnmatches)
    rndiff = pr_rnresults - qs_rnresults
    avgdiff = mean(rndiff)
    avgdiffreal = mean(real.(rndiff))
    avgdiffimag = mean(imag.(rndiff))
    mindiffreal = minimum(real.(rndiff))
    mindiffimag = minimum(imag.(rndiff))
    maxdiffreal = maximum(real.(rndiff))
    maxdiffimag = maximum(imag.(rndiff))
    maxabsdiffreal = maximum(maximum(abs.(x) for x in real.(rndiff)))
    maxabsdiffimag = maximum(maximum(abs.(x) for x in imag.(rndiff)))

    argmax(imag.(rndiff))
    argmin(imag.(rndiff))

    @test rnfraction == 1
    @test all(avgdiffreal .< 1e-10)
    @test all(avgdiffimag .< 1e-10)
    @test maxabsdiffreal < 1e-8
    @test maxabsdiffimag < 1e-8
end

@testset "complex, positive coefficient samples" begin
    pr_cresults = [sort(roots(collect(big.(coeffs))),
                        by=x->(round(real(x), digits=6), round(imag(x), digits=6))) for coeffs in eachcol(revccoeffs)]
    qs_cresults = [sort([QuarticSolver.solvequartic(coeffs...)...],
                        by=x->(round(real(x), digits=6), round(imag(x), digits=6))) for coeffs in eachcol(ccoeffs)]

    cmatches = isapprox.(pr_cresults, qs_cresults, atol=1e-8)
    cfraction = sum(cmatches)/length(cmatches)
    cdiff = pr_cresults - qs_cresults
    avgdiff = mean(cdiff)
    avgdiffreal = mean(real.(cdiff))
    avgdiffimag = mean(imag.(cdiff))
    mindiffreal = minimum(real.(cdiff))
    mindiffimag = minimum(imag.(cdiff))
    maxdiffreal = maximum(real.(cdiff))
    maxdiffimag = maximum(imag.(cdiff))
    maxabsdiffreal = maximum(maximum(abs.(x) for x in real.(cdiff)))
    maxabsdiffimag = maximum(maximum(abs.(x) for x in imag.(cdiff)))

    argmax(imag.(cdiff))
    argmin(imag.(cdiff))

    @test cfraction == 1
    @test all(avgdiffreal .< 1e-10)
    @test all(avgdiffimag .< 1e-10)
    @test maxabsdiffreal < 1e-8
    @test maxabsdiffimag < 1e-8
end

@testset "complex, positive and negative coefficient samples" begin
    pr_cnresults = [sort(roots(collect(big.(coeffs))),
                         by=x->(round(real(x), digits=6), round(imag(x), digits=6))) for coeffs in eachcol(revcncoeffs)]
    qs_cnresults = [sort([QuarticSolver.solvequartic(coeffs...)...],
                         by=x->(round(real(x), digits=6), round(imag(x), digits=6))) for coeffs in eachcol(cncoeffs)]

    cnmatches = isapprox.(pr_cnresults, qs_cnresults, atol=1e-8)
    cnfraction = sum(cnmatches)/length(cnmatches)
    cndiff = pr_cnresults - qs_cnresults
    avgdiff = mean(cndiff)
    avgdiffreal = mean(real.(cndiff))
    avgdiffimag = mean(imag.(cndiff))
    mindiffreal = minimum(real.(cndiff))
    mindiffimag = minimum(imag.(cndiff))
    maxdiffreal = maximum(real.(cndiff))
    maxdiffimag = maximum(imag.(cndiff))
    maxabsdiffreal = maximum(maximum(abs.(x) for x in real.(cndiff)))
    maxabsdiffimag = maximum(maximum(abs.(x) for x in imag.(cndiff)))

    argmax(imag.(cndiff))
    argmin(imag.(cndiff))

    @test cnfraction == 1
    @test all(avgdiffreal .< 1e-10)
    @test all(avgdiffimag .< 1e-10)
    @test maxabsdiffreal < 1e-8
    @test maxabsdiffimag < 1e-8
end

@testset "complex, positive coefficients larger than 1" begin
    pr_lcresults = [sort(roots(collect(big.(coeffs))),
                         by=x->(round(real(x), digits=6), round(imag(x), digits=6))) for coeffs in eachcol(largerevccoeffs)]
    qs_lcresults = [sort([QuarticSolver.solvequartic(coeffs...)...],
                         by=x->(round(real(x), digits=6), round(imag(x), digits=6))) for coeffs in eachcol(largeccoeffs)]

    lcmatches = isapprox.(pr_lcresults, qs_lcresults, atol=1e-8)
    lcfraction = sum(lcmatches)/length(lcmatches)
    lcdiff = pr_lcresults - qs_lcresults
    avgdiff = mean(lcdiff)
    avgdiffreal = mean(real.(lcdiff))
    avgdiffimag = mean(imag.(lcdiff))
    mindiffreal = minimum(real.(lcdiff))
    mindiffimag = minimum(imag.(lcdiff))
    maxdiffreal = maximum(real.(lcdiff))
    maxdiffimag = maximum(imag.(lcdiff))
    maxabsdiffreal = maximum(maximum(abs.(x) for x in real.(lcdiff)))
    maxabsdiffimag = maximum(maximum(abs.(x) for x in imag.(lcdiff)))

    argmax(imag.(lcdiff))
    argmin(imag.(lcdiff))

    @test lcfraction == 1
    @test all(avgdiffreal .< 1e-10)
    @test all(avgdiffimag .< 1e-10)
    @test maxabsdiffreal < 1e-8
    @test maxabsdiffimag < 1e-8
end

singleccoeffs = ccoeffs[:, 13]
singlerevccoeffs = reverse(singleccoeffs)
@benchmark roots($singlerevccoeffs)
@benchmark QuarticSolver.solvequartic($singleccoeffs[1],
                                      $singleccoeffs[2],
                                      $singleccoeffs[3],
                                      $singleccoeffs[4],
                                      $singleccoeffs[5])

max_stars_repo_name: EP-Guy/QuarticSolver

max_issues_repo_name: EP-Guy/QuarticSolver

max_forks_repo_name: EP-Guy/QuarticSolver


text:
export get_matrix_params, get_conv_params, get_example_network_params

"""
$(SIGNATURES)

Helper function to import the parameters for a layer carrying out matrix multiplication
    (e.g. fully connected layer / softmax layer) from `param_dict` as a
    [`Linear`](@ref) object.

The default format for the key is `'layer_name/weight'` and `'layer_name/bias'`;
    you can customize this by passing in the named arguments `matrix_name` and `bias_name`
    respectively. The expected parameter names will then be `'layer_name/matrix_name'`
    and `'layer_name/bias_name'`

# Arguments
* `param_dict::Dict{String}`: Dictionary mapping parameter names to array of weights
    / biases.
* `layer_name::String`: Identifies parameter in dictionary.
* `expected_size::NTuple{2, Int}`: Tuple of length 2 corresponding to the expected size
   of the weights of the layer.

"""
function get_matrix_params(
    param_dict::Dict{String},
    layer_name::String,
    expected_size::NTuple{2,Int};
    matrix_name::String = "weight",
    bias_name::String = "bias",
)::Linear

    params = Linear(
        param_dict["$layer_name/$matrix_name"],
        dropdims(param_dict["$layer_name/$bias_name"], dims = 1),
    )

    check_size(params, expected_size)

    return params
end

"""
$(SIGNATURES)

Helper function to import the parameters for a convolution layer from `param_dict` as a
    [`Conv2d`](@ref) object.

The default format for the key is `'layer_name/weight'` and `'layer_name/bias'`;
    you can customize this by passing in the named arguments `matrix_name` and `bias_name`
    respectively. The expected parameter names will then be `'layer_name/matrix_name'`
    and `'layer_name/bias_name'`

# Arguments
* `param_dict::Dict{String}`: Dictionary mapping parameter names to array of weights
    / biases.
* `layer_name::String`: Identifies parameter in dictionary.
* `expected_size::NTuple{4, Int}`: Tuple of length 4 corresponding to the expected size
    of the weights of the layer.

"""
function get_conv_params(
    param_dict::Dict{String},
    layer_name::String,
    expected_size::NTuple{4,Int};
    matrix_name::String = "weight",
    bias_name::String = "bias",
    expected_stride::Integer = 1,
    padding::Padding = SamePadding(),
)::Conv2d

    params = Conv2d(
        param_dict["$layer_name/$matrix_name"],
        dropdims(param_dict["$layer_name/$bias_name"], dims = 1),
        expected_stride,
        padding,
    )

    check_size(params, expected_size)

    return params
end


text:
include("modules/Processing.jl")

using BenchmarkTools
using .Processing

checkHeight = 100
sectors = [22, 12, 21, 13] # Top, Right, Bottom, Left
serialPort = "/dev/ttyACM0"
monitor = 2
frameLimiter = 60

config = init_processing(monitor, checkHeight, sectors, serialPort)
start_processing(config, frameLimiter)


text:
py"""
import torch.nn as nn
import torch
import torch.nn.functional as F
import numpy as np
from torch.nn import init

class Down_Conv(nn.Module):

    def __init__(self, in_channels, out_channels, mid_channels=None, mode = 'double_conv'):
        super().__init__()
        assert mode in {'double_conv', 'single_conv'}
        if not mid_channels:
            mid_channels = out_channels
        if mode == 'double_conv':
          self.conv = nn.Sequential(
              nn.Conv2d(in_channels, mid_channels, kernel_size=3, padding=1, bias = False),
              nn.BatchNorm2d(mid_channels),
              nn.LeakyReLU(0.1, inplace=True),
              nn.Conv2d(mid_channels, out_channels, kernel_size=3, padding=1, bias = False),
              nn.BatchNorm2d(out_channels),
              nn.LeakyReLU(0.1, inplace=True)
          )
        else:
          self.conv = nn.Sequential(
              nn.Conv2d(in_channels, out_channels, kernel_size=3, padding=1, bias = False),
              nn.BatchNorm2d(mid_channels),
              nn.LeakyReLU(0.1, inplace=True)
          )

    def forward(self, x):
        return self.conv(x)

class Up_Conv(nn.Module):

    def __init__(self, in_channels, out_channels, mid_channels=None, mode = 'double_conv'):
        super().__init__()
        assert mode in {'double_conv', 'single_conv'}
        if not mid_channels:
            mid_channels = out_channels
        if mode == 'double_conv':
          self.conv = nn.Sequential(
              nn.Conv2d(in_channels, mid_channels, kernel_size=3, padding=1, bias = False),
              nn.BatchNorm2d(mid_channels),
              nn.ReLU(inplace=True),
              nn.Conv2d(mid_channels, out_channels, kernel_size=3, padding=1, bias = False),
              nn.BatchNorm2d(out_channels),
              nn.ReLU(inplace=True)
          )
        else:
          self.conv = nn.Sequential(
              nn.Conv2d(in_channels, out_channels, kernel_size=3, padding=1, bias = False),
              nn.BatchNorm2d(out_channels),
              nn.ReLU(inplace=True)
          )

    def forward(self, x):
        return self.conv(x)

class Down(nn.Module):

    def __init__(self, in_channels, out_channels, mode = 'double_conv'):
        super().__init__()
        self.maxpool_conv = nn.Sequential(
            nn.MaxPool2d(2),
            Down_Conv(in_channels, out_channels, mode = mode)
        )

    def forward(self, x):
        return self.maxpool_conv(x)

class Up(nn.Module):

    def __init__(self, in_channels, out_channels, nearest=True, mode = 'double_conv'):
        super().__init__()

        # if nearest, use the normal convolutions to reduce the number of channels
        if nearest:
            self.up = nn.Upsample(scale_factor=2, mode='nearest')
            self.conv = Up_Conv(in_channels, out_channels // 2, in_channels // 2, mode = mode)
        else:
            self.up = nn.ConvTranspose2d(in_channels , in_channels // 2, kernel_size=2, stride=2)
            self.conv = Up_Conv(in_channels, out_channels, mode = mode)

    def forward(self, x1, x2):
        x1 = self.up(x1)
        # input is CHW
        diffX = x2.size()[2] - x1.size()[2]
        diffY = x2.size()[3] - x1.size()[3]
        # print('sizes',x1.size(),x2.size(),diffX // 2, diffX - diffX//2, diffY // 2, diffY - diffY//2)
        x1 = F.pad(x1, [diffX // 2, diffX - diffX // 2,
                        diffY // 2, diffY - diffY // 2])
        # if you have padding issues, see
        # https://github.com/HaiyongJiang/U-Net-Pytorch-Unstructured-Buggy/commit/0e854509c2cea854e247a9c615f175f76fbb2e3a
        # https://github.com/xiaopeng-liao/Pytorch-UNet/commit/8ebac70e633bac59fc22bb5195e513d5832fb3bd
        x = torch.cat([x2, x1], dim=1)
        return self.conv(x)

class OutConv(nn.Module):
    def __init__(self, in_channels, out_channels):
        super(OutConv, self).__init__()
        self.conv = nn.Conv2d(in_channels, out_channels, kernel_size=1)

    def forward(self, x):
        return self.conv(x)

class UNet(nn.Module):
    def __init__(self, n_channels, n_classes, init_filter=64, num_down=4, nearest=True, is_single=True):
        super(UNet, self).__init__()
        self.n_channels = n_channels
        self.n_classes = n_classes
        self.nearest = nearest
        self.num_down = num_down
        factor = 2 if nearest else 1
        filters = []
        for i in range(num_down + 1):
            filters.append(init_filter * (2 ** i))
        # print(filters)
        self.inc = Down_Conv(n_channels, filters[0])
        down_list = []
        up_list = []
        if is_single:
            for i in range(num_down):
                if i == num_down - 1:
                    down_list.append(Down(filters[i], filters[i + 1] // factor, mode='single_conv'))
                else:
                    down_list.append(Down(filters[i], filters[i + 1], mode='single_conv'))
            self.down = nn.ModuleList(down_list)
            for i in range(num_down):
                if i == num_down - 1:
                    up_list.append(
                        Up(filters[num_down - i], filters[num_down - i - 1] * factor, nearest, mode='double_conv'))
                else:
                    up_list.append(Up(filters[num_down - i], filters[num_down - i - 1], nearest, mode='single_conv'))
                self.up = nn.ModuleList(up_list)
        else:
            for i in range(num_down):
                if i == num_down - 1:
                    down_list.append(Down(filters[i], filters[i + 1] // factor, mode='double_conv'))
                else:
                    down_list.append(Down(filters[i], filters[i + 1], mode='double_conv'))
            self.down = nn.ModuleList(down_list)
            for i in range(num_down):
                if i == num_down - 1:
                    up_list.append(
                        Up(filters[num_down - i], filters[num_down - i - 1] * factor, nearest, mode='double_conv'))
                else:
                    up_list.append(Up(filters[num_down - i], filters[num_down - i - 1], nearest, mode='double_conv'))
                self.up = nn.ModuleList(up_list)

        self.outc = OutConv(filters[0], n_classes)

    def forward(self, x):
        down = []
        x = self.inc(x)
        down.append(x)
        for i, down_module in enumerate(self.down):
            x = down_module(x)
            down.append(x)
        for i, up_module in enumerate(self.up):
            if i == 0:
                x = up_module(down[-1], down[-2])
            else:
                x = up_module(x, down[self.num_down - 1 - i])
        logits = self.outc(x)
        return logits

def init_weights(net, init_type='normal', mean = 0.0, init_gain=0.02):
    def init_func(m):  # define the initialization function
        classname = m.__class__.__name__
        # print(classname)
        if hasattr(m, 'weight') and (classname.find('Conv') != -1 or classname.find('Linear') != -1):
            if init_type == 'normal':
                init.normal_(m.weight.data, mean, init_gain)
            elif init_type == 'xavier':
                init.xavier_normal_(m.weight.data, gain=init_gain)
            elif init_type == 'kaiming':
                init.kaiming_normal_(m.weight.data, a=0, mode='fan_in')
            elif init_type == 'orthogonal':
                init.orthogonal_(m.weight.data, gain=init_gain)
            else:
                raise NotImplementedError('initialization method [%s] is not implemented' % init_type)
            if hasattr(m, 'bias') and m.bias is not None:
                init.constant_(m.bias.data, mean)

        elif classname.find('BatchNorm') != -1:  # BatchNorm Layer's weight is not a matrix; only normal distribution applies.
                init.normal_(m.weight.data, 1.0, init_gain)
                init.constant_(m.bias.data, mean)

    print('initialize network with %s' % init_type)
    net.apply(init_func)  # apply the initialization function <init_func>

class My_model(nn.Module):
    def __init__(self, in_channels, out_channels):
        super(My_model, self).__init__()
        self.UNet = UNet(in_channels, out_channels, init_filter=8, num_down=3, nearest= True, is_single= True)
        init_weights(self.UNet, init_type= 'kaiming', mean= 0, init_gain= 0.02)
    def forward(self, img):
        out = self.UNet(img)
        return out

"""


text:
"""
    abstract type Quadrature{D,T} <: GridapType end

-[`get_coordinates(q::Quadrature)`](@ref)
-[`get_weights(q::Quadrature)`](@ref)
-[`test_quadrature`](@ref)
"""
abstract type Quadrature{D,T} <: GridapType end

# Abstract API

"""
    get_coordinates(q::Quadrature)
"""
function get_coordinates(q::Quadrature)
  @abstractmethod
end

"""
    get_weights(q::Quadrature)
"""
function get_weights(q::Quadrature)
  @abstractmethod
end

"""
    get_name(q::Quadrature)
"""
function get_name(q::Quadrature)
  @abstractmethod
end

# Default API

"""
    num_points(q::Quadrature)
"""
num_points(q::Quadrature) = length(get_weights(q))

"""
    num_point_dims(::Quadrature{D}) where D
    num_point_dims(::Type{<:Quadrature{D}}) where D
"""
num_point_dims(::Quadrature{D}) where D = D

num_point_dims(::Type{<:Quadrature{D}}) where D = D

"""
    num_dims(::Quadrature{D}) where D where D
    num_dims(::Type{<:Quadrature{D}}) where D
"""
num_dims(::Quadrature{D}) where D = D

num_dims(::Type{<:Quadrature{D}}) where D = D

# Tester

"""
    test_quadrature(q::Quadrature{D,T}) where {D,T}
"""
function test_quadrature(q::Quadrature{D,T}) where {D,T}
  x = get_coordinates(q)
  w = get_weights(q)
  @test isa(x,Vector{Point{D,T}})
  @test isa(w,Vector{T})
  @test length(x) == num_points(q)
  @test length(w) == num_points(q)
  @test D == num_dims(q)
  @test D == num_point_dims(q)
  @test D == num_dims(typeof(q))
  @test D == num_point_dims(typeof(q))
end

# Generic concrete implementation

"""
    struct GenericQuadrature{D,T} <: Quadrature{D,T}
      coordinates::Vector{Point{D,T}}
      weights::Vector{T}
      name::String
    end
"""
struct GenericQuadrature{D,T} <: Quadrature{D,T}
  coordinates::Vector{Point{D,T}}
  weights::Vector{T}
  name::String
end

function GenericQuadrature(
  coordinates::Vector{Point{D,T}}, weights::Vector{T}) where {D,T}
  name = "Unknown"
  GenericQuadrature(coordinates,weights,name)
end

get_coordinates(q::GenericQuadrature) = q.coordinates

get_weights(q::GenericQuadrature) = q.weights

get_name(q::GenericQuadrature) = q.name

function GenericQuadrature(a::Quadrature)
  GenericQuadrature(get_coordinates(a),get_weights(a),get_name(a))
end

function GenericQuadrature(a::GenericQuadrature)
  a
end

# Quadrature factory

abstract type QuadratureName end

@noinline function Quadrature(p::Polytope,name::QuadratureName,args...;kwargs...)
  @unreachable """\n
  Undefined factory function Quadrature for name $name and the given arguments.
  """
end

Quadrature(name::QuadratureName,args...;kwargs...) = (name, args, kwargs)

"""
    Quadrature(polytope::Polytope{D},degree) where D
"""
function Quadrature(p::Polytope,degree)
  if is_n_cube(p)
    quad = Quadrature(p,tensor_product,degree)
  elseif is_simplex(p)
    D = num_dims(p)
    #if (D==2 && degree in keys(_strang_tri_k2n)) ||
    #   (D==3 && degree in keys(_strang_tet_k2n))
    # For the moment strang only in 3d since
    # there are some 2d strang quadratures that are not accurate
    # as implemented here (to investigate why)
    if (D==3 && degree in keys(_strang_tet_k2n))
      quad = Quadrature(p,strang,degree)
    else
      quad = Quadrature(p,duffy,degree)
    end
  else
    @notimplemented "Quadratures only implemented for n-cubes and simplices"
  end
  quad
end