text:
"""
$(SIGNATURES)

Update solver for boundary cells
"""
function update_boundary!(
    KS::AbstractSolverSet,
    ctr::AV{TC},
    face::AV{TF},
    dt,
    residual;
    bc,
) where {TC<:Union{ControlVolume,ControlVolume1D},TF<:Union{Interface,Interface1D}}

    resL = zero(ctr[1].w)
    avgL = zero(ctr[1].w)
    resR = zero(ctr[1].w)
    avgR = zero(ctr[1].w)

    if bc[1] != :fix
        i = 1
        j = KS.pSpace.nx
        if KS.set.nSpecies == 1
            step!(
                face[i].fw,
                ctr[i].w,
                ctr[i].prim,
                face[i+1].fw,
                KS.gas.γ,
                KS.ps.dx[i],
                resL,
                avgL,
            )
        elseif KS.set.nSpecies == 2
            step!(
                face[i].fw,
                ctr[i].w,
                ctr[i].prim,
                face[i+1].fw,
                KS.gas.γ,
                KS.gas.mi,
                KS.gas.ni,
                KS.gas.me,
                KS.gas.ne,
                KS.gas.Kn[1],
                KS.ps.dx[i],
                dt,
                resL,
                avgL,
            )
        end
    end

    if bc[2] != :fix
        j = KS.pSpace.nx
        if KS.set.nSpecies == 1
            step!(
                face[j].fw,
                ctr[j].w,
                ctr[j].prim,
                face[j+1].fw,
                KS.gas.γ,
                KS.ps.dx[j],
                resR,
                avgR,
            )
        elseif KS.set.nSpecies == 2
            step!(
                face[j].fw,
                ctr[j].w,
                ctr[j].prim,
                face[j+1].fw,
                KS.gas.γ,
                KS.gas.mi,
                KS.gas.ni,
                KS.gas.me,
                KS.gas.ne,
                KS.gas.Kn[1],
                KS.ps.dx[j],
                dt,
                resR,
                avgR,
            )
        end
    end

    @. residual += sqrt((resL + resR) * 2) / (avgL + avgR + 1.e-7)

    ng = 1 - first(eachindex(KS.pSpace.x))
    if bc[1] == :period
        bc_period!(ctr, ng)
    elseif bc[1] == :extra
        bc_extra!(ctr, ng; dirc = :xl)
    elseif bc[1] == :balance
        bc_balance!(ctr[0], ctr[1], ctr[2])
    end
    if bc[2] == :extra
        bc_extra!(ctr, ng; dirc = :xr)
    elseif bc[2] == :balance
        bc_balance!(ctr[KS.ps.nx+1], ctr[KS.ps.nx], ctr[KS.ps.nx-1])
    end

    return nothing

end

function update_boundary!(
    KS::AbstractSolverSet,
    ctr::AV{TC},
    face::AV{TF},
    dt,
    residual;
    coll = symbolize(KS.set.collision),
    bc,
) where {TC<:Union{ControlVolume1F,ControlVolume1D1F},TF<:Union{Interface1F,Interface1D1F}}

    resL = zero(ctr[1].w)
    avgL = zero(ctr[1].w)
    resR = zero(ctr[1].w)
    avgR = zero(ctr[1].w)

    if bc[1] != :fix
        i = 1

        if KS.set.space[5:6] == "1v"
            step!(
                face[i].fw,
                face[i].ff,
                ctr[i].w,
                ctr[i].prim,
                ctr[i].f,
                face[i+1].fw,
                face[i+1].ff,
                KS.vSpace.u,
                KS.vSpace.weights,
                KS.gas.γ,
                KS.gas.μᵣ,
                KS.gas.ω,
                KS.gas.Pr,
                KS.ps.dx[i],
                dt,
                resL,
                avgL,
                coll,
            )
        elseif KS.set.space[5:6] == "3v"
            step!(
                face[i].fw,
                face[i].ff,
                ctr[i].w,
                ctr[i].prim,
                ctr[i].f,
                face[i+1].fw,
                face[i+1].ff,
                KS.vSpace.u,
                KS.vSpace.v,
                KS.vSpace.w,
                KS.vSpace.weights,
                KS.gas.γ,
                KS.gas.μᵣ,
                KS.gas.ω,
                KS.gas.Pr,
                KS.ps.dx[i],
                dt,
                resL,
                avgL,
                coll,
            )
        end
    end

    if bc[2] != :fix
        j = KS.pSpace.nx

        if KS.set.space[5:6] == "1v"
            step!(
                face[j].fw,
                face[j].ff,
                ctr[j].w,
                ctr[j].prim,
                ctr[j].f,
                face[j+1].fw,
                face[j+1].ff,
                KS.vSpace.u,
                KS.vSpace.weights,
                KS.gas.γ,
                KS.gas.μᵣ,
                KS.gas.ω,
                KS.gas.Pr,
                KS.ps.dx[j],
                dt,
                resR,
                avgR,
                coll,
            )
        elseif KS.set.space[5:6] == "3v"
            step!(
                face[j].fw,
                face[j].ff,
                ctr[j].w,
                ctr[j].prim,
                ctr[j].f,
                face[j+1].fw,
                face[j+1].ff,
                KS.vSpace.u,
                KS.vSpace.v,
                KS.vSpace.w,
                KS.vSpace.weights,
                KS.gas.γ,
                KS.gas.μᵣ,
                KS.gas.ω,
                KS.gas.Pr,
                KS.ps.dx[j],
                dt,
                resL,
                avgL,
                coll,
            )
        end
    end

    for i in eachindex(residual)
        residual[i] += sqrt((resL[i] + resR[i]) * 2) / (avgL[i] + avgR[i] + 1.e-7)
    end

    ng = 1 - first(eachindex(KS.pSpace.x))
    if bc[1] == :period
        bc_period!(ctr, ng)
    elseif bc[1] == :extra
        bc_extra!(ctr, ng; dirc = :xl)
    elseif bc[1] == :balance
        bc_balance!(ctr[0], ctr[1], ctr[2])
    end
    if bc[2] == :extra
        bc_extra!(ctr, ng; dirc = :xr)
    elseif bc[2] == :balance
        bc_balance!(ctr[KS.ps.nx+1], ctr[KS.ps.nx], ctr[KS.ps.nx-1])
    end

    return nothing

end

function update_boundary!(
    KS::AbstractSolverSet,
    ctr::AV{TC},
    face::AV{TF},
    dt,
    residual;
    coll = symbolize(KS.set.collision),
    bc,
) where {TC<:Union{ControlVolume2F,ControlVolume1D2F},TF<:Union{Interface2F,Interface1D2F}}

    resL = zero(ctr[1].w)
    avgL = zero(ctr[1].w)
    resR = zero(ctr[1].w)
    avgR = zero(ctr[1].w)

    if bc[1] != :fix
        i = 1
        if KS.set.nSpecies == 1
            step!(
                face[i].fw,
                face[i].fh,
                face[i].fb,
                ctr[i].w,
                ctr[i].prim,
                ctr[i].h,
                ctr[i].b,
                face[i+1].fw,
                face[i+1].fh,
                face[i+1].fb,
                KS.vSpace.u,
                KS.vSpace.weights,
                KS.gas.K,
                KS.gas.γ,
                KS.gas.μᵣ,
                KS.gas.ω,
                KS.gas.Pr,
                KS.ps.dx[i],
                dt,
                resL,
                avgL,
                coll,
            )
        else
            step!(
                face[i].fw,
                face[i].fh,
                face[i].fb,
                ctr[i].w,
                ctr[i].prim,
                ctr[i].h,
                ctr[i].b,
                face[i+1].fw,
                face[i+1].fh,
                face[i+1].fb,
                KS.vSpace.u,
                KS.vSpace.weights,
                KS.gas.K,
                KS.gas.γ,
                KS.gas.mi,
                KS.gas.ni,
                KS.gas.me,
                KS.gas.ne,
                KS.gas.Kn[1],
                KS.gas.Pr,
                KS.ps.dx[i],
                dt,
                resL,
                avgL,
                coll,
            )
        end
    end

    if bc[2] != :fix
        j = KS.pSpace.nx
        if KS.set.nSpecies == 1
            step!(
                face[j].fw,
                face[j].fh,
                face[j].fb,
                ctr[j].w,
                ctr[j].prim,
                ctr[j].h,
                ctr[j].b,
                face[j+1].fw,
                face[j+1].fh,
                face[j+1].fb,
                KS.vSpace.u,
                KS.vSpace.weights,
                KS.gas.K,
                KS.gas.γ,
                KS.gas.μᵣ,
                KS.gas.ω,
                KS.gas.Pr,
                KS.ps.dx[j],
                dt,
                resL,
                avgL,
                coll,
            )
        else
            step!(
                face[j].fw,
                face[j].fh,
                face[j].fb,
                ctr[j].w,
                ctr[j].prim,
                ctr[j].h,
                ctr[j].b,
                face[j+1].fw,
                face[j+1].fh,
                face[j+1].fb,
                KS.vSpace.u,
                KS.vSpace.weights,
                KS.gas.K,
                KS.gas.γ,
                KS.gas.mi,
                KS.gas.ni,
                KS.gas.me,
                KS.gas.ne,
                KS.gas.Kn[1],
                KS.gas.Pr,
                KS.ps.dx[j],
                dt,
                resL,
                avgL,
                coll,
            )
        end
    end

    for i in eachindex(residual)
        residual[i] += sqrt((resL[i] + resR[i]) * 2) / (avgL[i] + avgR[i] + 1.e-7)
    end

    ng = 1 - first(eachindex(KS.pSpace.x))
    if bc[1] == :period
        bc_period!(ctr, ng)
    elseif bc[1] == :extra
        bc_extra!(ctr, ng; dirc = :xl)
    elseif bc[1] == :balance
        bc_balance!(ctr[0], ctr[1], ctr[2])
    end
    if bc[2] == :extra
        bc_extra!(ctr, ng; dirc = :xr)
    elseif bc[2] == :balance
        bc_balance!(ctr[KS.ps.nx+1], ctr[KS.ps.nx], ctr[KS.ps.nx-1])
    end

    return nothing

end

function update_boundary!(
    KS::AbstractSolverSet,
    ctr::AV{TC},
    face::AV{TF},
    dt,
    residual;
    coll = symbolize(KS.set.collision)::Symbol,
    bc,
    isMHD = false::Bool,
) where {TC<:Union{ControlVolume3F,ControlVolume1D3F},TF<:Union{Interface3F,Interface1D3F}}

    resL = zero(ctr[1].w)
    avgL = zero(ctr[1].w)
    resR = zero(ctr[1].w)
    avgR = zero(ctr[1].w)

    if bc[1] != :fix
        i = 1
        step!(KS, face[i], ctr[i], face[i+1], KS.ps.dx[i], dt, resL, avgL, coll, isMHD)
    end

    if bc[2] != :fix
        j = KS.pSpace.nx
        step!(KS, face[j], ctr[j], face[j+1], KS.ps.dx[j], dt, resR, avgR, coll, isMHD)
    end

    for i in eachindex(residual)
        residual[i] += sqrt((resL[i] + resR[i]) * 2) / (avgL[i] + avgR[i] + 1.e-7)
    end

    ng = 1 - first(eachindex(KS.pSpace.x))
    if bc[1] == :period
        bc_period!(ctr, ng)
    elseif bc[1] == :extra
        bc_extra!(ctr, ng; dirc = :xl)
    elseif bc[1] == :balance
        bc_balance!(ctr[0], ctr[1], ctr[2])
    end
    if bc[2] == :extra
        bc_extra!(ctr, ng; dirc = :xr)
    elseif bc[2] == :balance
        bc_balance!(ctr[KS.ps.nx+1], ctr[KS.ps.nx], ctr[KS.ps.nx-1])
    end

    return nothing

end

function update_boundary!(
    KS::AbstractSolverSet,
    ctr::AV{TC},
    face::AV{TF},
    dt,
    residual;
    coll = symbolize(KS.set.collision)::Symbol,
    bc,
    isMHD = false::Bool,
) where {TC<:Union{ControlVolume4F,ControlVolume1D4F},TF<:Union{Interface4F,Interface1D4F}}

    resL = zero(ctr[1].w)
    avgL = zero(ctr[1].w)
    resR = zero(ctr[1].w)
    avgR = zero(ctr[1].w)

    if bc[1] != :fix
        i = 1
        step!(KS, face[i], ctr[i], face[i+1], KS.ps.dx[i], dt, resL, avgL, coll, isMHD)
    end

    if bc[2] != :fix
        j = KS.pSpace.nx
        step!(KS, face[j], ctr[j], face[j+1], KS.ps.dx[j], dt, resR, avgR, coll, isMHD)
    end

    for i in eachindex(residual)
        residual[i] += sqrt((resL[i] + resR[i]) * 2) / (avgL[i] + avgR[i] + 1.e-7)
    end

    ng = 1 - first(eachindex(KS.pSpace.x))
    if bc[1] == :period
        bc_period!(ctr, ng)
    elseif bc[1] == :extra
        bc_extra!(ctr, ng; dirc = :xl)
    elseif bc[1] == :balance
        bc_balance!(ctr[0], ctr[1], ctr[2])
    end
    if bc[2] == :extra
        bc_extra!(ctr, ng; dirc = :xr)
    elseif bc[2] == :balance
        bc_balance!(ctr[KS.ps.nx+1], ctr[KS.ps.nx], ctr[KS.ps.nx-1])
    end

    return nothing

end

function update_boundary!(
    KS::AbstractSolverSet,
    ctr::AM{TC},
    a1face::AM{TF},
    a2face::AM{TF},
    dt,
    residual;
    coll = symbolize(KS.set.collision)::Symbol,
    bc,
) where {TC<:Union{ControlVolume,ControlVolume2D},TF<:Union{Interface,Interface2D}}

    nx, ny, dx, dy = begin
        if KS.ps isa CSpace2D
            KS.ps.nr, KS.ps.nθ, KS.ps.dr, KS.ps.darc
        else
            KS.ps.nx, KS.ps.ny, KS.ps.dx, KS.ps.dy
        end
    end

    resL = zero(ctr[1].w)
    avgL = zero(ctr[1].w)
    resR = zero(ctr[1].w)
    avgR = zero(ctr[1].w)
    resU = zero(ctr[1].w)
    avgU = zero(ctr[1].w)
    resD = zero(ctr[1].w)
    avgD = zero(ctr[1].w)

    if bc[1] != :fix
        @inbounds for j = 1:ny
            step!(
                ctr[1, j].w,
                ctr[1, j].prim,
                a1face[1, j].fw,
                a1face[2, j].fw,
                a2face[1, j].fw,
                a2face[1, j+1].fw,
                KS.gas.γ,
                dx[1, j] * dy[1, j],
                resL,
                avgL,
                coll,
            )
        end
    end

    if bc[2] != :fix
        @inbounds for j = 1:ny
            step!(
                ctr[nx, j].w,
                ctr[nx, j].prim,
                a1face[nx, j].fw,
                a1face[nx+1, j].fw,
                a2face[nx, j].fw,
                a2face[nx, j+1].fw,
                KS.gas.γ,
                dx[nx, j] * dy[nx, j],
                resR,
                avgR,
                coll,
            )
        end
    end

    if bc[3] != :fix
        @inbounds for i = 2:nx-1 # skip overlap
            step!(
                ctr[i, 1].w,
                ctr[i, 1].prim,
                a1face[i, 1].fw,
                a1face[i+1, 1].fw,
                a2face[i, 1].fw,
                a2face[i, 2].fw,
                KS.gas.γ,
                dx[i, 1] * dy[i, 1],
                resD,
                avgD,
                coll,
            )
        end
    end

    if bc[4] != :fix
        @inbounds for i = 2:nx-1 # skip overlap
            step!(
                ctr[i, ny].w,
                ctr[i, ny].prim,
                a1face[i, ny].fw,
                a1face[i+1, ny].fw,
                a2face[i, ny].fw,
                a2face[i, ny+1].fw,
                KS.gas.γ,
                dx[i, ny] * dy[i, ny],
                resU,
                avgU,
                coll,
            )
        end
    end

    for i in eachindex(residual)
        residual[i] +=
            sqrt((resL[i] + resR[i] + resU[i] + resD[i]) * 2) /
            (avgL[i] + avgR[i] + avgU[i] + avgD[i] + 1.e-7)
    end

    ngx = 1 - first(eachindex(KS.pSpace.x[:, 1]))
    if bc[1] == :period
        bc_period!(ctr, ngx; dirc = :x)
    elseif bc[1] in (:extra, :mirror)
        bcfun = eval(Symbol("bc_" * string(bc[1]) * "!"))
        bcfun(ctr, ngx; dirc = :xl)
    end
    if bc[2] in (:extra, :mirror)
        bcfun = eval(Symbol("bc_" * string(bc[2]) * "!"))
        bcfun(ctr, ngx; dirc = :xr)
    end

    ngy = 1 - first(eachindex(KS.pSpace.y[1, :]))
    if bc[3] == :period
        bc_period!(ctr, ngy; dirc = :y)
    elseif bc[3] in (:extra, :mirror)
        bcfun = eval(Symbol("bc_" * string(bc[3]) * "!"))
        bcfun(ctr, ngy; dirc = :yl)
    end
    if bc[4] in (:extra, :mirror)
        bcfun = eval(Symbol("bc_" * string(bc[4]) * "!"))
        bcfun(ctr, ngy; dirc = :yr)
    end

    return nothing

end

function update_boundary!(
    KS::AbstractSolverSet,
    ctr::AM{TC},
    a1face::AM{TF},
    a2face::AM{TF},
    dt,
    residual;
    coll = symbolize(KS.set.collision)::Symbol,
    bc,
) where {TC<:Union{ControlVolume1F,ControlVolume2D1F},TF<:Union{Interface1F,Interface2D1F}}

    nx, ny, dx, dy = begin
        if KS.ps isa CSpace2D
            KS.ps.nr, KS.ps.nθ, KS.ps.dr, KS.ps.darc
        else
            KS.ps.nx, KS.ps.ny, KS.ps.dx, KS.ps.dy
        end
    end

    resL = zero(ctr[1].w)
    avgL = zero(ctr[1].w)
    resR = zero(ctr[1].w)
    avgR = zero(ctr[1].w)
    resU = zero(ctr[1].w)
    avgU = zero(ctr[1].w)
    resD = zero(ctr[1].w)
    avgD = zero(ctr[1].w)

    if bc[1] != :fix
        @inbounds for j = 1:ny
            step!(
                ctr[1, j].w,
                ctr[1, j].prim,
                ctr[1, j].f,
                a1face[1, j].fw,
                a1face[1, j].ff,
                a1face[2, j].fw,
                a1face[2, j].ff,
                a2face[1, j].fw,
                a2face[1, j].ff,
                a2face[1, j+1].fw,
                a2face[1, j+1].ff,
                KS.vSpace.u,
                KS.vSpace.v,
                KS.vSpace.weights,
                KS.gas.γ,
                KS.gas.μᵣ,
                KS.gas.ω,
                KS.gas.Pr,
                dx[1, j] * dy[1, j],
                dt,
                resL,
                avgL,
                coll,
            )
        end
    end
    if bc[2] != :fix
        @inbounds for j = 1:ny
            step!(
                ctr[nx, j].w,
                ctr[nx, j].prim,
                ctr[nx, j].f,
                a1face[nx, j].fw,
                a1face[nx, j].ff,
                a1face[nx+1, j].fw,
                a1face[nx+1, j].ff,
                a2face[nx, j].fw,
                a2face[nx, j].ff,
                a2face[nx, j+1].fw,
                a2face[nx, j+1].ff,
                KS.vSpace.u,
                KS.vSpace.v,
                KS.vSpace.weights,
                KS.gas.γ,
                KS.gas.μᵣ,
                KS.gas.ω,
                KS.gas.Pr,
                dx[nx, j] * dy[nx, j],
                dt,
                resR,
                avgR,
                coll,
            )
        end
    end
    if bc[3] != :fix
        @inbounds for i = 2:nx-1 # skip overlap
            step!(
                ctr[i, 1].w,
                ctr[i, 1].prim,
                ctr[i, 1].f,
                a1face[i, 1].fw,
                a1face[i, 1].ff,
                a1face[i+1, 1].fw,
                a1face[i+1, 1].ff,
                a2face[i, 1].fw,
                a2face[i, 1].ff,
                a2face[i, 2].fw,
                a2face[i, 2].ff,
                KS.vSpace.u,
                KS.vSpace.v,
                KS.vSpace.weights,
                KS.gas.γ,
                KS.gas.μᵣ,
                KS.gas.ω,
                KS.gas.Pr,
                dx[i, 1] * dy[i, 1],
                dt,
                resD,
                avgD,
                coll,
            )
        end
    end
    if bc[4] != :fix
        @inbounds for i = 2:nx-1 # skip overlap
            step!(
                ctr[i, ny].w,
                ctr[i, ny].prim,
                ctr[i, ny].f,
                a1face[i, ny].fw,
                a1face[i, ny].ff,
                a1face[i+1, ny].fw,
                a1face[i+1, ny].ff,
                a2face[i, ny].fw,
                a2face[i, ny].ff,
                a2face[i, ny+1].fw,
                a2face[i, ny+1].ff,
                KS.vSpace.u,
                KS.vSpace.v,
                KS.vSpace.weights,
                KS.gas.γ,
                KS.gas.μᵣ,
                KS.gas.ω,
                KS.gas.Pr,
                dx[i, ny] * dy[i, ny],
                dt,
                resU,
                avgU,
                coll,
            )
        end
    end

    for i in eachindex(residual)
        residual[i] +=
            sqrt((resL[i] + resR[i] + resU[i] + resD[i]) * 2) /
            (avgL[i] + avgR[i] + avgU[i] + avgD[i] + 1.e-7)
    end

    ngx = 1 - first(eachindex(KS.pSpace.x[:, 1]))
    if bc[1] == :period
        bc_period!(ctr, ngx; dirc = :x)
    elseif bc[1] in (:extra, :mirror)
        bcfun = eval(Symbol("bc_" * string(bc[1]) * "!"))
        bcfun(ctr, ngx; dirc = :xl)
    end
    if bc[2] in (:extra, :mirror)
        bcfun = eval(Symbol("bc_" * string(bc[2]) * "!"))
        bcfun(ctr, ngx; dirc = :xr)
    end

    ngy = 1 - first(eachindex(KS.pSpace.y[1, :]))
    if bc[3] == :period
        bc_period!(ctr, ngy; dirc = :y)
    elseif bc[3] in (:extra, :mirror)
        bcfun = eval(Symbol("bc_" * string(bc[3]) * "!"))
        bcfun(ctr, ngy; dirc = :yl)
    end
    if bc[4] in (:extra, :mirror)
        bcfun = eval(Symbol("bc_" * string(bc[4]) * "!"))
        bcfun(ctr, ngy; dirc = :yr)
    end

    return nothing

end

function update_boundary!(
    KS::AbstractSolverSet,
    ctr::AM{TC},
    a1face::AM{TF},
    a2face::AM{TF},
    dt,
    residual;
    coll = symbolize(KS.set.collision)::Symbol,
    bc,
) where {TC<:Union{ControlVolume2F,ControlVolume2D2F},TF<:Union{Interface2F,Interface2D2F}}

    nx, ny, dx, dy = begin
        if KS.ps isa CSpace2D
            KS.ps.nr, KS.ps.nθ, KS.ps.dr, KS.ps.darc
        else
            KS.ps.nx, KS.ps.ny, KS.ps.dx, KS.ps.dy
        end
    end

    resL = zero(ctr[1].w)
    avgL = zero(ctr[1].w)
    resR = zero(ctr[1].w)
    avgR = zero(ctr[1].w)
    resU = zero(ctr[1].w)
    avgU = zero(ctr[1].w)
    resD = zero(ctr[1].w)
    avgD = zero(ctr[1].w)

    if bc[1] != :fix
        @inbounds for j = 1:ny
            step!(
                ctr[1, j].w,
                ctr[1, j].prim,
                ctr[1, j].h,
                ctr[1, j].b,
                a1face[1, j].fw,
                a1face[1, j].fh,
                a1face[1, j].fb,
                a1face[2, j].fw,
                a1face[2, j].fh,
                a1face[2, j].fb,
                a2face[1, j].fw,
                a2face[1, j].fh,
                a2face[1, j].fb,
                a2face[1, j+1].fw,
                a2face[1, j+1].fh,
                a2face[1, j+1].fb,
                KS.vSpace.u,
                KS.vSpace.v,
                KS.vSpace.weights,
                KS.gas.K,
                KS.gas.γ,
                KS.gas.μᵣ,
                KS.gas.ω,
                KS.gas.Pr,
                dx[1, j] * dy[1, j],
                dt,
                resL,
                avgL,
                coll,
            )
        end
    end
    if bc[2] != :fix
        @inbounds for j = 1:ny
            step!(
                ctr[nx, j].w,
                ctr[nx, j].prim,
                ctr[nx, j].h,
                ctr[nx, j].b,
                a1face[nx, j].fw,
                a1face[nx, j].fh,
                a1face[nx, j].fb,
                a1face[nx+1, j].fw,
                a1face[nx+1, j].fh,
                a1face[nx+1, j].fb,
                a2face[nx, j].fw,
                a2face[nx, j].fh,
                a2face[nx, j].fb,
                a2face[nx, j+1].fw,
                a2face[nx, j+1].fh,
                a2face[nx, j+1].fb,
                KS.vSpace.u,
                KS.vSpace.v,
                KS.vSpace.weights,
                KS.gas.K,
                KS.gas.γ,
                KS.gas.μᵣ,
                KS.gas.ω,
                KS.gas.Pr,
                dx[nx, j] * dy[nx, j],
                dt,
                resR,
                avgR,
                coll,
            )
        end
    end
    if bc[3] != :fix
        @inbounds for i = 2:nx-1 # skip overlap
            step!(
                ctr[i, 1].w,
                ctr[i, 1].prim,
                ctr[i, 1].h,
                ctr[i, 1].b,
                a1face[i, 1].fw,
                a1face[i, 1].fh,
                a1face[i, 1].fb,
                a1face[i+1, 1].fw,
                a1face[i+1, 1].fh,
                a1face[i+1, 1].fb,
                a2face[i, 1].fw,
                a2face[i, 1].fh,
                a2face[i, 1].fb,
                a2face[i, 2].fw,
                a2face[i, 2].fh,
                a2face[i, 2].fb,
                KS.vSpace.u,
                KS.vSpace.v,
                KS.vSpace.weights,
                KS.gas.K,
                KS.gas.γ,
                KS.gas.μᵣ,
                KS.gas.ω,
                KS.gas.Pr,
                dx[i, 1] * dy[i, 1],
                dt,
                resD,
                avgD,
                coll,
            )
        end
    end
    if bc[4] != :fix
        @inbounds for i = 2:nx-1 # skip overlap
            step!(
                ctr[i, ny].w,
                ctr[i, ny].prim,
                ctr[i, ny].h,
                ctr[i, ny].b,
                a1face[i, ny].fw,
                a1face[i, ny].fh,
                a1face[i, ny].fb,
                a1face[i+1, ny].fw,
                a1face[i+1, ny].fh,
                a1face[i+1, ny].fb,
                a2face[i, ny].fw,
                a2face[i, ny].fh,
                a2face[i, ny].fb,
                a2face[i, ny+1].fw,
                a2face[i, ny+1].fh,
                a2face[i, ny+1].fb,
                KS.vSpace.u,
                KS.vSpace.v,
                KS.vSpace.weights,
                KS.gas.K,
                KS.gas.γ,
                KS.gas.μᵣ,
                KS.gas.ω,
                KS.gas.Pr,
                dx[i, ny] * dy[i, ny],
                dt,
                resU,
                avgU,
                coll,
            )
        end
    end

    for i in eachindex(residual)
        residual[i] +=
            sqrt((resL[i] + resR[i] + resU[i] + resD[i]) * 2) /
            (avgL[i] + avgR[i] + avgU[i] + avgD[i] + 1.e-7)
    end

    ngx = 1 - first(eachindex(KS.pSpace.x[:, 1]))
    if bc[1] == :period
        bc_period!(ctr, ngx; dirc = :x)
    elseif bc[1] in (:extra, :mirror)
        bcfun = eval(Symbol("bc_" * string(bc[1]) * "!"))
        bcfun(ctr, ngx; dirc = :xl)
    end
    if bc[2] in (:extra, :mirror)
        bcfun = eval(Symbol("bc_" * string(bc[2]) * "!"))
        bcfun(ctr, ngx; dirc = :xr)
    end

    ngy = 1 - first(eachindex(KS.pSpace.y[1, :]))
    if bc[3] == :period
        bc_period!(ctr, ngy; dirc = :y)
    elseif bc[3] in (:extra, :mirror)
        bcfun = eval(Symbol("bc_" * string(bc[3]) * "!"))
        bcfun(ctr, ngy; dirc = :yl)
    end
    if bc[4] in (:extra, :mirror)
        bcfun = eval(Symbol("bc_" * string(bc[4]) * "!"))
        bcfun(ctr, ngy; dirc = :yr)
    end

    return nothing

end

function update_boundary!(
    KS::AbstractSolverSet,
    ctr::AV{TC},
    face::AV{TF},
    dt,
    residual;
    coll::Symbol,
    bc,
) where {TC<:ControlVolumeUS,TF<:Interface2D}

    for i in eachindex(KS.ps.cellType)
        if KS.ps.cellType[i] == 3
            ids = KS.ps.cellNeighbors[i, :]
            deleteat!(ids, findall(x -> x == -1, ids))
            id1, id2 = ids
            ctr[i].w .= 0.5 .* (ctr[id1].w .+ ctr[id2].w)
            ctr[i].prim .= KitBase.conserve_prim(ctr[i].w, KS.gas.γ)
        end
    end

    return nothing

end

function update_boundary!(
    KS::AbstractSolverSet,
    ctr::AV{TC},
    face::AV{TF},
    dt,
    residual;
    coll::Symbol,
    bc,
) where {TC<:ControlVolumeUS1F,TF<:Interface2D1F}

    for i in eachindex(KS.ps.cellType)
        if KS.ps.cellType[i] == 3
            ids = KS.ps.cellNeighbors[i, :]
            deleteat!(ids, findall(x -> x == -1, ids))
            id1, id2 = ids
            ctr[i].w .= 0.5 .* (ctr[id1].w .+ ctr[id2].w)
            ctr[i].f .= 0.5 .* (ctr[id1].f .+ ctr[id2].f)
            ctr[i].prim .= KitBase.conserve_prim(ctr[i].w, KS.gas.γ)
        end
    end

    return nothing

end

function update_boundary!(
    KS::AbstractSolverSet,
    ctr::AV{TC},
    face::AV{TF},
    dt,
    residual;
    coll::Symbol,
    bc,
) where {TC<:ControlVolumeUS2F,TF<:Interface2D2F}

    for i in eachindex(KS.ps.cellType)
        if KS.ps.cellType[i] == 3
            ids = KS.ps.cellNeighbors[i, :]
            deleteat!(ids, findall(x -> x == -1, ids))
            id1, id2 = ids
            ctr[i].w .= 0.5 .* (ctr[id1].w .+ ctr[id2].w)
            ctr[i].h .= 0.5 .* (ctr[id1].h .+ ctr[id2].h)
            ctr[i].b .= 0.5 .* (ctr[id1].b .+ ctr[id2].b)
            ctr[i].prim .= KitBase.conserve_prim(ctr[i].w, KS.gas.γ)
        end
    end

    return nothing

end


text:
# We shall say that an n-digit number is pandigital if it makes use of all the
# digits 1 to n exactly once. For example, 2143 is a 4-digit pandigital and is
# also prime.
#
# What is the largest n-digit pandigital prime that exists?

using ProjectEulerSolutions

# Helper function for testing primality
function testprime(vals::Array)
    return isprime(sum(vals .* (10 .^ (0:length(vals)-1))))
end

# Creates pandigital permutations, and tests for primality.  Pandigitals of
# length 8 and 9 are divisible by 3 and thus cannot be prime.
function p041solution_permutation()::Integer
    for n in 7:-1:1
        matrix = zeros(Integer, 0)
        permutation_matrix!(collect(1:n), matrix, testprime)
        primes = transpose(10 .^ (0:n-1)) * reshape(matrix, (n, fld(length(matrix), n)))
        return maximum(primes)
    end
    return -1
end

# Brute force solution, checking every prime up to 7,654,321, since all
# pandigitals of length 8 and 9 are divisible by 3 and thus cannot be prime.
function p041solution_sieve()::Integer
    n = 7_654_321
    primes = sieve_eratosthenes(n)
    max_pandigital_prime = 2143
    lower_bound = 1e6

    for p in primes
        digs = digits(p)
        # max digit = length, min digit = 1, and every digit unique
        if length(digs) == maximum(digs) && minimum(digs) == 1 && length(Set(digs)) == length(digs)
            max_pandigital_prime = p
        end
    end

    return max_pandigital_prime
end

p041 = Problems.Problem(Dict("Permutation" => p041solution_permutation,
                             "Sieve" => p041solution_sieve))

Problems.benchmark(p041)

ext: jl

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
# ---
# title: 724. Find Pivot Index
# id: problem724
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/find-pivot-index/description/>
# hidden: true
# ---
#
# Given an array of integers `nums`, write a method that returns the "pivot"
# index of this array.
#
# We define the pivot index as the index where the sum of all the numbers to the
# left of the index is equal to the sum of all the numbers to the right of the
# index.
#
# If no such index exists, we should return -1. If there are multiple pivot
# indexes, you should return the left-most pivot index.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,7,3,6,5,6]
#     Output: 3
#     Explanation:
#     The sum of the numbers to the left of index 3 (nums[3] = 6) is equal to the sum of numbers to the right of index 3.
#     Also, 3 is the first index where this occurs.
#
#
# **Example 2:**
#
#
#
#     Input: nums = [1,2,3]
#     Output: -1
#     Explanation:
#     There is no index that satisfies the conditions in the problem statement.
#
#
#
#
# **Constraints:**
#
#   * The length of `nums` will be in the range `[0, 10000]`.
#   * Each element `nums[i]` will be an integer in the range `[-1000, 1000]`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
export ## Types
       Cooling,
       ## Methods
       label,
       rate,
       liouvillian_left,
       liouvillian_right,
       liouvillian_bilat

label(d::Dissipation) = d.label
rate(d::Dissipation) = d.rate

# NOTE: Only single body dissipators are currently supported

struct Cooling <: Dissipation
  label::AbstractString
  rate::Real
  system::QSystem
end

function liouvillian_left(d::Cooling, t::Real)
  # Note these return the DUAL to the time propagator generator
  local a, n, r
  a = lowering(d.system)
  n = a'*a
  r = rate(d)
  return -1/2*r*n
end

function liouvillian_right(d::Cooling, t::Real)
  # Note these return the DUAL to the time propagator generator
  local a, n, r
  a = lowering(d.system)
  n = a'*a
  r = rate(d)
  return -1/2*r*n
end

function liouvillian_bilat(d::Cooling, t::Real)
  # Note these return the DUAL to the time propagator generator
  local a, n, r
  a = lowering(d.system)
  r = rate(d)
  return r*kron(transpose(a),a')
end


text:
using DifferentialEquations, Optim, DiffEqFlux, DiffEqSensitivity, Plots
using Zygote, ForwardDiff, ReverseDiff
# using BlackBoxOptim#,GalacticOptim
using BenchmarkTools
using RecursiveArrayTools
# -------------------------------------------------------------------------

mutable struct F_t_Data
    F # forces vector
    t # lifetimes vector
    kd # kernel density of F
end

mutable struct M
    model   # model function
    u₀      # initial conditions
    u₀_type # string defining initial condtion type
    tspan   # time span
    dens    # site densities
    cons    # constant parameters
    f       # force
    solver  # ode solver
end

# function solve_prob(M,p;tᵢ=[])
#     f = (du,u,p,t) -> M.model(du,u,p,t,M.dens,M.cons)
#
#     # k = p[1:end-1]
#     # a = p[end]
#     u₀ = initial_u(p,M)
#
#     prob = ODEProblem(f,u₀,M.tspan,p)
#     sol = solve(prob,Rodas5(),tstops=tᵢ,dtmax=0.5,
#                             # callback=PositiveDomain(),
#                             isoutofdomain=(u,p,t)->any(x->x<0,u),
#                             abstol=1e-10,
#                             reltol=1e-10,
#                             verbose=false)
# end
function solve_prob(M,p;tᵢ=[])

    f = (du,u,p,t) -> M.model(du,u,p,t,M.dens,M.cons;f=M.f)
    u₀ = initial_u(p,M)

    prob = ODEProblem(f,u₀,M.tspan,p)
    sol = solve(prob,M.solver(),tstops=tᵢ,dtmax=0.5,
                            # callback=PositiveDomain(),
                            isoutofdomain=(u,p,t)->any(x->x<0,u),
                            abstol=1e-10,
                            reltol=1e-10,
                            verbose=false)
end

function mle_loss(p,data,M)

    L = []
    for (Fᵢ,tᵢ) in zip(data.F,data.t)

        M.f = Fᵢ
        sol = solve_prob(M,p;tᵢ)

        #check whether integration is successful
        if sol.retcode != :Success || sol.t[end] < tᵢ
            push!(L,0.0) # Likelihood of model is 0
        else

            ŷᵢ = sum(sol(tᵢ))

            mean_t = mean_lifetime_linear(sol)
            # Likelihood
            Lᵢ = ŷᵢ / mean_t
            if Lᵢ < 0
                @show Lᵢ, mean_t, ŷᵢ, tᵢ
            end

            push!(L,Lᵢ)
        end
    end
    # logL = logLikelihood(L)
    # logL = wLogLikelihood(L,1 ./ pdf(data.kd,data.F))
    # logL = wLogLikelihood(L,data.t)
    # logL = wLogLikelihood(L,log.(data.t .+ 1.0) ./ pdf(data.kd,data.F))
    logL = wLogLikelihood(L,data.t .^(0.7) ./ pdf(data.kd,data.F))
    # logL = wLogLikelihood(L,data.t .^(1))
    return logL
end

function logLikelihood(L)
    logL = -sum(log.(L))
end

function wLogLikelihood(L,w)
    wlogL = -sum(w .* log.(L))
end

I(P) = log(1/P)

using QuadGK
function mean_lifetime(sol) # numerical integration to be changed
    # <t> = ∫p(t) between [0,∞), in this case <t> = ∫sol
    # stop inegration when solution reaches 0
    endpt_id = findfirst(iszero,sol[:])
    if endpt_id == nothing
        endpt = sol.t[end]
    else
        endpt = sol.t[endpt_id]
    end

    t_meanᵤ ,err = quadgk(sol,0.0,endpt,atol=1e-10,rtol=1e-10)
    # if sum(t_meanᵤ) == 0.0
    #     @show sol
    # end
    t_mean = sum(t_meanᵤ)
end

using NumericalIntegration
function mean_lifetime_linear(sol)
    t = sol.t
    u = sum(sol,dims=1)'
    sum(integrate(t,u))

end

function optimization(M,data,options)

    # loss_in = (p) -> mle_loss(p,data,M)
    loss_in = (p) -> mle_loss(10 .^ p,data,M)

    bounds = options
    loss_in(pinit)

    # BBO optimization (use log bounds for easier search across orders of magnitude)
    # #= -----------------------------------------------------------------------
    log_bounds = map(x->log10.(x .+ 1e-5), bounds)
    res_bbo = bboptimize(loss_in; SearchRange = log_bounds,
    # res_bbo = bboptimize(loss_in; SearchRange = bounds,
                                NumDimensions = length(pinit),
                                Method = :adaptive_de_rand_1_bin,
                                # Method = :adaptive_de_rand_1_bin_radiuslimited,
                                # Method = :separable_nes,
                                NThreads=Threads.nthreads()-1,
                                MaxSteps = 1500)#,
    #
    p₀_bbo = best_candidate(res_bbo)
    p₀_bbo = 10 .^ p₀_bbo
    #-----------------------------------------------------------------------
    # =#

    # @show p₀_bbo
    # p₀_bbo = log10.(pinit)
    # p₀_bbo = pinit

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

    # then hone in on minimum with gradient based method
    loss_in = (p) -> mle_loss(p,data,M)

    # bounds = map(x->log10.(x), options)
    lb,ub = boundaries(bounds)

    # @show so.minimize(x ->loss_in(x), pinit,method="BFGS")
    # @show so.minimize(loss_in, pinit,method="BFGS")
    # @show so.newton(loss_in, pinit)#,method="Newton")
    # #=
    od = OnceDifferentiable(loss_in, p₀_bbo; autodiff = :forward);
    res = optimize(od,lb,ub,p₀_bbo,Fminbox(BFGS()),#,
    # res = optimize(loss_in,lb,ub,p₀_bbo,Fminbox(NelderMead()),#
                                        Optim.Options(show_trace=true,
                                        iterations = 50) )#,
                                        # f_tol = 1e-5,
                                        outer_iterations = 1))
    # =#
    # res = optimize(loss_in,p₀_bbo,lb,ub,Fminbox(BFGS()),
    #                                     iterations = 1,
    #                                     Optim.Options(iterations = 10))

    #= -----------------------------------------------------------------------
    optprob = GalacticOptim.OptimizationFunction((x,p) -> loss_in(x),p₀_bbo,GalacticOptim.AutoForwardDiff())
    # optprob = GalacticOptim.OptimizationFunction((x,p) -> loss_in(x))
    # # prob = GalacticOptim.OptimizationProblem(optprob,p₀_bbo)
    # # res = GalacticOptim.solve(prob,NelderMead(),cb = callback)
    #
    prob = GalacticOptim.OptimizationProblem(optprob,p₀_bbo;lb=lb,ub=ub)
    # )#
    res = GalacticOptim.solve(prob,BFGS(initial_stepnorm=0.01),
                                Fminbox(BFGS()),
                                cb = callback,
                                maxiters = 10)
    # res = GalacticOptim.solve(prob,BFGS(),cb = callback,maxiters = 20)
    =#

    rates = res.minimizer
    l = res.minimum

    hes = ForwardDiff.hessian(loss_in,rates)

    @show M.model
    @show l
    @show rates
    return rates, l, hes
end

callback = function (p, l)
  display(l)
  return false
end

function boundaries(opt)
    lb = [b[1] for b in opt]
    ub = [b[2] for b in opt]
    return lb,ub
end

using LinearAlgebra
hessian2σ(H) = sqrt.(abs.( diag(inv(H)) ))


text:
# using Turing
# using Stan
#
# include(Pkg.dir("Turing")*"/benchmarks/benchmarkhelper.jl")
# include(Pkg.dir("Turing")*"/example-models/stan-models/MoC-stan.data.jl")
#
# @model nbmodel(K, V, M, N, z, w, doc, alpha, beta) = begin
#   theta ~ Dirichlet(alpha)
#   phi = Array{Any}(undef, K)
#   for k = 1:K
#     phi[k] ~ Dirichlet(beta)
#   end
#   for m = 1:M
#     z[m] ~ Categorical(theta)
#   end
#   for n = 1:N
#     w[n] ~ Categorical(phi[z[doc[n]]])
#   end
#   phi
# end
#
#
# bench_res = tbenchmark("NUTS(1000, 0.65)", "nbmodel", "data=nbstandata[1]")
# bench_res[4].names = ["phi[1]", "phi[2]", "phi[3]", "phi[4]"]
# logd = build_logd("Naive Bayes", bench_res...)
#
# include(Pkg.dir("Turing")*"/benchmarks/"*"MoC-stan.run.jl")
# logd["stan"] = stan_d
# logd["time_stan"] = nb_time
#
# print_log(logd)

include("../utility.jl")
using Turing
using Random
Random.seed!(129)
# Test for vectorize UnivariateDistribution
@model vdemo(x) = begin
  s ~ InverseGamma(2,3)
  m ~ Normal(0, sqrt(s))
  x ~ [Normal(m, sqrt(s))]
  # for i = 1:length(x)
  #   x[i] ~ Normal(m, sqrt.(s))
  # end
  return s, m
end

alg = HMC(250, 0.01, 5)
x = randn(1000)
res = sample(vdemo(x), alg)

# check_numerical(res, [:s, :m], [1, sum(x) / (1 + length(x))])

# Test for vectorize MultivariateDistribution

D = 2
@model vdemo2(x) = begin
  μ ~ MvNormal(zeros(D), ones(D))
  x ~ [MvNormal(μ, ones(D))]
end

alg = HMC(250, 0.01, 5)
res = sample(vdemo2(randn(D,1000)), alg)

# TODO: Test for vectorize MatrixDistribution


text:
module CTMCEnsemble

using StatsBase

export
    average,
    product,
    vote,
    sink,
    bordacount,
    dempstershafer_rogova,
    powermethod,
    svdmethod,
    ctmc,
    top1,
    top5,
    softmax!

"""
    average(preds, weights=nothing; multiplicity=true)

Compute the average. `multiplicity` for arithmetic mean based on mulitpiliciity.

# Example

```jldoctest
julia> average([([0.5, 0.5], [1, 2]), ([0.5, 0.5], [2, 3])])
3-element Array{Float64,1}:
 0.333333
 0.333333
 0.333333

julia> average([([0.5, 0.5], [1, 2]), ([0.5, 0.5], [2, 3])], multiplicity=false)
3-element Array{Float64,1}:
 0.25
 0.5
 0.25

julia> A = [0.5 0.2 0.1 0.3; 0.5 0.8 0.9 0.7]
2×4 Array{Float64,2}:
 0.5  0.2  0.1  0.3
 0.5  0.8  0.9  0.7

julia> B = [0.5 0.1 0.2 0.7; 0.5 0.9 0.8 0.3]
2×4 Array{Float64,2}:
 0.5  0.1  0.2  0.7
 0.5  0.9  0.8  0.3

julia> average([(A, [1, 2]), (B, [2, 3])])
3×4 Array{Float64,2}:
 0.333333  0.129032  0.0689655  0.230769
 0.333333  0.290323  0.37931    0.538462
 0.333333  0.580645  0.551724   0.230769

julia> average([(A, [1, 2]), (B, [2, 3])], multiplicity=false)
3×4 Array{Float64,2}:
 0.25  0.1   0.05  0.15
 0.5   0.45  0.55  0.7
 0.25  0.45  0.4   0.15
```
"""
function average(preds, weights=nothing; multiplicity=true)
    ndims(preds[1][1]) == 1 &&
        return _average(preds, weights, multiplicity=multiplicity)
    ans = []
    for i = 1:size(preds[1][1], 2)
        push!(ans, _average(map(x -> (x[1][:, i], x[2]), preds), weights, multiplicity=multiplicity))
    end
    hcat(ans...)
end

function _average(preds, weights=nothing; multiplicity=true)
    nclass = maximum(maximum.(getindex.(preds, 2)))
    weights == nothing && (weights = ones(length(preds)))
    v = zeros(nclass)
    m = zeros(nclass)
    for ((pred, label), w) in zip(preds, weights)
        v[label] .+= w .* pred
        m[label] .+= w
    end
    if multiplicity
        assert(all(m .> 0))
        v ./= m
    end
    normalize!(v, 1)
end

"""
    product(preds, weights=nothing; multiplicity=true)

Compute the product. `multiplicity` for geometric mean based on mulitpiliciity.

# Example

```jldoctest
julia> product([([0.5, 0.5], [1, 2]), ([0.5, 0.5], [2, 3])])
3-element Array{Float64,1}:
 0.333333
 0.333333
 0.333333

julia> product([([0.5, 0.5], [1, 2]), ([0.5, 0.5], [2, 3])], multiplicity=false)
3-element Array{Float64,1}:
 0.4
 0.2
 0.4

julia> A = [0.5 0.2 0.1 0.3; 0.5 0.8 0.9 0.7]
2×4 Array{Float64,2}:
 0.5  0.2  0.1  0.3
 0.5  0.8  0.9  0.7

julia> B = [0.5 0.1 0.2 0.7; 0.5 0.9 0.8 0.3]
2×4 Array{Float64,2}:
 0.5  0.1  0.2  0.7
 0.5  0.9  0.8  0.3

julia> product([(A, [1, 2]), (B, [2, 3])])
3×4 Array{Float64,2}:
 0.333333  0.14463   0.0755136  0.230769
 0.333333  0.204537  0.320377   0.538462
 0.333333  0.650833  0.604109   0.230769

julia> product([(A, [1, 2]), (B, [2, 3])], multiplicity=false)
3×4 Array{Float64,2}:
 0.4  0.169492   0.0925926  0.275229
 0.2  0.0677966  0.166667   0.449541
 0.4  0.762712   0.740741   0.275229
```
"""
function product(preds, weights=nothing; multiplicity=true)
    ndims(preds[1][1]) == 1 &&
        return _product(preds, weights, multiplicity=multiplicity)
    ans = []
    for i = 1:size(preds[1][1], 2)
        push!(ans, _product(map(x -> (x[1][:, i], x[2]), preds), weights, multiplicity=multiplicity))
    end
    hcat(ans...)
end

function _product(preds, weights=nothing; multiplicity=true)
    nclass = maximum(maximum.(getindex.(preds, 2)))
    weights == nothing && (weights = ones(length(preds)))
    v = ones(nclass)
    m = zeros(nclass)
    for ((pred, label), w) in zip(preds, weights)
        v[label] .*= w .* pred
        m[label] .+= w
    end
    if multiplicity
        assert(all(m .> 0))
        v .^= (1 ./ m)
    end
    normalize!(v, 1)
end

"""
    dempstershafer_rogova(preds)

Compute the dempstershafer_rogova.

# Example

```jldoctest
julia> dempstershafer_rogova([([0.2, 0.8], [1, 2])])
2-element Array{Float64,1}:
 0.0588235
 0.941176

julia> dempstershafer_rogova([([0.5, 0.5], [1, 2]), ([0.5, 0.5], [2, 3])])
3-element Array{Float64,1}:
 0.428571
 0.142857
 0.428571

julia> A = [0.5 0.2 0.1 0.3; 0.5 0.8 0.9 0.7]
2×4 Array{Float64,2}:
 0.5  0.2  0.1  0.3
 0.5  0.8  0.9  0.7

julia> B = [0.5 0.1 0.2 0.7; 0.5 0.9 0.8 0.3]
2×4 Array{Float64,2}:
 0.5  0.1  0.2  0.7
 0.5  0.9  0.8  0.3

julia> dempstershafer_rogova([(A, [1, 2]), (B, [2, 3])])
3×4 Array{Float64,2}:
 0.428571  0.0503319   0.0134788  0.18598
 0.142857  0.00884956  0.0519897  0.628041
 0.428571  0.940819    0.934531   0.18598
```
"""
function dempstershafer_rogova(preds)
    ndims(preds[1][1]) == 1 &&
        return _dempstershafer_rogova(preds)
    ans = []
    for i = 1:size(preds[1][1], 2)
        push!(ans, _dempstershafer_rogova(map(x -> (x[1][:, i], x[2]), preds)))
    end
    hcat(ans...)
end

function _dempstershafer_rogova(preds)
    nclass = maximum(maximum.(getindex.(preds, 2)))
    v = ones(nclass)
    for (pred, label) in preds
        p = prod(1 .- pred)
        q = pred .* p ./ (1 .- pred)
        e = q ./ (1 .- pred + q)
        v[label] .*= e
    end
    normalize!(v, 1)
end

"""
    vote(preds, weights=nothing)

Compute the vote.

# Example

```jldoctest
julia> vote([([0.5, 0.5], [1, 2]), ([0.5, 0.5], [2, 3])])
3-element Array{Float64,1}:
 0.5
 0.5
 0.0

julia> A = [0.5 0.2 0.1 0.3; 0.5 0.8 0.9 0.7]
2×4 Array{Float64,2}:
 0.5  0.2  0.1  0.3
 0.5  0.8  0.9  0.7

julia> B = [0.5 0.1 0.2 0.7; 0.5 0.9 0.8 0.3]
2×4 Array{Float64,2}:
 0.5  0.1  0.2  0.7
 0.5  0.9  0.8  0.3

julia> vote([(A, [1, 2]), (B, [2, 3])])
3×4 Array{Float64,2}:
 0.5  0.0  0.0  0.0
 0.5  0.5  0.5  1.0
 0.0  0.5  0.5  0.0
```
"""
function vote(preds, weights=nothing)
    ndims(preds[1][1]) == 1 &&
        return _vote(preds, weights)
    ans = []
    for i = 1:size(preds[1][1], 2)
        push!(ans, _vote(map(x -> (x[1][:, i], x[2]), preds), weights))
    end
    hcat(ans...)
end

function _vote(preds, weights=nothing)
    nclass = maximum(maximum.(getindex.(preds, 2)))
    weights == nothing && (weights = ones(length(preds)))
    v = zeros(nclass)
    for ((pred, label), w) in zip(preds, weights)
        i = label[indmax(pred)]
        v[i] .+= w
    end
    normalize!(v, 1)
end

"""
    bordacount(preds, weights=nothing)

Compute the Borda count vote.

# Example

```jldoctest
julia> bordacount([([0.5, 0.5], [1, 2]), ([0.5, 0.5], [2, 3])])
3-element Array{Float64,1}:
 0.25
 0.5
 0.25

julia> A = [0.5 0.2 0.1 0.3; 0.5 0.8 0.9 0.7]
2×4 Array{Float64,2}:
 0.5  0.2  0.1  0.3
 0.5  0.8  0.9  0.7

julia> B = [0.5 0.1 0.2 0.7; 0.5 0.9 0.8 0.3]
2×4 Array{Float64,2}:
 0.5  0.1  0.2  0.7
 0.5  0.9  0.8  0.3

julia> bordacount([(A, [1, 2]), (B, [2, 3])])
3×4 Array{Float64,2}:
 0.25  0.166667  0.166667  0.166667
 0.5   0.5       0.5       0.666667
 0.25  0.333333  0.333333  0.166667
```
"""
function bordacount(preds, weights=nothing)
    ndims(preds[1][1]) == 1 &&
        return _bordacount(preds, weights)
    ans = []
    for i = 1:size(preds[1][1], 2)
        push!(ans, _bordacount(map(x -> (x[1][:, i], x[2]), preds), weights))
    end
    hcat(ans...)
end

function _bordacount(preds, weights=nothing)
    nclass = maximum(maximum.(getindex.(preds, 2)))
    weights == nothing && (weights = ones(length(preds)))
    v = zeros(nclass)
    for ((pred, label), w) in zip(preds, weights)
        v[label] .+= w .* competerank(pred)
    end
    normalize!(v, 1)
end

assemble(p::Vector) = repmat(p, 1, length(p)) - eye(length(p))

"""
    build(preds, weights=nothing; nclass=0)

Build the generator matrix.

# Example

```jldoctest
julia> CTMCEnsemble.build([([0.5, 0.5], [1, 2]), ([0.5, 0.5], [2, 3])])
3×3 Array{Float64,2}:
 -0.5   0.5   0.0
  0.5  -1.0   0.5
  0.0   0.5  -0.5
```
"""
function build(preds, weights=nothing; nclass=0)
    nclass == 0 && (nclass = maximum(maximum.(getindex.(preds, 2))))
    weights == nothing && (weights = ones(length(preds)))
    G = zeros(nclass, nclass)
    for ((pred, label), w) in zip(preds, weights)
        G[label, label] .+= w .* assemble(pred)
    end
    G
end

"""
    powermethod(preds, weights=nothing; maxiter=16)

Compute stationary distribution by power method.

# Example

```jldoctest
julia> powermethod([([0.5, 0.5], [1, 2]), ([0.5, 0.5], [2, 3])], maxiter=16)
3-element Array{Float64,1}:
 0.333333
 0.333333
 0.333333

julia> A = [0.5 0.2 0.1 0.3; 0.5 0.8 0.9 0.7]
2×4 Array{Float64,2}:
 0.5  0.2  0.1  0.3
 0.5  0.8  0.9  0.7

julia> B = [0.5 0.1 0.2 0.7; 0.5 0.9 0.8 0.3]
2×4 Array{Float64,2}:
 0.5  0.1  0.2  0.7
 0.5  0.9  0.8  0.3

julia> powermethod([(A, [1, 2]), (B, [2, 3])])
3×4 Array{Float64,2}:
 0.333333  0.0244257  0.0217406  0.230769
 0.333333  0.0976007  0.195655   0.538462
 0.333333  0.877974   0.782604   0.230769
```
"""
function powermethod(preds, weights=nothing; maxiter=24)
    ndims(preds[1][1]) == 1 &&
        return _powermethod(preds, weights)
    ans = []
    for i = 1:size(preds[1][1], 2)
        push!(ans, _powermethod(map(x -> (x[1][:, i], x[2]), preds), weights, maxiter=maxiter))
    end
    hcat(ans...)
end

function _powermethod(preds, weights=nothing; maxiter=24)
    weights == nothing && (weights = ones(length(preds)))
    A = build(preds, weights)
    A .+= eye(A) .* sum(weights)
    v = average(preds, weights)

    for i in 1:maxiter
        v = normalize!(A * v, 1)
    end

    v / sum(v)
end

"""
    svdmethod(preds, weights=nothing)

Compute stationary distribution by svd method.

# Example

```jldoctest
julia> svdmethod([([0.5, 0.5], [1, 2]), ([0.5, 0.5], [2, 3])])
3-element Array{Float64,1}:
 0.333333
 0.333333
 0.333333

julia> A = [0.5 0.2 0.1 0.3; 0.5 0.8 0.9 0.7]
2×4 Array{Float64,2}:
 0.5  0.2  0.1  0.3
 0.5  0.8  0.9  0.7

julia> B = [0.5 0.1 0.2 0.7; 0.5 0.9 0.8 0.3]
2×4 Array{Float64,2}:
 0.5  0.1  0.2  0.7
 0.5  0.9  0.8  0.3

julia> svdmethod([(A, [1, 2]), (B, [2, 3])])
3×4 Array{Float64,2}:
 0.333333  0.0243902  0.0217391  0.230769
 0.333333  0.097561   0.195652   0.538462
 0.333333  0.878049   0.782609   0.230769
```
"""
function svdmethod(preds, weights=nothing)
    ndims(preds[1][1]) == 1 &&
        return _svdmethod(preds, weights)
    ans = []
    for i = 1:size(preds[1][1], 2)
        push!(ans, _svdmethod(map(x -> (x[1][:, i], x[2]), preds), weights))
    end
    hcat(ans...)
end

function _svdmethod(preds, weights=nothing;)
    λ, ϕ = eig(build(preds, weights))
    v = real(ϕ[:, indmin(abs.(λ))])
    v / sum(v)
end

"""
    stationdist(G)

Compute stationary probability distribution given generator ``G``.

# Example

```jldoctest
julia> G = CTMCEnsemble.build([([0.5, 0.5], [1, 2]), ([0.5, 0.5], [2, 3])]);

julia> CTMCEnsemble.stationdist(G)
3-element Array{Float64,1}:
 0.333333
 0.333333
 0.333333
```
"""
function stationdist(G)
    λ, v = eigs(G, nev=1, which=:LR)
    assert(all(abs.(λ[1]) < 1e-5))
    assert(isreal(v))
    v = vec(real(v))
    v / sum(v)
end

"""
    ctmc(preds, weights=nothing)

Compute stationary distribution by CTMC method.

# Example

```jldoctest
julia> ctmc([([0.5, 0.5], [1, 2]), ([0.5, 0.5], [2, 3])])
3-element Array{Float64,1}:
 0.333333
 0.333333
 0.333333

julia> A = [0.5 0.2 0.1 0.3; 0.5 0.8 0.9 0.7]
2×4 Array{Float64,2}:
 0.5  0.2  0.1  0.3
 0.5  0.8  0.9  0.7

julia> B = [0.5 0.1 0.2 0.7; 0.5 0.9 0.8 0.3]
2×4 Array{Float64,2}:
 0.5  0.1  0.2  0.7
 0.5  0.9  0.8  0.3

julia> ctmc([(A, [1, 2]), (B, [2, 3])])
3×4 Array{Float64,2}:
 0.333333  0.0243902  0.0217391  0.230769
 0.333333  0.097561   0.195652   0.538462
 0.333333  0.878049   0.782609   0.230769
```
"""
function ctmc(preds, weights=nothing)
    ndims(preds[1][1]) == 1 &&
        return stationdist(build(preds, weights))
    ans = []
    for i = 1:size(preds[1][1], 2)
        push!(ans, stationdist(build(map(x -> (x[1][:, i], x[2]), preds), weights)))
    end
    hcat(ans...)
end

"""
    sink(preds, weights=nothing)

Compute stationary distribution by sink method.

# Example

```jldoctest
julia> sink([([0.5, 0.5], [1, 2]), ([0.5, 0.5], [2, 3])])
3-element Array{Float64,1}:
 1.0
 0.0
 0.0

julia> A = [0.5 0.2 0.1 0.3; 0.5 0.8 0.9 0.7]
2×4 Array{Float64,2}:
 0.5  0.2  0.1  0.3
 0.5  0.8  0.9  0.7

julia> B = [0.5 0.1 0.2 0.7; 0.5 0.9 0.8 0.3]
2×4 Array{Float64,2}:
 0.5  0.1  0.2  0.7
 0.5  0.9  0.8  0.3

julia> sink([(A, [1, 2]), (B, [2, 3])])
3×4 Array{Float64,2}:
 1.0  0.0  0.0  0.0
 0.0  0.0  0.0  1.0
 0.0  1.0  1.0  0.0

julia> sink([(A, [1, 2]), (B, [2, 3])], multiple_src=true)
3×4 Array{Float64,2}:
 0.333333  0.0  0.0  0.0
 0.333333  0.0  0.0  1.0
 0.333333  1.0  1.0  0.0
```
"""
function sink(preds, weights=nothing; multiple_src=false)
    ndims(preds[1][1]) == 1 &&
        return chase(build(preds, weights), 1)
    src = 1
    if multiple_src
        nclass = maximum(maximum.(getindex.(preds, 2)))
        src = 1:nclass
    end
    ans = []
    for i = 1:size(preds[1][1], 2)
        push!(ans, chase(build(map(x -> (x[1][:, i], x[2]), preds), weights), src))
    end
    hcat(ans...)
end

"""
    chase(Q::Array, src)

Chase to the highest points from src points in the graph with weights Q.
\$Q_{i,j}\$ is the weight going from point \$j\$ to \$i\$.
Returns an array indicates the number of each point being end points.

# Example

```jldoctest
julia> G = CTMCEnsemble.build([([0.5, 0.5], [1, 2]), ([0.5, 0.5], [2, 3])]);

julia> CTMCEnsemble.chase(G, 1)
3-element Array{Float64,1}:
 1.0
 0.0
 0.0

julia> CTMCEnsemble.chase(G, [1, 2])
3-element Array{Float64,1}:
 0.5
 0.5
 0.0

julia> G = CTMCEnsemble.build([([0.4, 0.6], [1, 2]), ([0.4, 0.6], [2, 3])]);

julia> CTMCEnsemble.chase(G, [1, 2, 3])
3-element Array{Float64,1}:
 0.0
 0.0
 1.0
```
"""
function chase(Q::Array, src)
    nclass = size(Q, 1)
    assert(nclass == size(Q, 2) && maximum(src) <= nclass)
    p = zeros(nclass)
    d = Dict()
    for i in src
        at = i
        inds = []
        while true
            at in inds && break
            if haskey(d, at)
                at = d[at]
                break
            end
            push!(inds, at)
            v = Q[:, at] - Q[at, :]
            maximum(v) <= 1e-7 && break
            at = indmax(v)
        end
        for j in inds
            d[j] = at
        end
        p[at] += 1
    end
    normalize!(p, 1)
end

"""
    top1(data, label)

Top-1 accuracy.

# Example

```jldoctest
julia> p = [0.333333  0.0243902  0.0217391  0.230769
            0.333333  0.097561   0.195652   0.538462
            0.333333  0.878049   0.782609   0.230769];

julia> top1(p, [1, 3, 3, 2])
1.0
```
"""
top1(data, label) = mean(indmax.(view.([data], [:], 1:size(data, 2))) .== vec(label))

"""
    top5(data, label)

Top-5 accuracy.

# Example

```jldoctest
julia> p = [0.49  0.09  0.71  0.07  0.28
            0.73  0.48  0.01  0.96  0.51
            0.87  0.09  0.76  0.63  0.39
            0.37  0.65  0.89  0.31  0.42
            0.6   0.49  0.19  0.21  0.77
            0.56  0.32  0.27  1.0   0.92
            0.5   0.83  0.99  0.4   0.81
            0.34  0.03  0.83  0.07  0.62
            0.93  0.75  0.15  0.37  0.21
            0.25  0.19  0.83  0.69  0.64];

julia> top5(p, 1:5)
0.6
```
"""
top5(data, label) = mean(in.(vec(label), getindex.(sortperm.(view.([data], [:], 1:size(data, 2)), rev=true), [1:5])))

"""
    softmax!(data)

Inplace softmax col-wise.

# Example

```jldoctest
julia> p = [0.49  0.09  0.71  0.07  0.28
            0.73  0.48  0.01  0.96  0.51
            0.87  0.09  0.76  0.63  0.39
            0.37  0.65  0.89  0.31  0.42
            0.6   0.49  0.19  0.21  0.77
            0.56  0.32  0.27  1.0   0.92
            0.5   0.83  0.99  0.4   0.81
            0.34  0.03  0.83  0.07  0.62
            0.93  0.75  0.15  0.37  0.21
            0.25  0.19  0.83  0.69  0.64];

julia> softmax!(p)
10×5 Array{Float64,2}:
 0.0907496  0.0711453  0.109493   0.063526   0.0739451
 0.115365   0.10508    0.0543729  0.154694   0.0930673
 0.132702   0.0711453  0.115107   0.111213   0.0825433
 0.0804877  0.124552   0.131088   0.0807574  0.0850571
 0.101302   0.106136   0.0650961  0.0730723  0.120702
 0.0973297  0.0895435  0.0705178  0.161007   0.140235
 0.0916616  0.149116   0.144874   0.0883626  0.125628
 0.0781089  0.0670021  0.123454   0.063526   0.103889
 0.140908   0.137651   0.0625437  0.0857511  0.0689459
 0.0713862  0.0786277  0.123454   0.11809    0.105988
```
"""
function softmax!(data)
    m = maximum(data)
    data .= exp.(data .- m)
    normalize!.(view.([data], [:], 1:size(data, 2)), 1)
    data
end

end # module


text:
# https://fenicsproject.org/docs/dolfin/latest/python/demos/stokes-iterative/demo_stokes-iterative.py.html

module stokes_iterative

using FenicsPy

# Test for PETSc or Tpetra
@assert has_linear_algebra_backend("PETSc") || has_linear_algebra_backend("Tpetra") "DOLFIN has not been configured with Trilinos or PETSc. Exiting."

@assert has_krylov_solver_preconditioner("amg") "Sorry, this demo is only available when DOLFIN is compiled with AMG  \n preconditioner, Hypre or ML."

if has_krylov_solver_method("minres")
    krylov_method = "minres"
elseif has_krylov_solver_method("tfqmr")
    krylov_method = "tfqmr"
else
    println("Default linear algebra backend was not compiled with MINRES or TFQMR \n Krylov subspace method. Terminating.")
end

# Load mesh
mesh = UnitCubeMesh.create(16, 16, 16, CellType.Type.hexahedron)

# Build function space
P2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
TH = P2 * P1
W = FunctionSpace(mesh, TH)

# Boundaries
right = "x[0] > (1.0 - $DOLFIN_EPS)"
left = "x[0] < $DOLFIN_EPS"
top_bottom = "x[1] > 1.0 - $DOLFIN_EPS || x[1] < $DOLFIN_EPS"

# No-slip boundary condition for velocity
noslip = Constant((0.0, 0.0, 0.0))
bc0 = DirichletBC(W.sub(0), noslip, top_bottom)

# Inflow boundary condition for velocity
inflow = Expression(("-sin(x[1]*pi)", "0.0", "0.0"), degree=2)
bc1 = DirichletBC(W.sub(0), inflow, right)

# Collect boundary conditions
bcs = [bc0, bc1]

# Define variational problem
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
f = Constant((0.0, 0.0, 0.0))
a = inner(grad(u), grad(v))*dx + div(v)*p*dx + q*div(u)*dx
L = inner(f, v)*dx

# Form for use in constructing preconditioner matrix
b = inner(grad(u), grad(v))*dx + p*q*dx

# Assemble system
A, bb = assemble_system(a, L, bcs)

# Assemble preconditioner system
P, btmp = assemble_system(b, L, bcs)

# Create Krylov solver and AMG preconditioner
solver = KrylovSolver(krylov_method, "amg")

# Associate operator (A) and preconditioner matrix (P)
solver.set_operators(A, P)

# Solve
U = FeFunction(W)
solver.solve(U.vector(), bb)

# Get sub-functions
u, p = U.split()

# Save solution in VTK format
ufile_pvd = File("stokes_iterative/velocity.pvd")
ufile_pvd << u
pfile_pvd = File("stokes_iterative/pressure.pvd")
pfile_pvd << p

end # module stokes_iterative


text:
module JacobiFunc

# Jacobi elliptic functions and complete elliptic integral of the first kind.
#
# References:
#      iK: Abramowitz, Stegun: Handbook of Mathematical Functions, Sec. 17.6.
#  jacobi: Abramowitz, Stegun: Handbook of Mathematical Functions, Sec. 16.12 & 16.13.
#          ellipjc.m in Toby Driscoll's Schwarz-Christoffel Toolbox
#

function iK(m::AbstractFloat)
    a,b = one(m),sqrt(m)
    while abs(a-b) > eps(a+b)
        a,b = 0.5*(a+b),sqrt(a*b)
    end
    return π/(2*a)
end
iK(m::Integer) = iK(float(m))
K(m) = iK(one(m)-m)

function landaudescend(m::AbstractFloat)
    f = one(m)
    κ = Vector{typeof(m)}()
    while true
        if m > 1e-3
            push!(κ, (1-sqrt(1-m))/(1+sqrt(1-m)))
        else
            # Improve accuracy by evaluating a Taylor expansion
            # instead of original formula for small m
            push!(κ, @evalpoly m/4 0 1 2 5 14 42 132)
        end
        m = κ[end]^2
        f *= 1+κ[end]
        if m < 4*eps(typeof(m)) break end
    end
    return f,κ
end
landaudescend(m::Integer) = landaudescend(float(m))

function jacobiκ(u::Number,κ)
    sn,cn,dn = sin(u) - κ[end]^2/4 * (u - sin(u)*cos(u)) * cos(u),
               cos(u) - κ[end]^2/4 * (u - sin(u)*cos(u)) * sin(u),
               1 - κ[end]^2/2 * sin(u)^2;
    for i = length(κ):-1:1
        sn,cn,dn = (1+κ[i])*sn/(1+κ[i]*sn^2),
                         cn*dn/(1+κ[i]*sn^2),
                 (1-κ[i]*sn^2)/(1+κ[i]*sn^2)
    end
    return sn,cn,dn
end
jacobiκ(u::AbstractArray, κ) = [jacobiκ(u,κ) for u in u]

function jacobi(u,m)
    f,κ = landaudescend(m)
    return jacobiκ(u/f,κ)
end

end # module


text:
"""
	potree2trie(potree::String) -> trie.children

Trie data structures for Potree hierarchy:
 - each file in octree directory is stored in a node,
 - first node is `r`, root of potree.
"""
function potree2trie(potree::String)
	metadata = CloudMetadata(potree) # metadata of potree
	tree = joinpath(potree,metadata.octreeDir,"r") # path to octree files
	trie = DataStructures.Trie{String}()

	# search all files, si potrebbe fare con un'unica function
	if ( metadata.pointAttributes[1] == "LAS" )
		files = searchfile(tree,".las")
	elseif ( metadata.pointAttributes[1] == "LAZ" )
		files = searchfile(tree,".laz")
	elseif ( metadata.pointAttributes[1] != "LAZ" )	# modo poco elegante di individuare != las/laz
		files = searchfile(tree,".bin")
	else
		throw(DomainError(metadata.pointAttributes,"Format not yet allowed"))
	end

	# build trie from filename
	for file in files
		name = rsplit(splitdir(file)[2],".")[1]
		trie[name] = file
	end

	# root of potree
	return trie.children['r']
end

"""
	max_depth(trie::DataStructures.Trie{String}) -> Int

Maximum depth of trie.
"""
function max_depth(trie::DataStructures.Trie{String})::Int
	if length(trie.children) == 0
		return 0
	else
		depth = Int[]
		for key in collect(keys(trie.children))
			push!(depth,max_depth(trie.children[key]))
		end
		return max(depth...) + 1
	end
end

"""
	cut_trie!(trie::DataStructures.Trie{String}, LOD::Int, l = 0::Int)

Return truncated trie.
"""
function cut_trie!(trie::DataStructures.Trie{String}, LOD::Int, l = 0::Int)
	if l >= LOD
		empty!(trie.children) #delete all children
	end
	for key in collect(keys(trie.children))
		cut_trie!(trie.children[key], LOD, l+1) # next level
	end
end

"""
	get_all_values(trie::DataStructures.Trie{String}) -> Array{String,1}

Return a collection of all values in trie.
"""
function get_all_values(trie::DataStructures.Trie{String})::Array{String,1}
	key_s = collect(keys(trie))
	data = Array{String,1}(undef,length(key_s))
	for i in 1:length(key_s)
		data[i] = trie[key_s[i]]
	end
	return data
end

"""
	sub_trie(t::DataStructures.Trie{String}, root::AbstractString) -> DataStructures.Trie

Return the subtrie with defined `root` node.
"""
function sub_trie(t::DataStructures.Trie{String}, root::AbstractString)::DataStructures.Trie
	return subtrie(t, root[2:end])
end

# =====================
"""
get_files_in_potree_folder(potree::String, LOD::Int, all_prev=true::Bool) -> Array{String,1}

Return all files at that level of potree
"""
function get_files_in_potree_folder(potree::String, LOD::Int, all_prev=true::Bool)::Array{String,1}
	trie = potree2trie(potree)
	return get_files(trie, LOD, String[], 0, all_prev)
end

"""
get_files(trie::DataStructures.Trie{String}, LOD::Int, data=String[]::Array{String,1}, l = 0::Int, all_prev = true::Bool) -> Array{String,1}

Accumulate all values from root to defined level `LOD`.
"""
function get_files(trie::DataStructures.Trie{String}, LOD::Int, data=String[]::Array{String,1}, l = 0::Int, all_prev = true::Bool)::Array{String,1}
	if all_prev
		if l<=LOD
			push!(data,trie.value)
			for key in collect(keys(trie.children))
				get_files(trie.children[key],LOD,data,l+1,all_prev)
			end
		end
	else
		if l==LOD
			push!(data,trie.value)
		end

		if l<LOD
			for key in collect(keys(trie.children))
				get_files(trie.children[key],LOD,data,l+1,all_prev)
			end
		end
	end
	return data
end

"""
readhrc(potree::String)

Reads the hrc file of the potree
"""
function readhrc(potree::String)

	cloudmetadata = CloudMetadata(potree) # useful parameters togli quelli che non usi
	tree = joinpath(potree,cloudmetadata.octreeDir,"r") # path to directory "r"
	hrcs = searchfile(tree,".hrc")

	for hrc in hrcs
		raw = read(hrc)
		treehrc = reshape(raw, (5, div(length(raw), 5)))

		for i in 1:size(treehrc,2)
			children = bitstring(UInt8(treehrc[1,i]))
			npoints = parse(Int, bitstring(UInt8(treehrc[5,i]))*bitstring(UInt8(treehrc[4,i]))*bitstring(UInt8(treehrc[3,i]))*bitstring(UInt8(treehrc[2,i])); base=2)
			#struct da finire
			@show children, npoints
		end
	end

end


text:
import Base.GMP: MPZ
import MixedSubdivisions: MixedCell

struct BinomialSystemSolver
    A::Matrix{Int32} # binomial system lhs
    b::Vector{ComplexF64} # binomial system rhs
    X::ElasticArray{ComplexF64,2,1} # Solutions
    # HNF
    H::Matrix{Int64} # hermite normal form
    U::Matrix{Int64} # trafo matrix for hnf
    H_big::Matrix{BigInt} # hermite normal form
    U_big::Matrix{BigInt} # trafo matrix for hnf
    # Solve data structures
    γ::Vector{Float64} # angle of the elements in b
    μ::Vector{Float64} # γ^U
    αs::Vector{Float64} # partial angles in the triangular solve
    Aᵀ::Matrix{Float64} # transposed of A to solve coords
    unit_roots_table::ElasticArray{Int32,2,1}
end

function BinomialSystemSolver(A::Matrix, b::Vector)
    BSS = BinomialSystemSolver(size(A, 1))
    BSS.A .= A
    BSS.b .= b
    BSS
end

function BinomialSystemSolver(n::Int)
    A = zeros(Int32, n, n)
    b = zeros(ComplexF64, n)
    X = ElasticArray{ComplexF64}(undef, n, 0)
    H = zeros(Int64, n, n)
    U = zeros(Int64, n, n)
    H_big = Matrix{BigInt}(undef, n, n)
    H_big .= zero.(BigInt)
    U_big = Matrix{BigInt}(undef, n, n)
    U_big .= zero.(BigInt)
    γ = zeros(Float64, n)
    μ = zeros(Float64, n)
    αs = zeros(Float64, n)
    Aᵀ = zeros(Float64, n, n)
    S = ElasticArray{Int32}(undef, n, 0)
    BinomialSystemSolver(A, b, X, H, U, H_big, U_big, γ, μ, αs, Aᵀ, S)
end

function init!(BSS::BinomialSystemSolver, support, coeffs, cell::MixedCell)
    for (i, (aᵢ, bᵢ)) in enumerate(cell.indices)
        for j = 1:size(BSS.A, 1)
            BSS.A[j, i] = support[i][j, aᵢ] - support[i][j, bᵢ]
        end
        BSS.b[i] = -coeffs[i][bᵢ] / coeffs[i][aᵢ]
    end
    BSS
end

function fill_unit_roots_combinations!(S::ElasticArray, H, d̂)
    n = size(H, 1)
    resize!(S, n, d̂)
    d, e = d̂, 1
    @inbounds for i = 1:n
        dᵢ = Int64(H[i, i])
        d = d ÷ dᵢ
        k = 1
        for _ = 1:e, j = 0:(dᵢ-1), _ = 1:d
            S[i, k] = j
            k += 1
        end
        e *= dᵢ
    end
    S
end

function compute_angular_part!(
    BSS::BinomialSystemSolver,
    H::Matrix{Int64},
    U::Matrix{Int64},
    d̂::Int,
)
    @unpack X, γ, b, μ, unit_roots_table, αs = BSS
    n = size(H, 1)

    # apply coordinate change to b
    γ .= angle.(b) ./ 2π
    μ .= zero(eltype(μ))
    for j = 1:n
        for i = 1:n
            μ[j] += rem(U[i, j] * γ[i], 2.0, RoundNearest)
        end
        μ[j] = rem(μ[j], 2.0, RoundNearest)
    end
    # solve triangular system
    @inbounds for i = 1:d̂, j = n:-1:1
        α = μ[j] + unit_roots_table[j, i]
        α = rem((μ[j] + unit_roots_table[j, i]) / H[j, j], 2.0, RoundNearest)
        for k = n:-1:(j+1)
            α -= rem(αs[k] * H[k, j] / H[j, j], 2.0, RoundNearest)
            α = rem(α, 2.0, RoundNearest)
        end
        # α = rem(α, 2.0, RoundNearest)
        X[j, i] = cis(2π * α)
        αs[j] = α
    end
    BSS
end

function compute_angular_part!(
    BSS::BinomialSystemSolver,
    H::Matrix{BigInt},
    U::Matrix{BigInt},
    d̂::Int,
)
    # @unpack X, γ, b, μ, unit_roots_table = BSS
    @unpack X, unit_roots_table, b = BSS
    n = size(H, 1)

    # compute precision necessary
    p = max(maximum(x -> MPZ.sizeinbase(x, 2), H), maximum(x -> MPZ.sizeinbase(x, 2), U))
    prec = max((ceil(Int, (p + 53) / 32)) * 32, 64)

    γ = angle.(b) ./ 2π

    μ = [BigFloat(0.0; precision = prec) for i = 1:n]
    α = γᵢ = BigFloat(0.0; precision = prec)
    αk = γᵢⱼ = BigFloat(0.0; precision = prec)
    m = BigFloat(2.0; precision = prec)

    for i = 1:n
        set!(γᵢ, γ[i])
        for j = 1:n
            mul!(γᵢⱼ, γᵢ, U[i, j])
            if γᵢⱼ < -1 || γᵢⱼ > 1
                rem!(γᵢⱼ, γᵢⱼ, m, RoundNearest)
            end
            add!(μ[j], μ[j], γᵢⱼ)
        end
    end
    αs = [BigFloat(0.0; precision = prec) for i = 1:n]
    # solve triangular system
    for i = 1:d̂, j = n:-1:1
        add!(α, μ[j], Int64(unit_roots_table[j, i]))
        div!(α, α, H[j, j])
        for k = n:-1:(j+1)
            mul!(αk, αs[k], H[k, j])
            div!(αk, αk, H[j, j])
            sub!(α, α, αk)
        end
        if α < -1 || α > 1
            rem!(α, α, m, RoundNearest)
        end
        X[j, i] = cis(2π * Float64(α))
        set!(αs[j], α)
    end
    BSS
end

function solve(BSS::BinomialSystemSolver, support, coeffs, cell::MixedCell)
    init!(BSS, support, coeffs, cell)
    solve!(BSS)
end

function solve(BSS::BinomialSystemSolver, A, b)
    BSS.A .= A
    BSS.b .= b
    solve!(BSS)
end

function solve!(BSS::BinomialSystemSolver)
    try
        hnf!(BSS.H, BSS.U, BSS.A)
        solve!(BSS, BSS.H, BSS.U)
    catch e
        isa(e, OverflowError) || rethrow(e)
        hnf!(BSS.H_big, BSS.U_big, BSS.A)
        solve!(BSS, BSS.H_big, BSS.U_big)
    end
    BSS.X
end

function solve!(BSS::BinomialSystemSolver, H::Matrix{<:Integer}, U::Matrix{<:Integer})
    n = size(H, 1)
    d̂ = 1
    for i = 1:n
        d̂ *= Int64(H[i, i])
    end
    resize!(BSS.X, n, d̂)
    fill_unit_roots_combinations!(BSS.unit_roots_table, H, d̂)

    # We split the computation of the solution in 2 stages
    # 1) Compute solutions for the angular part -> d̂ many solutions
    compute_angular_part!(BSS, H, U, d̂)

    @unpack A, Aᵀ, b, μ, X = BSS
    # 2) Solve for the absolute value
    μ .= log.(fast_abs.(b))
    LA.transpose!(Aᵀ, A)
    LA.ldiv!(LA.lu!(Aᵀ), μ)
    μ .= exp.(μ)
    for j = 1:d̂, i = 1:n
        X[i, j] *= μ[i]
    end
    X
end

#########################
## Hermite Normal Form ##
#########################

"""
    hnf(A, T=elytpe(A))
Compute the hermite normal form `H` of `A` by overwriting `A` and the corresponding transformation
matrix `U` using precision T. This is `A*U == H` and `H` is a lower triangular matrix.
The implementation follows the algorithm described in [1].

[1] Kannan, Ravindran, and Achim Bachem. "Polynomial algorithms for computing the Smith and Hermite normal forms of an integer matrix." SIAM Journal on Computing 8.4 (1979): 499-507.
"""
function hnf(A, T = eltype(A))
    H = similar(A, T)
    H .= A
    U = similar(A, T)
    U .= zero.(T)
    hnf!(H, U)
    H, U
end

# use checked arithmethic
⊡(x, y) = Base.checked_mul(x, y)
⊞(x, y) = Base.checked_add(x, y)

"""
    hnf!(H, U, A)
Inplace version of [hnf](@ref).

    hnf!(A, U)
Inplace version of [hnf](@ref) overwriting `A` with `H`.
"""
hnf!(H, U, A) = hnf!(copyto!(H, A), U)
function hnf!(A, U)
    n = size(A, 1)
    U .= 0
    @inbounds for i = 1:n
        U[i, i] = one(eltype(U))
    end
    @inbounds for i = 1:(n-1)
        ii = i ⊞ 1
        for j = 1:i
            if !iszero(A[j, j]) || !iszero(A[j, ii])
                # 4.1
                r, p, q = gcdx(A[j, j], A[j, ii])
                # 4.2
                d_j = -A[j, ii] ÷ r
                d_ii = A[j, j] ÷ r
                for k = 1:n
                    a_kj, a_kii = A[k, j], A[k, ii]
                    A[k, j] = a_kj ⊡ p ⊞ a_kii ⊡ q
                    A[k, ii] = a_kj ⊡ d_j ⊞ a_kii ⊡ d_ii

                    u_kj, u_kii = U[k, j], U[k, ii]
                    U[k, j] = u_kj ⊡ p ⊞ u_kii ⊡ q
                    U[k, ii] = u_kj ⊡ d_j ⊞ u_kii ⊡ d_ii
                end
            end
            # 4.3
            if j > 1
                reduce_off_diagonal!(A, U, j)
            end
        end
        #5
        reduce_off_diagonal!(A, U, ii)
    end
    # Special case for 1 × 1 matrices to guarantee that the diagonal is positive
    # This comes up in polyhedral homotopy for univariate polynomials
    if n == 1
        U[1, 1] = flipsign(U[1, 1], A[1, 1])
        A[1, 1] = flipsign(A[1, 1], A[1, 1])
    end

    nothing
end

function cdiv(x, y)
    d = ceil(typeof(x), x / y)
    if 0 ≥ x - d * y > -y
        return d
    else
        d, r = divrem(x, y)
        d + (r > 0)
    end
end

@inline function reduce_off_diagonal!(A, U, k)
    n = size(A, 1)
    @inbounds if A[k, k] < 0
        for i = 1:n
            A[i, k] = -A[i, k]
            U[i, k] = -U[i, k]
        end
    end
    @inbounds for z = 1:(k-1)
        if !iszero(A[z, z])
            # r = -ceil(eltype(A), big(A[k, z]) / A[z, z])
            r = -cdiv(A[k, z], A[z, z])
            for i = 1:n
                U[i, z] = U[i, z] ⊞ r ⊡ U[i, k]
                A[i, z] = A[i, z] ⊞ r ⊡ A[i, k]
            end
        end
    end
    nothing
end

function hnf!(A::AbstractMatrix{BigInt}, U::AbstractMatrix{BigInt})
    n = size(A, 1)
    for j = 1:n, i = 1:n
        if i == j
            MPZ.set_si!(U[i, j], 1)
        else
            MPZ.set_si!(U[i, j], 0)
        end
    end
    r, p, q = zero(BigInt), zero(BigInt), zero(BigInt)
    d_j, d_ii = zero(BigInt), zero(BigInt)
    t1, t2, t3 = zero(BigInt), zero(BigInt), zero(BigInt)
    @inbounds for i = 1:(n-1)
        ii = i + 1
        for j = 1:i
            if !iszero(A[j, j]) || !iszero(A[j, ii])
                # 4.1
                MPZ.gcdext!(r, p, q, A[j, j], A[j, ii])
                # 4.2
                # d_j = -A[j, ii] ÷ r
                MPZ.mul_si!(d_j, A[j, ii], -1)
                MPZ.tdiv_q!(d_j, d_j, r)

                # d_ii = A[j, j] ÷ r
                MPZ.tdiv_q!(d_ii, A[j, j], r)
                for k = 1:n
                    # a_kj, a_kii = A[k, j], A[k, ii]
                    a_kj = MPZ.set!(t1, A[k, j])
                    a_kii = MPZ.set!(t2, A[k, ii])
                    # A[k, j] = a_kj ⊡ p ⊞ a_kii ⊡ q
                    MPZ.mul!(A[k, j], a_kj, p)
                    MPZ.mul!(t3, a_kii, q)
                    MPZ.add!(A[k, j], A[k, j], t3)
                    # A[k, ii] = a_kj ⊡ d_j ⊞ a_kii ⊡ d_ii
                    MPZ.mul!(A[k, ii], a_kj, d_j)
                    MPZ.mul!(t3, a_kii, d_ii)
                    MPZ.add!(A[k, ii], A[k, ii], t3)

                    # u_kj, u_kii = U[k, j], U[k, ii]
                    u_kj = MPZ.set!(t1, U[k, j])
                    u_kii = MPZ.set!(t2, U[k, ii])
                    # U[k, j] = u_kj ⊡ p ⊞ u_kii ⊡ q
                    MPZ.mul!(U[k, j], u_kj, p)
                    MPZ.mul!(t3, u_kii, q)
                    MPZ.add!(U[k, j], U[k, j], t3)

                    # U[k, ii] = u_kj ⊡ d_j ⊞ u_kii ⊡ d_ii
                    MPZ.mul!(U[k, ii], u_kj, d_j)
                    MPZ.mul!(t3, u_kii, d_ii)
                    MPZ.add!(U[k, ii], U[k, ii], t3)
                end
            end
            # 4.3
            if j > 1
                reduce_off_diagonal!(A, U, j, t1, t2)
            end
        end
        #5
        reduce_off_diagonal!(A, U, ii, t1, t2)
    end
    # Special case for 1 × 1 matrices to guarantee that the diagonal is positive
    # This comes up in polyhedral homotopy for univariate polynomials
    if n == 1
        MPZ.mul_si!(U[1, 1], A[1, 1] < 0 ? -1 : 1)
        MPZ.mul_si!(A[1, 1], A[1, 1] < 0 ? -1 : 1)
    end

    nothing
end

cdiv_q!(x::BigInt, a::BigInt, b::BigInt) =
    (ccall((:__gmpz_cdiv_q, :libgmp), Cvoid, (MPZ.mpz_t, MPZ.mpz_t, MPZ.mpz_t), x, a, b);
    x)

function reduce_off_diagonal!(
    A::AbstractMatrix{BigInt},
    U::AbstractMatrix{BigInt},
    k,
    t1,
    t2,
)
    n = size(A, 1)
    @inbounds if A[k, k] < 0
        for i = 1:n
            MPZ.mul_si!(A[i, k], -1)
            MPZ.mul_si!(U[i, k], -1)
        end
    end

    @inbounds for z = 1:(k-1)
        if !iszero(A[z, z])
            r = cdiv_q!(t1, A[k, z], A[z, z])
            r = MPZ.mul_si!(r, r, -1)
            # r = -divceil(A[k, z], A[z,z])
            for i = 1:n
                MPZ.mul!(t2, r, U[i, k])
                MPZ.add!(U[i, z], U[i, z], t2)
                # U[i, z] = U[i, z] ⊞ r ⊡ U[i, k]
                MPZ.mul!(t2, r, A[i, k])
                MPZ.add!(A[i, z], A[i, z], t2)
                # A[i, z] = A[i, z] ⊞ r ⊡ A[i, k]
            end
        end
    end
    nothing
end


text:
using DelimitedFiles

function read_adiabatic_test_data(fi::String, have_low_cut::Bool=false)

    d = readdlm(fi)

    if have_low_cut
        Nbins = length(d[1:end-3,1])
        pmin = d[1,1]
        pmax = d[1+Nbins, 1]

        low_cut  = d[end-1,1]
        CRpC     = d[end,1]

        CRpS = d[1:Nbins,2]
        CRpN = d[1:Nbins,3]
    else
        Nbins = length(d[1:end-2,1])
        pmin = d[1,1]
        pmax = d[1+Nbins, 1]

        low_cut  = pmin
        CRpC     = d[end,1]

        CRpS = d[1:Nbins,2]
        CRpN = d[1:Nbins,3]
    end

    par = CRMomentumDistributionConfig(pmin, pmax, Nbins)#, true)

    cr = CRMomentumDistribution(par.Nbins)

    # get zeroth bin
    cr.CRp_bound[1] = low_cut
    cr.CRp_dis[1] = CRpN[1]

    # all other bins
    j = 2
    for i = 1:Nbins-1

        # upper boundary of bin
        cr.CRp_bound[j] = pmin * 10.0^(par.bin_width*i)
        cr.CRp_dis[j] = CRpN[i] * ( cr.CRp_bound[j]/cr.CRp_bound[j-1])^(-CRpS[i])
        if cr.CRp_bound[j] > CRpC/par.mc_p
            cr.CRp_bound[j] = CRpC/par.mc_p
        end

        # lower bound of next bin
        cr.CRp_bound[j+1] = cr.CRp_bound[j]
        cr.CRp_dis[j+1] = CRpN[i+1]
        if cr.CRp_bound[j] == CRpC/par.mc_p
            cr.CRp_bound[j+1] = CRpC/par.mc_p
        end

        j += 2

    end

    # last boundary
    cr.CRp_bound[j] = pmax
    if cr.CRp_bound[j-1] < CRpC/par.mc_p
        cr.CRp_bound[j] = CRpC/par.mc_p
    end
    cr.CRp_dis[j] = CRpN[Nbins] * ( cr.CRp_bound[j]/cr.CRp_bound[j-1])^(-CRpS[Nbins])
    cr.CRp_bound[j+1] = cr.CRp_bound[j]

    return cr
end


text:
abstract type AbstractObject{T <: SSDFloat} end

@inline function in(pt::AbstractCoordinatePoint{T}, c::AbstractObject{T})::Bool where {T}
    return in(pt, c.geometry)
end

# function in(p::AbstractCoordinatePoint{T}, c::AbstractObject{T}, rs::Vector{T}) where {T <: SSDFloat}
#     rv = false
#     for g in c.geometry_positive
#         if typeof(g) in []#[ConeMantle{T}]
#             in(p, g, rs) ? rv = true : nothing
#         else
#             (p in g) ? rv = true : nothing
#         end
#     end
#     return rv
# end

function in(p::AbstractCoordinatePoint{T}, v::AbstractVector{<:AbstractObject{T}}) where {T <: SSDFloat}
    rv = false
    for object in v
        if p in object
            rv = true
        end
    end
    return rv
end

function find_closest_gridpoint(point::CylindricalPoint{T}, grid::CylindricalGrid{T}) where {T <: SSDFloat}
    return Int[searchsortednearest( grid.r.ticks, point.r), searchsortednearest(grid.φ.ticks, point.φ), searchsortednearest(grid.z.ticks, point.z)]
end
function find_closest_gridpoint(point::CartesianPoint{T}, grid::CylindricalGrid{T}) where {T <: SSDFloat}
    find_closest_gridpoint(CylindricalPoint(point),grid)
end

function find_closest_gridpoint(point::CartesianPoint{T}, grid::CartesianGrid{T}) where {T <: SSDFloat}
    return Int[searchsortednearest( grid[:x].ticks, point.x), searchsortednearest(grid[:y].ticks, point.y), searchsortednearest(grid.z.ticks, point.z)]
end
function find_closest_gridpoint(point::CylindricalPoint{T}, grid::CartesianGrid{T}) where {T <: SSDFloat}
    find_closest_gridpoint(CartesianPoint(point),grid)
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

@recipe function f(solution::EstimationSolution, vars=nothing)
  # retrieve underlying domain
  sdomain = solution.domain

  # valid variables
  validvars = sort(collect(keys(solution.mean)))

  # plot all variables by default
  isnothing(vars) && (vars = validvars)
  @assert vars ⊆ validvars "invalid variable name"

  # plot mean and variance for each variable
  layout --> (length(vars), 2)
  legend --> false

  for (i, var) in enumerate(vars)
    @series begin
      subplot := 2i - 1
      title --> string(var, " mean")
      sdomain, solution.mean[var]
    end
    @series begin
      subplot := 2i
      title --> string(var, " variance")
      sdomain, solution.variance[var]
    end
  end
end


text:
#!/usr/local/bin/julia

using JuMP

include("Types.jl")
include("recipe_protos_cheap.jl")

macro minsec(n) :((60($n), ($n), 60 * $n * mins)) end
macro PerSec(R)
	return quote
		r = get(Recipes, $R, nothing)
		if r == nothing
			1
		else
			Recipes[$R].factory[end:end][1].speed / Recipes[$R].time
		end
	end
end
macro Recipes() :(collect(keys(Recipes))) end
macro FOUT(R, I) :(Factories[$R] * get(Recipes[$R].outs, $I, 0) * @PerSec($R)) end
macro FIN(R, I) :(Factories[$R] * get(Recipes[$R].ins, $I, 0) * @PerSec($R))  end
macro FOUTS(K) :(sum([@FOUT(F, $K) for F in @Recipes])) end
macro FINS(K) :(sum([@FIN(F, $K) for F in @Recipes])) end

macro FCON(R, I) :(getvalue(Factories[$R]) * get(Recipes[$R].ins, $I, 0)) end
macro FCONS(K) :(sum([@FCON(F, $K) for F in @Recipes])) end

function production(r::Recipe, nfacs, mins)
	if nfacs == 0 return end

	function prnt(igs)
		for (k,v) in igs
			@printf "\t\t \"%s\" %0.2f / %0.2f per min/sec : Total %0.2f in \n" k @minsec(@PerSec(k) * v * nfacs)...
		end
	end
	@printf "%s: %0.2f factories in %d mins\n" r.name nfacs mins
	@printf "\tProduces\n"
	prnt(r.outs)
	@printf "\tConsumes\n"
	prnt(r.ins)
end

function report(status, mins)
	if status == :Optimal
		for r in @Recipes
			production(Recipes[r], getvalue(Factories[r]), mins)
		end
		function prnt(i)
			n = @FCONS(i)
			if n > 0
				@printf "\"%s\" %0.2f / %0.2f per min/sec : Total %0.2f\n" i @minsec(n)...
			end
		end
		@printf "Raws in %d mins\n" mins
		for k in ["iron-plate" "copper-plate" "water" "coal" "petroleum-gas" "stone"]
			prnt(k)
		end
	elseif status == :Unbounded
		print("Problem is unbounded")
	elseif status == :Infeasible
		print("Problem is infeasible")
	elseif status == :UserLimit
		print("Iteration limit or timeout")
	elseif status == :Error
		print("Solver exited with an error")
	elseif status == :NotSolved
		print("Model built in memory but not optimized")
	end
end

m = Model()
@variable(m, Factories[@Recipes] >= 0)
for R in @Recipes
	@constraint(m, @FOUTS(R) >= @FINS(R))
end

@objective(m, Min, sum([Factories[R] for R in @Recipes]))

const PerMin = 1 / 60
const Minutes = 1

const YellowBelt = 800PerMin
const RedBelt = 1600PerMin
const BlueBelt = 2400PerMin

@constraint(m, @FOUTS("rocket-part") >= RedBelt / Minutes)

report(solve(m), Minutes)

max_stars_repo_name: lawless-m/factorio-utils

max_issues_repo_name: lawless-m/factorio-utils

max_forks_repo_name: lawless-m/factorio-utils


text:
module PolygonClipping

import Base.show
import Base.push!
import Base.length
using StaticArrays

export Vertex, Polygon, push!, intersection, isinside, show, unprocessed,
       VertexException, EdgeException, DegeneracyException, length, remove, infill

mutable struct Vertex
    location::SVector{2,Float64}
    next::Union{Vertex, Nothing}
    prev::Union{Vertex, Nothing}
    nextpoly::Union{Vertex, Nothing}
    intersect::Bool
    entry::Bool
    neighbor::Union{Vertex, Nothing}
    visited::Bool
    alpha::Float64

end

Vertex(x) = Vertex(x, nothing, nothing, nothing, false, true, nothing, false, 0.0)
Vertex(x::Real,y::Real) = Vertex(SVector(x,y), nothing, nothing, nothing, false, true, nothing, false, 0.0)
Vertex(x::AbstractVector, a::Vertex, b::Vertex) = Vertex(x, a, b, nothing, false, true, nothing, false, 0.0)

struct VertexException <: Exception end
struct EdgeException <: Exception end
struct DegeneracyException <: Exception end

mutable struct Polygon
    start::Union{Vertex, Nothing}

    Polygon() = new(nothing)
end

Base.iterate(m::Polygon) = m.start === nothing ? nothing : (m.start, m.start) # (initial, first?)
function Base.iterate(m::Polygon, state)
    if m.start === state.next
        return nothing
    else
        return (state.next, state.next)
    end
end

function length(p::Polygon)
    n = 0
    for i in p
        n += 1
    end
    n
end

function push!(p::Polygon, v::Vertex)
    if p.start == nothing
        p.start = v
        v.prev = v
        v.next = v
    else
        v.next = p.start
        v.prev = p.start.prev
        p.start.prev.next = v
        p.start.prev = v
    end
end

function Vertex(s::Vertex, c::Vertex, alpha::Float64)
    # Insert a vertex between s and c at alpha from s
    location = s.location + alpha*(c.location-s.location)
    a = Vertex(location, c, s)
    a.alpha = alpha
    s.next = a
    c.prev = a
    return a
end

function Vertex(s::Vertex, c::Vertex, location::SVector{2,Float64})
    # Insert a vertex between s and c at location
    a = Vertex(location, c, s)
    s.next = a
    c.prev = a
    return a
end

function unprocessed(p::Polygon)
    for v in p
        if !v.visited && v.intersect
            return true
        end
    end
    return false
end

function remove(v::Vertex, poly::Polygon)
    if v === poly.start
        if v.next === v
            poly.start = nothing
            v.next = nothing
            v.prev = nothing
            return
        else
            poly.start = v.next
        end
    end
    v.next.prev = v.prev
    v.prev.next = v.next
    v.next = nothing
    v.prev = nothing
    return
end

function detq(q1,q2,r)
    (q1[1]-r[1])*(q2[2]-r[2])-(q2[1]-r[1])*(q1[2]-r[2])
end

function isinside(v::Vertex, poly::Polygon)
    # See: http://www.sciencedirect.com/science/article/pii/S0925772101000128
    # "The point in polygon problem for arbitrary polygons"
    # An implementation of Hormann-Agathos (2001) Point in Polygon algorithm
    c = false
    r = v.location
    for q1 in poly
        q2 = q1.next
        if q1.location == r
            throw(VertexException())
        end
        if q2.location[2] == r[2]
            if q2.location[1] == r[1]
                throw(VertexException())
            elseif (q1.location[2] == r[2]) && ((q2.location[1] > r[1]) == (q1.location[1] < r[1]))
                throw(EdgeException())
            end
        end
        if (q1.location[2] < r[2]) != (q2.location[2] < r[2]) # crossing
            if q1.location[1] >= r[1]
                if q2.location[1] > r[1]
                    c = !c
                elseif ((detq(q1.location,q2.location,r) > 0) == (q2.location[2] > q1.location[2])) # right crossing
                    c = !c
                end
            elseif q2.location[1] > r[1]
                if ((detq(q1.location,q2.location,r) > 0) == (q2.location[2] > q1.location[2])) # right crossing
                    c = !c
                end
            end
        end
    end
    return c
end

function phase1!(p1,p2,p3)
    phase1!(p1,p2)
    phase1!(p1,p3)
end

function phase1!(subject::Polygon, clip::Polygon)
    sv = subject.start
    svn = sv.next
    while true
        cv = clip.start
        cvn = cv
        while true
            # Skip ahead to next non-inserted vertex
            while true
                cvn = cvn.next
                if !cvn.intersect
                    break
                end
            end

            intersect, a, b = intersection(sv, svn, cv, cvn)
            if intersect
                # Find where to insert vertices
                av = sv
                bv = cv
                while av.alpha <= a && av !== svn
                    av = av.next
                end
                while bv.alpha <= b && bv !== cvn
                    bv = bv.next
                end

                location = sv.location + a*(svn.location-sv.location)
                i1 = Vertex(av.prev, av, location)
                i2 = Vertex(bv.prev, bv, location)
                i1.alpha = a
                i2.alpha = b
                i1.intersect = true
                i2.intersect = true
                i1.neighbor = i2
                i2.neighbor = i1
            end
            if cvn === clip.start
                break
            end
            cv = cvn
        end
        if svn === subject.start
            break
        end
        sv = svn
        # Skip ahead to next non-inserted vertex
        while true
            svn = svn.next
            if !svn.intersect
                break
            end
        end
    end
end

function phase2!(subject::Polygon, clip::Polygon)
    status = false
    sv = subject.start
    if isinside(sv, clip)
        status = false
    else
        status = true
    end
    for sv in subject
        if sv.intersect
            sv.entry = status
            status = !status
        else
            sv.entry = status
        end
    end
end

function intersection(subject::Polygon, clip::Polygon)
    phase1!(subject, clip)

    phase2!(subject, clip)
    phase2!(clip, subject)

    results = Polygon[]
    numpoly = 1
    search_start = subject.start
    while true
        current = search_start
        while true
            if current.intersect
                push!(results, Polygon())
                push!(results[numpoly], Vertex(current.location))
                break
            end
            current = current.next
            if current == subject.start
                return results
            end
        end
        start = current
        first_change = true
        while true
            if current.entry
                while true
                    current = current.next
                    push!(results[numpoly], Vertex(current.location))
                    if current.intersect
                        current.intersect = false
                        break
                    end
                end
            else
                while true
                    current = current.prev
                    push!(results[numpoly], Vertex(current.location))
                    if current.intersect
                        current.intersect = false
                        break
                    end
                end
            end
            if current.neighbor == start
                break
            elseif first_change
                first_change = false
                search_start = current
            end
            current = current.neighbor
            current.intersect = false
        end
        numpoly = numpoly + 1
    end
    return results
end

function infill(subject::Polygon, clip::Polygon)
    if subject.start == nothing || clip.start == nothing
        return
    end
    phase1!(subject, clip)

    phase2!(subject, clip)
    phase2!(clip, subject)

    results = Polygon[]
    numpoly = 1
    while unprocessed(subject)
        current = subject.start
        nonintersections = 1
        while true
            if current.intersect
                if !current.visited
                    current.visited = true
                    push!(results, Polygon())
                    newvert = Vertex(current.location)
                    newvert.intersect = true
                    push!(results[numpoly], newvert)
                    #current = current.neighbor
                    break
                end
            else
                nonintersections += 1
            end
            current.visited = true
            current = current.next
        end
        # we are going down the right-hand side of an intersection so our jogs are flipped
        flipped = false
        if nonintersections % 3 == 0
            flipped = true
        end
        start = current.location
        crossings = 1 # count first intersection
        traversing = true
        while traversing
            if current.entry
                while true
                    current = current.next
                    newvert = Vertex(current.location)
                    push!(results[numpoly], newvert)
                    current.visited = true
                    if current.intersect
                        newvert.intersect = true
                        break
                    end
                end
            else
                while true
                    current = current.prev
                    newvert = Vertex(current.location)
                    push!(results[numpoly], newvert)
                    current.visited = true
                    if current.intersect
                        newvert.intersect = true
                        break
                    end
                end
            end
            if current.neighbor.visited && current.intersect
                traversing = false
                # remove vertices that we inserted while traversing
                remove(results[numpoly].start.prev, results[numpoly])
                vert = results[numpoly].start.prev
                while !vert.intersect
                    vert_prev = vert.prev
                    remove(vert, results[numpoly])
                    vert = vert_prev
                end
                break
            end
            current = current.neighbor
            current.visited = true
            crossings = crossings + 1
            if flipped && crossings == 2
                current.entry = !current.entry
            end
            if !flipped && crossings == 4
                current.entry = !current.entry
                crossings = 0
            end
        end
        numpoly = numpoly + 1
    end
    return results
end

function intersection(sv::Vertex, svn::Vertex, cv::Vertex, cvn::Vertex)
    s1 = sv.location
    s2 = svn.location
    c1 = cv.location
    c2 = cvn.location

    den = (c2[2] - c1[2]) * (s2[1] - s1[1]) - (c2[1] - c1[1]) * (s2[2] - s1[2])
    if den == 0.0
        return false, 0.0, 0.0
    end
    a = ((c2[1] - c1[1]) * (s1[2] - c1[2]) - (c2[2] - c1[2]) * (s1[1] - c1[1])) / den
    b = ((s2[1] - s1[1]) * (s1[2] - c1[2]) - (s2[2] - s1[2]) * (s1[1] - c1[1])) / den

    if (0.0 < a < 1.0) && (0.0 < b < 1.0)
        return true, a, b
    elseif ((a in [0.0,1.0]) && (0.0 <= b <= 1.0)) || ((b in [0.0,1.0]) && (0.0 <= a <= 1.0))
        throw(DegeneracyException())
    else
        return false, 0.0, 0.0
    end
end

end # module


text:
mironov = function(f,mv,cf)

    eps_0 = 8.854*10^(-12)
    eps_winf = 4.9
    fHz = f

    #Initializing the GRMDM spectroscopic parameters with clay fraction
    #RI & NAC of dry soils
    znd = 1.634 - 0.539 * cf + 0.2748 * cf^2
    zkd = 0.03952 - 0.04038 * cf

    #Maximum bound water fraction
    zxmvt = 0.02863 + 0.30673 * cf
    #zxmvt = 0.02863 + 0.0030673 * cf

    #Bound water parameters
    zep0b = 79.8 - 85.4 * cf + 32.7 * cf^2
    ztaub = 1.062*10^(-11) + 3.450*10^(-12) * cf
    zsigmab = 0.3112 + 0.467 * cf

    #Unbound (free) water parameters
    zep0u = 100
    ztauu = 8.5*10^(-12)
    zsigmau = 0.3631 + 1.217 * cf

    #Computation of epsilon water (bound & unbound)
    zcxb = (zep0b - eps_winf) / (1 + (2*pi*fHz*ztaub)^2)
    zepwbx = eps_winf + zcxb
    zepwby = zcxb * (2*pi*fHz*ztaub) + zsigmab / (2*pi*eps_0*fHz)
    zcxu = (zep0u - eps_winf) / (1 + (2*pi*fHz*ztauu)^2)
    zepwux = eps_winf + zcxu
    zepwuy = zcxu * (2*pi*fHz*ztauu) + zsigmau / (2*pi*eps_0*fHz)

    #Computation of refractive index of water (bound & unbound)
    znb = sqrt( sqrt( zepwbx^2 + zepwby^2) + zepwbx ) / sqrt(2)
    zkb = sqrt( sqrt( zepwbx^2 + zepwby^2) - zepwbx ) / sqrt(2)
    znu = sqrt( sqrt( zepwux^2 + zepwuy^2) + zepwux ) / sqrt(2)
    zku = sqrt( sqrt( zepwux^2 + zepwuy^2) - zepwux ) / sqrt(2)

    #Computation of soil refractive index (nm & km): xmv can be a vector
    zxmvt2 = min(mv,zxmvt)
    zflag = float(mv >= zxmvt)
    znm = znd + (znb - 1) * zxmvt2 + (znu - 1) * (mv-zxmvt) * zflag
    zkm = zkd + zkb * zxmvt2 + zku * (mv-zxmvt) * zflag

    #Computation of soil dielectric constant:
    zepmx = znm^2 - zkm^2
    zepmy = znm * zkm * 2
    #eps = zepmx + i*zepmy
    return zepmx, zepmy

end