text:
"""
$(TYPEDSIGNATURES)

Update flow variables over control volumes
"""
function update!(
    KS::AbstractSolverSet,
    ctr::AV{TC},
    face::AV{TF},
    dt,
    residual::FN; # 1D / 2D
    coll = symbolize(KS.set.collision),
    bc = symbolize(KS.set.boundary),
) where {TC<:Union{ControlVolume,ControlVolume1D},TF<:Union{Interface,Interface1D}}

    sumRes = zero(ctr[1].w)
    sumAvg = zero(ctr[1].w)

    @inbounds @threads for i = 1:KS.pSpace.nx
        ctr[i].w, sumRes, sumAvg = step!(
            face[i].fw,
            ctr[i].w,
            ctr[i].prim,
            face[i+1].fw,
            KS.gas.a,
            KS.ps.dx[i],
            sumRes,
            sumAvg,
        )
    end

    residual = sqrt(sumRes * KS.pSpace.nx) / (sumAvg + 1.e-7)

    if bc[1] == :period
        ctr[0].w = ctr[KS.pSpace.nx].w
        ctr[0].sw = ctr[KS.pSpace.nx].sw
        ctr[KS.pSpace.nx+1].w = ctr[1].w
        ctr[KS.pSpace.nx+1].sw = ctr[1].sw
    elseif bc[1] == :extra
        ctr[0].w = ctr[1].w
    end

    if bc[2] == :extra
        ctr[KS.pSpace.nx+1].w = ctr[KS.pSpace.nx].w
    end

    return residual

end

function update!(
    KS::AbstractSolverSet,
    ctr::AV{TC},
    face::AV{TF},
    dt,
    residual::AA{<:FN}; # 1D / 2D
    coll = symbolize(KS.set.collision),
    bc = symbolize(KS.set.boundary),
) where {TC<:Union{ControlVolume,ControlVolume1D},TF<:Union{Interface,Interface1D}}

    sumRes = zero(ctr[1].w)
    sumAvg = zero(ctr[1].w)

    if ndims(sumRes) == 1
        @inbounds @threads for i = 2:KS.pSpace.nx-1
            step!(
                face[i].fw,
                ctr[i].w,
                ctr[i].prim,
                face[i+1].fw,
                KS.gas.γ,
                KS.ps.dx[i],
                sumRes,
                sumAvg,
            )
        end
    elseif ndims(sumRes) == 2
        @inbounds @threads for i = 2:KS.pSpace.nx-1
            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,
                sumRes,
                sumAvg,
            )
        end
    end

    for i in eachindex(residual)
        residual[i] = sqrt(sumRes[i] * KS.pSpace.nx) / (sumAvg[i] + 1.e-7)
    end

    if bc isa Symbol
        update_boundary!(KS, ctr, face, dt, residual; bc = [bc, bc])
    else
        update_boundary!(KS, ctr, face, dt, residual; bc = bc)
    end

    return nothing

end

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

    sumRes = zero(ctr[1].w)
    sumAvg = zero(ctr[1].w)

    phase = KS.set.space[3:end]
    if phase == "1f1v"
        @inbounds @threads for i = 2:KS.pSpace.nx-1
            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,
                sumRes,
                sumAvg,
                coll,
            )
        end
    elseif phase == "1f3v"
        if KS.set.collision == "fsm"
            @inbounds @threads for i = 2:KS.pSpace.nx-1
                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.gas.γ,
                    KS.gas.fsm.Kn,
                    KS.gas.fsm.nm,
                    KS.gas.fsm.ϕ,
                    KS.gas.fsm.ψ,
                    KS.gas.fsm.χ,
                    KS.ps.dx[i],
                    dt,
                    sumRes,
                    sumAvg,
                    coll,
                )
            end
        else
            @inbounds @threads for i = 2:KS.pSpace.nx-1
                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,
                    sumRes,
                    sumAvg,
                    coll,
                )
            end
        end
    end

    for i in eachindex(residual)
        residual[i] = sqrt(sumRes[i] * KS.pSpace.nx) / (sumAvg[i] + 1.e-7)
    end

    if bc isa Symbol
        update_boundary!(KS, ctr, face, dt, residual; coll = coll, bc = [bc, bc])
    else
        update_boundary!(KS, ctr, face, dt, residual; coll = coll, bc = bc)
    end

    return nothing

end

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

    sumRes = zero(ctr[1].w)
    sumAvg = zero(ctr[1].w)

    if ndims(sumRes) == 1
        @inbounds @threads for i = 2:KS.pSpace.nx-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,
                sumRes,
                sumAvg,
                coll,
            )
        end
    elseif ndims(sumRes) == 2
        @inbounds @threads for i = 2:KS.pSpace.nx-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.mi,
                KS.gas.ni,
                KS.gas.me,
                KS.gas.ne,
                KS.gas.Kn[1],
                KS.gas.Pr,
                KS.ps.dx[i],
                dt,
                sumRes,
                sumAvg,
                coll,
            )
        end
    end

    for i in eachindex(residual)
        residual[i] = sqrt(sumRes[i] * KS.pSpace.nx) / (sumAvg[i] + 1.e-7)
    end

    if bc isa Symbol
        update_boundary!(KS, ctr, face, dt, residual; coll = coll, bc = [bc, bc])
    else
        update_boundary!(KS, ctr, face, dt, residual; coll = coll, bc = bc)
    end

    return nothing

end

function update!(
    KS::AbstractSolverSet,
    ctr::AV{TC},
    face::AV{TF},
    dt,
    residual; # 1D / 2D
    coll = symbolize(KS.set.collision),
    bc = symbolize(KS.set.boundary),
    isMHD = true::Bool,
) where {TC<:ControlVolume1D3F,TF<:Interface1D3F}

    sumRes = zero(ctr[1].w)
    sumAvg = zero(ctr[1].w)

    @inbounds @threads for i = 2:KS.pSpace.nx-1
        step!(KS, face[i], ctr[i], face[i+1], KS.ps.dx[i], dt, sumRes, sumAvg, coll, isMHD)
    end

    for i in eachindex(residual)
        residual[i] = sqrt(sumRes[i] * KS.pSpace.nx) / (sumAvg[i] + 1.e-7)
    end

    if bc isa Symbol
        update_boundary!(
            KS,
            ctr,
            face,
            dt,
            residual;
            coll = coll,
            bc = [bc, bc],
            isMHD = isMHD,
        )
    else
        update_boundary!(KS, ctr, face, dt, residual; coll = coll, bc = bc, isMHD = isMHD)
    end

    #=
    ng = 1 - first(eachindex(KS.pSpace.x))
    if bc == :extra
        for i in 1:ng
            ctr[1-i].w .= ctr[1].w
            ctr[1-i].prim .= ctr[1].prim
            ctr[1-i].h0 .= ctr[1].h0
            ctr[1-i].h1 .= ctr[1].h1
            ctr[1-i].h2 .= ctr[1].h2
            ctr[1-i].E .= ctr[1].E
            ctr[1-i].B .= ctr[1].B
            ctr[1-i].ϕ = deepcopy(ctr[1].ϕ)
            ctr[1-i].ψ = deepcopy(ctr[1].ψ)
            ctr[1-i].lorenz .= ctr[1].lorenz

            ctr[KS.pSpace.nx+i].w .= ctr[KS.pSpace.nx].w
            ctr[KS.pSpace.nx+i].prim .= ctr[KS.pSpace.nx].prim
            ctr[KS.pSpace.nx+i].h0 .= ctr[KS.pSpace.nx].h0
            ctr[KS.pSpace.nx+i].h1 .= ctr[KS.pSpace.nx].h1
            ctr[KS.pSpace.nx+i].h2 .= ctr[KS.pSpace.nx].h2
            ctr[KS.pSpace.nx+i].E .= ctr[KS.pSpace.nx].E
            ctr[KS.pSpace.nx+i].B .= ctr[KS.pSpace.nx].B
            ctr[KS.pSpace.nx+i].ϕ = deepcopy(ctr[KS.pSpace.nx].ϕ)
            ctr[KS.pSpace.nx+i].ψ = deepcopy(ctr[KS.pSpace.nx].ψ)
            ctr[KS.pSpace.nx+i].lorenz .= ctr[KS.pSpace.nx].lorenz
        end
    elseif bc == :period
        for i in 1:ng
            ctr[1-i].w .= ctr[KS.pSpace.nx+1-i].w
            ctr[1-i].prim .= ctr[KS.pSpace.nx+1-i].prim
            ctr[1-i].h0 .= ctr[KS.pSpace.nx+1-i].h0
            ctr[1-i].h1 .= ctr[KS.pSpace.nx+1-i].h1
            ctr[1-i].h2 .= ctr[KS.pSpace.nx+1-i].h2
            ctr[1-i].E .= ctr[KS.pSpace.nx+1-i].E
            ctr[1-i].B .= ctr[KS.pSpace.nx+1-i].B
            ctr[1-i].ϕ = deepcopy(ctr[KS.pSpace.nx+1-i].ϕ)
            ctr[1-i].ψ = deepcopy(ctr[KS.pSpace.nx+1-i].ψ)
            ctr[1-i].lorenz .= ctr[KS.pSpace.nx+1-i].lorenz

            ctr[KS.pSpace.nx+i].w .= ctr[i].w
            ctr[KS.pSpace.nx+i].prim .= ctr[i].prim
            ctr[KS.pSpace.nx+i].h0 .= ctr[i].h0
            ctr[KS.pSpace.nx+i].h1 .= ctr[i].h1
            ctr[KS.pSpace.nx+i].h2 .= ctr[i].h2
            ctr[KS.pSpace.nx+i].E .= ctr[i].E
            ctr[KS.pSpace.nx+i].B .= ctr[i].B
            ctr[KS.pSpace.nx+i].ϕ = deepcopy(ctr[i].ϕ)
            ctr[KS.pSpace.nx+i].ψ = deepcopy(ctr[i].ψ)
            ctr[KS.pSpace.nx+i].lorenz .= ctr[i].lorenz
        end
    elseif bc == :balance
        @. ctr[0].w = 0.5 * (ctr[-1].w + ctr[1].w)
        @. ctr[0].prim = 0.5 * (ctr[-1].prim + ctr[1].prim)
        @. ctr[0].h0 = 0.5 * (ctr[-1].h0 + ctr[1].h0)
        @. ctr[0].h1 = 0.5 * (ctr[-1].h1 + ctr[1].h1)
        @. ctr[0].h2 = 0.5 * (ctr[-1].h2 + ctr[1].h2)
        @. ctr[0].E = 0.5 * (ctr[-1].E + ctr[1].E)
        @. ctr[0].B = 0.5 * (ctr[-1].B + ctr[1].B)
        ctr[0].ϕ = 0.5 * (ctr[-1].ϕ + ctr[1].ϕ)
        ctr[0].ψ = 0.5 * (ctr[-1].ψ + ctr[1].ψ)
        @. ctr[0].lorenz = 0.5 * (ctr[-1].lorenz + ctr[1].lorenz)

        @. ctr[KS.pSpace.nx+1].w = 0.5 * (ctr[KS.pSpace.nx].w + ctr[KS.pSpace.nx+2].w)
        @. ctr[KS.pSpace.nx+1].prim = 0.5 * (ctr[KS.pSpace.nx].prim + ctr[KS.pSpace.nx+2].prim)
        @. ctr[KS.pSpace.nx+1].h0 = 0.5 * (ctr[KS.pSpace.nx].h0 + ctr[KS.pSpace.nx+2].h0)
        @. ctr[KS.pSpace.nx+1].h1 = 0.5 * (ctr[KS.pSpace.nx].h1 + ctr[KS.pSpace.nx+2].h1)
        @. ctr[KS.pSpace.nx+1].h2 = 0.5 * (ctr[KS.pSpace.nx].h2 + ctr[KS.pSpace.nx+2].h2)
        @. ctr[KS.pSpace.nx+1].E = 0.5 * (ctr[KS.pSpace.nx].E + ctr[KS.pSpace.nx+2].E)
        @. ctr[KS.pSpace.nx+1].B = 0.5 * (ctr[KS.pSpace.nx].B + ctr[KS.pSpace.nx+2].B)
        ctr[KS.pSpace.nx+1].ϕ = 0.5 * (ctr[KS.pSpace.nx].ϕ + ctr[KS.pSpace.nx+2].ϕ)
        ctr[KS.pSpace.nx+1].ψ = 0.5 * (ctr[KS.pSpace.nx].ψ + ctr[KS.pSpace.nx+2].ψ)
        @. ctr[KS.pSpace.nx+1].lorenz = 0.5 * (ctr[KS.pSpace.nx].lorenz + ctr[KS.pSpace.nx+2].lorenz)
    else
    end
    =#
    return nothing

end

function update!(
    KS::AbstractSolverSet,
    ctr::AV{TC},
    face::AV{TF},
    dt,
    residual; # 1D / 2D
    coll = symbolize(KS.set.collision),
    bc = symbolize(KS.set.boundary),
    isMHD = true::Bool,
) where {TC<:ControlVolume1D4F,TF<:Interface1D4F}

    sumRes = zero(ctr[1].w)
    sumAvg = zero(ctr[1].w)

    @inbounds @threads for i = 2:KS.pSpace.nx-1
        step!(KS, face[i], ctr[i], face[i+1], KS.ps.dx[i], dt, sumRes, sumAvg, coll, isMHD)
    end

    for i in eachindex(residual)
        residual[i] = sqrt(sumRes[i] * KS.pSpace.nx) / (sumAvg[i] + 1.e-7)
    end

    if bc isa Symbol
        update_boundary!(
            KS,
            ctr,
            face,
            dt,
            residual;
            coll = coll,
            bc = [bc, bc],
            isMHD = isMHD,
        )
    else
        update_boundary!(KS, ctr, face, dt, residual; coll = coll, bc = bc, isMHD = isMHD)
    end

    return nothing

end

function update!(
    KS::AbstractSolverSet,
    ctr::AM{TC},
    a1face::AM{TF},
    a2face::AM{TF},
    dt,
    residual;
    coll = symbolize(KS.set.collision),
    bc = symbolize(KS.set.boundary),
) 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

    sumRes = zero(ctr[1].w)
    sumAvg = zero(ctr[1].w)

    if ndims(sumRes) == 1

        @inbounds @threads for j = 2:ny-1
            for i = 2:nx-1
                step!(
                    ctr[i, j].w,
                    ctr[i, j].prim,
                    a1face[i, j].fw,
                    a1face[i+1, j].fw,
                    a2face[i, j].fw,
                    a2face[i, j+1].fw,
                    KS.gas.γ,
                    dx[i, j] * dy[i, j],
                    sumRes,
                    sumAvg,
                    coll,
                )
            end
        end

    end

    for i in eachindex(residual)
        residual[i] = sqrt(sumRes[i] * nx * ny) / (sumAvg[i] + 1.e-7)
    end

    if bc isa Symbol
        update_boundary!(
            KS,
            ctr,
            a1face,
            a2face,
            dt,
            residual;
            coll = coll,
            bc = [bc, bc, bc, bc],
        )
    else
        update_boundary!(KS, ctr, a1face, a2face, dt, residual; coll = coll, bc = bc)
    end

    return nothing

end

function update!(
    KS::AbstractSolverSet,
    ctr::AM{TC},
    a1face::AM{TF},
    a2face::AM{TF},
    dt,
    residual;
    coll = symbolize(KS.set.collision),
    bc = symbolize(KS.set.boundary),
) 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

    sumRes = zero(ctr[1].w)
    sumAvg = zero(ctr[1].w)

    if ndims(sumRes) == 1

        @inbounds @threads for j = 2:ny-1
            for i = 2:nx-1
                step!(
                    ctr[i, j].w,
                    ctr[i, j].prim,
                    ctr[i, j].f,
                    a1face[i, j].fw,
                    a1face[i, j].ff,
                    a1face[i+1, j].fw,
                    a1face[i+1, j].ff,
                    a2face[i, j].fw,
                    a2face[i, j].ff,
                    a2face[i, j+1].fw,
                    a2face[i, j+1].ff,
                    KS.vSpace.u,
                    KS.vSpace.v,
                    KS.vSpace.weights,
                    KS.gas.γ,
                    KS.gas.μᵣ,
                    KS.gas.ω,
                    KS.gas.Pr,
                    dx[i, j] * dy[i, j],
                    dt,
                    sumRes,
                    sumAvg,
                    coll,
                )
            end
        end

    end

    for i in eachindex(residual)
        residual[i] = sqrt(sumRes[i] * nx * ny) / (sumAvg[i] + 1.e-7)
    end

    if bc isa Symbol
        update_boundary!(
            KS,
            ctr,
            a1face,
            a2face,
            dt,
            residual;
            coll = coll,
            bc = [bc, bc, bc, bc],
        )
    else
        update_boundary!(KS, ctr, a1face, a2face, dt, residual; coll = coll, bc = bc)
    end

    return nothing

end

#--- 2d2f ---#
function update!(
    KS::AbstractSolverSet,
    ctr::AM{TC},
    a1face::AM{TF},
    a2face::AM{TF},
    dt,
    residual;
    coll = symbolize(KS.set.collision),
    bc = symbolize(KS.set.boundary),
) 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

    sumRes = zero(ctr[1].w)
    sumAvg = zero(ctr[1].w)

    if ndims(sumRes) == 1

        @inbounds @threads for j = 2:ny-1
            for i = 2:nx-1
                step!(
                    ctr[i, j].w,
                    ctr[i, j].prim,
                    ctr[i, j].h,
                    ctr[i, j].b,
                    a1face[i, j].fw,
                    a1face[i, j].fh,
                    a1face[i, j].fb,
                    a1face[i+1, j].fw,
                    a1face[i+1, j].fh,
                    a1face[i+1, j].fb,
                    a2face[i, j].fw,
                    a2face[i, j].fh,
                    a2face[i, j].fb,
                    a2face[i, j+1].fw,
                    a2face[i, j+1].fh,
                    a2face[i, 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[i, j] * dy[i, j],
                    dt,
                    sumRes,
                    sumAvg,
                    coll,
                )
            end
        end

    end

    for i in eachindex(residual)
        residual[i] = sqrt(sumRes[i] * nx * ny) / (sumAvg[i] + 1.e-7)
    end

    if bc isa Symbol
        update_boundary!(
            KS,
            ctr,
            a1face,
            a2face,
            dt,
            residual;
            coll = coll,
            bc = [bc, bc, bc, bc],
        )
    else
        update_boundary!(KS, ctr, a1face, a2face, dt, residual; coll = coll, bc = bc)
    end

    return nothing

end

function update!(
    KS::AbstractSolverSet,
    ctr::AV{ControlVolumeUS},
    face::AV{T},
    dt,
    residual; # 1D / 2D
    coll = symbolize(KS.set.collision),
    bc = symbolize(KS.set.boundary),
) where {T<:Union{Interface,Interface2D}}

    sumRes = zero(ctr[1].w)
    sumAvg = zero(ctr[1].w)

    @inbounds @threads for i in eachindex(ctr)
        if KS.ps.cellType[i] in (0, 2)
            dirc = [sign(dot(ctr[i].n[j], face[KS.ps.cellFaces[i, j]].n)) for j = 1:3]

            step!(
                ctr[i].w,
                ctr[i].prim,
                face[KS.ps.cellFaces[i, 1]].fw,
                face[KS.ps.cellFaces[i, 2]].fw,
                face[KS.ps.cellFaces[i, 3]].fw,
                KS.gas.γ,
                KS.pSpace.cellArea[i],
                dirc,
                sumRes,
                sumAvg,
            )
        end
    end

    for i in eachindex(residual)
        residual[i] = sqrt(sumRes[i] * size(KS.ps.cellid, 1)) / (sumAvg[i] + 1.e-7)
    end

    update_boundary!(KS, ctr, face, dt, residual; coll = coll, bc = bc)

    return nothing

end

function update!(
    KS::AbstractSolverSet,
    ctr::AV{ControlVolumeUS1F},
    face::AV{T},
    dt,
    residual; # 1D / 2D
    coll = symbolize(KS.set.collision),
    bc = symbolize(KS.set.boundary),
) where {T<:Union{Interface1F,Interface2D1F}}

    sumRes = zero(ctr[1].w)
    sumAvg = zero(ctr[1].w)

    @inbounds @threads for i in eachindex(ctr)
        if KS.ps.cellType[i] in (0, 2)
            dirc = [sign(dot(ctr[i].n[j], face[KS.ps.cellFaces[i, j]].n)) for j = 1:3]

            step!(
                ctr[i].w,
                ctr[i].prim,
                ctr[i].f,
                face[KS.ps.cellFaces[i, 1]].fw,
                face[KS.ps.cellFaces[i, 1]].ff,
                face[KS.ps.cellFaces[i, 2]].fw,
                face[KS.ps.cellFaces[i, 2]].ff,
                face[KS.ps.cellFaces[i, 3]].fw,
                face[KS.ps.cellFaces[i, 3]].ff,
                KS.vSpace.u,
                KS.vSpace.v,
                KS.vSpace.weights,
                KS.gas.K,
                KS.gas.γ,
                KS.gas.μᵣ,
                KS.gas.ω,
                KS.gas.Pr,
                KS.pSpace.cellArea[i],
                dirc,
                dt,
                sumRes,
                sumAvg,
                coll,
            )
        end
    end

    for i in eachindex(residual)
        residual[i] = sqrt(sumRes[i] * size(KS.ps.cellid, 1)) / (sumAvg[i] + 1.e-7)
    end

    update_boundary!(KS, ctr, face, dt, residual; coll = coll, bc = bc)

    return nothing

end

function update!(
    KS::AbstractSolverSet,
    ctr::AV{ControlVolumeUS2F},
    face::AV{T},
    dt,
    residual; # 1D / 2D
    coll = symbolize(KS.set.collision),
    bc = symbolize(KS.set.boundary),
) where {T<:Union{Interface2F,Interface2D2F}}

    sumRes = zero(ctr[1].w)
    sumAvg = zero(ctr[1].w)

    @inbounds @threads for i in eachindex(ctr)
        if KS.ps.cellType[i] in (0, 2)
            dirc = [sign(dot(ctr[i].n[j], face[KS.ps.cellFaces[i, j]].n)) for j = 1:3]

            step!(
                ctr[i].w,
                ctr[i].prim,
                ctr[i].h,
                ctr[i].b,
                face[KS.ps.cellFaces[i, 1]].fw,
                face[KS.ps.cellFaces[i, 1]].fh,
                face[KS.ps.cellFaces[i, 1]].fb,
                face[KS.ps.cellFaces[i, 2]].fw,
                face[KS.ps.cellFaces[i, 2]].fh,
                face[KS.ps.cellFaces[i, 2]].fb,
                face[KS.ps.cellFaces[i, 3]].fw,
                face[KS.ps.cellFaces[i, 3]].fh,
                face[KS.ps.cellFaces[i, 3]].fb,
                KS.vSpace.u,
                KS.vSpace.v,
                KS.vSpace.weights,
                KS.gas.K,
                KS.gas.γ,
                KS.gas.μᵣ,
                KS.gas.ω,
                KS.gas.Pr,
                KS.pSpace.cellArea[i],
                dirc,
                dt,
                sumRes,
                sumAvg,
                coll,
            )
        end
    end

    for i in eachindex(residual)
        residual[i] = sqrt(sumRes[i] * size(KS.ps.cellid, 1)) / (sumAvg[i] + 1.e-7)
    end

    update_boundary!(KS, ctr, face, dt, residual; coll = coll, bc = bc)

    return nothing

end


text:
ENV["MOCHA_USE_CUDA"] = "true"

using Mocha
using JLD
srand(12345678)

dataLayer = AsyncHDF5DataLayer(name="train-data", source="train.txt", batch_size=64, shuffle=false)

shallowLayers = [
	ConvolutionLayer(name="conv1", n_filter=20, kernel=(5,5), bottoms=[:data], tops=[:conv])
	PoolingLayer(name="pool1", kernel=(2,2), stride=(2,2), bottoms=[:conv], tops=[:pool])
	InnerProductLayer(name="ip1", output_dim=18, bottoms=[:pool], tops=[:ip1])
	SoftmaxLossLayer(name="loss", bottoms=[:ip1,:label])
]

backend = DefaultBackend()
init(backend)

net = Net("MNIST-train", backend, [dataLayer, shallowLayers...])
println(net)
open("net.dot", "w") do out net2dot(out, net) end

f = jldopen("shallowNet.jld", "r")
load_network(f, net)

array = Array{Float32, 2}(size(net.output_blobs[:ip1])...)
@show net.output_blobs[:ip1]
copy!(net.output_blobs[:ip1], array)
@show array

exp_dir = "snapshots-$(Mocha.default_backend_type)"

method = SGD()
params = make_solver_parameters(method, max_iter=100, regu_coef=0.0005,
                                mom_policy=MomPolicy.Fixed(0.9),
                                lr_policy=LRPolicy.Inv(0.01, 0.0001, 0.75),
                                load_from=exp_dir)
solver = Solver(method, params)

setup_coffee_lounge(solver, save_into="$exp_dir/statistics.jld", every_n_iter=1000)

# report training progress every 100 iterations
add_coffee_break(solver, TrainingSummary(), every_n_iter=100)

# save snapshots every 5000 iterations
#add_coffee_break(solver, Snapshot(exp_dir), every_n_iter=5000)

# show performance on test data every 1000 iterations
data_layer_test = HDF5DataLayer(name="test-data", source="test.txt", batch_size=100)
acc_layer = SquareLossLayer(name="test-accuracy", bottoms=[:ip1, :label])
test_net = Net("MNIST-test", backend, [data_layer_test, shallowLayers...])
add_coffee_break(solver, ValidationPerformance(test_net), every_n_iter=1000)

solve(solver, net)

@show net.output_blobs[:ip1]

f = jldopen("shallowNet.jld", "w")
save_network(f, net)
close(f)

#Profile.init(int(1e8), 0.001)
#@profile solve(solver, net)
#open("profile.txt", "w") do out
#  Profile.print(out)
#end

destroy(net)
destroy(test_net)
shutdown(backend)

ext: jl


text:
module NnetTest

using Learn
# using StatPlots; glvisualize(leg=false, show=true)

function doit()

    nin = 10
    nout = 1
    τ = 1000
    w1 = randn(nout, nin)
    w2 = randn(nout, nin)
    # f(x) = tanh(w1*x + w2*sin(20x) + 0.05randn(nout))[1] > 0.0 ? 1.0 : 0.0
    f(x) = tanh(10w1*x)[1] > 0.0 ? 1.0 : 0.0
    inputs = randn(nin,τ)
    targets = zeros(nout,τ)
    for (i,x) in enumerate(DataSubset(inputs, 1:nobs(inputs)))
        targets[i] = f(x)
    end
    # scatter(inputs', targets', layout=nin, m=(3,0.2,stroke(0)), smooth=true)

    # build our model
    nh = 10
    chain = Chain(
        Affine(nin,nh),
        Activation(:relu,nh),
        Affine(nh,nout),
        Activation(:logistic,nout)
    )
    obj = objective(chain, CrossentropyLoss(), L1Penalty(1e-6))

    learner = make_learner(
        GradientLearner(FixedLR(1e-4), SGD()),
        maxiter = 10000,
        oniter = (model,i) -> begin
            @show i
        end
    )
    learn!(obj, learner, MiniBatches((inputs, targets), 5))

    # eval
    for (i,(input,target)) in enumerate(eachobs((inputs,targets)))
        @show i,transform!(chain,input),target
    end

    # store for inspection
    obj,learner
end
obj,learner = doit()
end #module

max_stars_repo_name: JuliaML/StochasticOptimization

max_issues_repo_name: JuliaML/StochasticOptimization

max_forks_repo_name: JuliaML/StochasticOptimization


text:
t_start = time_ns()

import Mercury as Hg
using ZMQ
using Sockets

using StaticArrays
using LinearAlgebra

if !isdefined(@__MODULE__, :TestMsg)
    include(joinpath(pkgdir(Hg), "test", "jlout", "test_msg_pb.jl"))
end

mutable struct SimpleNode <: Hg.Node
    ctx::ZMQ.Context   # NOTE: can we get rid of needing to store the ZMQ context?
    nodeio::Hg.NodeIO
    test_msg::TestMsg
end

function SimpleNode(ctx; rate = 10)
    test_msg = TestMsg(x = 0, y = 0, z = 0)
    nodeio = Hg.NodeIO(rate = rate)
    SimpleNode(ctx, nodeio, test_msg)
end

function Hg.setupIO!(node::SimpleNode, nodeio::Hg.NodeIO)
    # Create a publisher
    ctx = node.ctx
    addr = ip"127.0.0.1"
    port = 5555
    pub = Hg.ZmqPublisher(ctx, addr, port, name = "test_pub")

    # Register the publisher to publish the `TestMsg` stored in the node
    Hg.add_publisher!(nodeio, node.test_msg, pub)
end

function Hg.compute(node::SimpleNode)
    # Update the internal message
    A = @SMatrix rand(3, 3)
    A = A'A
    Achol = cholesky(A)
    b = SA[node.test_msg.x, node.test_msg.y, node.test_msg.z]
    x = Achol \ b
    node.test_msg.x += 1
    node.test_msg.y = x[2]
    node.test_msg.z = x[3]

    # Publish the message
    #   Mercury will automatically encode the message as a string of bytes and send it over ZMQ
    Hg.publish(Hg.getIO(node).pubs[1])  # NOTE: can we make this easier, maybe by using a Dict?

    if (node.test_msg.x % node.nodeio.opts.rate) == 0
        println("Sent x value of ", node.test_msg.x)
    end
end

function launch_simple_node(; rate = 10, time = 5)
    node = SimpleNode(ZMQ.context(), rate = rate)
    Hg.setupIO!(node, Hg.getIO(node))
    @time Hg.compute(node)
    t_since_start = (time_ns() - t_start) / 1e9
    println("Took $t_since_start seconds to precompile.")
    task = @async Hg.launch(node)
    sleep(5)
    Hg.stopnode(node)
    node
end

launch_simple_node(rate = 100)


text:
"""
The R NAMED property, represented by 2 bits in the info field. This can take
values 0,1 or 2, corresponding to whether it is bound to 0,1 or 2 or more
symbols. See
http://cran.r-project.org/doc/manuals/r-patched/R-exts.html#Named-objects-and-copying
"""
function bound(s::Ptr{S}) where S<:Sxp
    u = unsafe_load(convert(Ptr{UnknownSxp},s))
    (u.info >>> 6) & 0x03
end

"""
Sxp methods for `length` return the R length.

`Rf_xlength` handles Sxps that are not vector-like and R's
"long vectors", which have a negative value for the `length` member.
"""
length(s::Ptr{S}) where S<:Sxp = Int(ccall((:Rf_xlength,libR),Cptrdiff_t,(Ptr{S},),s))
length(r::RObject) = length(r.p)

## Predicates applied to an Sxp
##
## Many of these are unneeded but a few extra definitions is not a big deal
for sym in (:isArray,:isComplex,:isEnvironment,:isExpression,:isFactor,
            :isFrame,:isFree,:isFunction,:isInteger,:isLanguage,:isList,
            :isLogical,:isSymbol,:isMatrix,:isNewList,:isNull,:isNumeric,
            :isNumber,:isObject,:isOrdered,:isPairList,:isPrimitive,
            :isReal,:isS4,:isString,:isTs,:isUnordered,:isUnsorted,
            :isUserBinop,:isValidString,:isValidStringF,:isVector,
            :isVectorAtomic,:isVectorizable,:isVectorList)
    @eval begin
        $sym(s::Ptr{S}) where S<:Sxp = ccall(($(string("Rf_",sym)),libR),Bool,(Ptr{SxpPtrInfo},),s)
        $sym(r::RObject) = $sym(r.p)
    end
end

for (S, J) in ((:LglSxp, "LOGICAL"), (:IntSxp, "INTEGER"),
                (:RealSxp, "REAL"), (:CplxSxp, "COMPLEX"), (:RawSxp, "RAW"))
    @eval dataptr(s::Ptr{$S}) = convert(Ptr{eltype($S)}, ccall(($J,libR), Ptr{Cvoid}, (Ptr{UnknownSxp},), s))
end

"""
Represent the contents of a VectorSxp type as a `Vector`.

This does __not__ copy the contents.  If the argument is not named (in R) or
otherwise protected from R's garbage collection (e.g. by keeping the
containing RObject in scope) the contents of this vector can be modified or
could cause a memory error when accessed.

The contents are as stored in R.  Missing values (NA's) are represented
in R by sentinels.  Missing data values in RealSxp and CplxSxp show
up as `NaN` and `NaN + NaNim`, respectively.  Missing data in IntSxp show up
as `-2147483648`, the minimum 32-bit integer value.  Internally a `LglSxp` is
represented as `Vector{Int32}`.  The convention is that `0` is `false`,
`-2147483648` is `NA` and all other values represent `true`.
"""
unsafe_vec(s::Ptr{S}) where S<:VectorSxp = unsafe_wrap(Array, dataptr(s), length(s))
unsafe_vec(r::RObject{S}) where S<:VectorSxp = unsafe_vec(r.p)

"""
The same as `unsafe_vec`, except returns an appropriately sized array.
"""
unsafe_array(s::Ptr{S}) where S<:VectorSxp =  unsafe_wrap(Array, dataptr(s), size(s))
unsafe_array(r::RObject{S}) where S<:VectorSxp = unsafe_array(r.p)

# Sxp iterator
IteratorSize(x::Ptr{S}) where S<:Sxp = Base.SizeUnknown()
IteratorEltype(x::Ptr{S}) where S<:Sxp = Base.EltypeUnknown()
pairs(s::Ptr{S}) where S<:Sxp = Pairs(s, Base.OneTo(length(s)))

# RObject iterator
@inline iterate(x::RObject) = iterate(x.p)
@inline iterate(x::RObject, state) = iterate(x.p, state)
IteratorSize(x::RObject) = IteratorSize(x.p)
IteratorEltype(x::RObject) = IteratorEltype(x.p)
pairs(r::RObject{S}) where S<:Sxp = Pairs(r, Base.OneTo(length(r)))

# NilSxp

IteratorSize(x::Ptr{NilSxp}) = Base.HasLength()
iterate(s::Ptr{NilSxp}) = nothing
iterate(s::Ptr{NilSxp}, state) = nothing

# VectorSxp

getindex(r::RObject{S}, I...) where S<:VectorSxp = getindex(sexp(r), I...)
getindex(r::RObject{S}, I::AbstractArray) where S<:VectorSxp = getindex(sexp(r), I)
setindex!(r::RObject{S}, value, keys...) where S<:VectorSxp = setindex!(sexp(r), value, keys...)
setindex!(r::RObject{S}, ::Missing, keys...) where S<:VectorSxp = setindex!(sexp(r), naeltype(S), keys...)

IteratorSize(x::Ptr{S}) where S<:VectorSxp = Base.HasLength()
IteratorEltype(x::Ptr{S}) where S<:VectorSxp = Base.HasEltype()
iterate(s::Ptr{S}) where S<:VectorSxp = iterate(s, 0)
function iterate(s::Ptr{S}, state) where S<:VectorSxp
    state ≥ length(s) && return nothing
    state += 1
    (s[state], state)
end

"""
Set element of a VectorSxp by a label.
"""
function getindex(s::Ptr{S}, label::AbstractString) where S<:VectorSxp
    ls = getnames(s)
    for (i,l) in enumerate(ls)
        if rcopy(String, l) == label
            return s[i]
        end
    end
    throw(BoundsError())
end
getindex(s::Ptr{S}, label::Symbol) where S<:VectorSxp = getindex(s,string(label))

"""
Set element of a VectorSxp by a label.
"""
function setindex!(s::Ptr{S}, value::Ptr{T}, label::AbstractString) where {S<:VectorSxp, T<:Sxp}
    ls = getnames(s)
    for (i,l) in enumerate(ls)
        if rcopy(String, l) == label
            s[i] = value
            return
        end
    end
    throw(BoundsError())
end
setindex!(s::Ptr{S}, value, label::AbstractString) where S<:VectorSxp = setindex!(s, sexp(value), label)
setindex!(s::Ptr{S}, value, label::Symbol) where S<:VectorSxp = setindex!(s, value, string(label))

# VectorAtomicSxp
getindex(s::Ptr{S}, I::Integer) where S<:VectorAtomicSxp = getindex(unsafe_vec(s),I)
getindex(s::Ptr{S}, I::Integer...) where S<:VectorAtomicSxp = getindex(unsafe_array(s),I...)
getindex(s::Ptr{S}, I::AbstractVector) where S<:VectorAtomicSxp = getindex(unsafe_vec(s),I)

function setindex!(s::Ptr{S}, value, I::Integer...) where S<:VectorAtomicSxp
    setindex!(unsafe_array(s), value, I...)
end
function setindex!(s::Ptr{S}, value, I::Integer) where S<:VectorAtomicSxp
    setindex!(unsafe_vec(s), value, I)
end

# VectorList

IteratorEltype(x::RObject{S}) where S<:VectorListSxp = Base.EltypeUnknown()
iterate(s::RObject{S}) where S<:VectorListSxp = iterate(s, 0)
function iterate(s::RObject{S}, state) where S<:VectorListSxp
    state ≥ length(s) && return nothing
    state += 1
    (RObject(s[state]), state)
end
getindex(r::RObject{S}, I...) where S<:VectorListSxp = RObject(getindex(sexp(r), I...))

# StrSxp

function getindex(s::Ptr{StrSxp}, key::Integer)
    c = ccall((:STRING_ELT, libR), Ptr{CharSxp}, (Ptr{StrSxp}, Cint), s, key-1)
end

function setindex!(s::Ptr{StrSxp}, value::Ptr{CharSxp}, key::Integer)
    1 <= key <= length(s) || throw(BoundsError())
    ccall((:SET_STRING_ELT,libR), Nothing,
          (Ptr{StrSxp},Cptrdiff_t, Ptr{CharSxp}),
          s, key-1, value)
    value
end
function setindex!(s::Ptr{StrSxp}, value::AbstractString, key::Integer)
    setindex!(s,sexp(CharSxp,value),key)
end

# VecSxp and ExprSxp

function getindex(s::Ptr{S}, key::Integer) where S<:Union{VecSxp,ExprSxp}
    sexp(ccall((:VECTOR_ELT, libR), Ptr{UnknownSxp}, (Ptr{S}, Cint), s, key-1))
end

function setindex!(s::Ptr{S}, value::Ptr{T}, key::Integer) where {S<:Union{VecSxp,ExprSxp}, T<:Sxp}
    1 <= key <= length(s) || throw(BoundsError())
    ccall((:SET_VECTOR_ELT,libR), Ptr{T},
          (Ptr{S},Cptrdiff_t, Ptr{T}),
          s, key-1, value)
end
function setindex!(s::Ptr{S}, value, key::Integer) where {S<:Union{VecSxp,ExprSxp}}
    setindex!(s,sexp(value),key)
end

# PairListSxps

cdr(s::Ptr{S}) where S<:PairListSxp = sexp(ccall((:CDR,libR),Ptr{UnknownSxp},(Ptr{S},),s))
car(s::Ptr{S}) where S<:PairListSxp = sexp(ccall((:CAR,libR),Ptr{UnknownSxp},(Ptr{S},),s))
tag(s::Ptr{S}) where S<:PairListSxp = sexp(ccall((:TAG,libR),Ptr{UnknownSxp},(Ptr{S},),s))

function setcar!(s::Ptr{S}, c::Ptr{T}) where {S<:PairListSxp, T<:Sxp}
    ccall((:SETCAR,libR),Ptr{Cvoid},(Ptr{S},Ptr{T}),s,c)
    nothing
end
setcar!(s::Ptr{S}, c::RObject{T}) where {S<:PairListSxp, T<:Sxp} = setcar!(s,sexp(c))

function settag!(s::Ptr{S}, c::Ptr{T}) where {S<:PairListSxp, T<:Sxp}
    ccall((:SET_TAG,libR),Nothing,(Ptr{S},Ptr{T}),s,c)
    nothing
end
settag!(s::Ptr{S}, c::RObject{T}) where {S<:PairListSxp, T<:Sxp} = settag!(s,sexp(c))

function setcdr!(s::Ptr{S}, c::Ptr{T}) where {S<:PairListSxp, T<:Sxp}
    ccall((:SETCDR,libR),Ptr{Cvoid},(Ptr{S},Ptr{T}),s,c)
    nothing
end
setcdr!(s::Ptr{S}, c::RObject{T}) where {S<:PairListSxp, T<:Sxp} = setcdr!(s,sexp(c))

iterate(s::Ptr{S}) where S<:PairListSxp = iterate(s, s)
function iterate(s::Ptr{S}, state) where S<:PairListSxp
    state == sexp(Const.NilValue) && return nothing
    car(state), cdr(state)
end

iterate(s::RObject{S}) where S<:PairListSxp = iterate(s, s.p)
function iterate(s::RObject{S}, state) where S<:PairListSxp
    state == sexp(Const.NilValue) && return nothing
    RObject(car(state)), cdr(state)
end

# iterator for PairListSxp
IteratorSize(x::Pairs{K, V, I, Ptr{S}}) where {K, V, I, S<:PairListSxp} = Base.SizeUnknown()
IteratorEltype(x::Pairs{K, V, I, Ptr{S}}) where {K, V, I, S<:PairListSxp} = Base.EltypeUnknown()
@inline iterate(x::Pairs{K, V, I, Ptr{S}}) where {K, V, I, S<:PairListSxp} = iterate(x, x.data)
@inline function iterate(x::Pairs{K, V, I, Ptr{S}}, state) where {K, V, I, S<:PairListSxp}
    state == sexp(Const.NilValue) && return nothing
    (tag(state), car(state)), cdr(state)
end

IteratorSize(x::Pairs{K, V, I, RObject{S}}) where {K, V, I, S<:PairListSxp} = Base.SizeUnknown()
IteratorEltype(x::Pairs{K, V, I, RObject{S}}) where {K, V, I, S<:PairListSxp} = Base.EltypeUnknown()
@inline iterate(x::Pairs{K, V, I, RObject{S}}) where {K, V, I, S<:PairListSxp} = iterate(x, x.data.p)
@inline function iterate(x::Pairs{K, V, I, RObject{S}}, state) where {K, V, I, S<:PairListSxp}
    state == sexp(Const.NilValue) && return nothing
    (RObject(tag(state)), RObject(car(state))), cdr(state)
end

"extract the i-th element of a PairListSxp"
function getindex(l::Ptr{S},I::Integer) where S<:PairListSxp
    1 ≤ I ≤ length(l) || throw(BoundsError())
    for i in 2:I
        l = cdr(l)
    end
    car(l)
end
getindex(r::RObject{S},I::Integer) where S<:PairListSxp = RObject(getindex(sexp(r),I))

"extract an element from a PairListSxp by label"
function getindex(s::Ptr{S}, label::AbstractString) where S<:PairListSxp
    ls = getnames(s)
    for (i,l) in enumerate(ls)
        if rcopy(String, l) == label
            return s[i]
        end
    end
    throw(BoundsError())
end
getindex(s::Ptr{S}, label::Symbol) where S<:PairListSxp = getindex(s,string(label))
getindex(s::RObject{S}, label) where S<:PairListSxp = RObject(getindex(s.p,label))

"assign value v to the i-th element of a PairListSxp"
function setindex!(l::Ptr{S},v::Ptr{T},I::Integer) where {S<:PairListSxp, T<:Sxp}
    1 ≤ I ≤ length(l) || throw(BoundsError())
    for i in 2:I
        l = cdr(l)
    end
    setcar!(l,v)
end
function setindex!(s::Ptr{S}, value, key::Integer) where S<:PairListSxp
    setindex!(s,sexp(value),key)
end

"""
Set element of a PairListSxp by a label.
"""
function setindex!(s::Ptr{S}, value::Ptr{T}, label::AbstractString) where {S<:PairListSxp, T<:Sxp}
    ls = getnames(s)
    for (i,l) in enumerate(ls)
        if rcopy(String, l) == label
            s[i] = value
            return
        end
    end
    throw(BoundsError())
end
setindex!(s::Ptr{S}, value, label::AbstractString) where S<:PairListSxp = setindex!(s, sexp(value), label)
setindex!(s::Ptr{S}, value, label::Symbol) where S<:PairListSxp = setindex!(s, value, string(label))

# for RObjects
setindex!(r::RObject{S}, value, key) where S<:PairListSxp = setindex!(sexp(r), value, key)

# S4Sxp
"extract an element from a S4Sxp by label"
function getindex(s::Ptr{S4Sxp}, sym::Ptr{SymSxp})
    if ccall((:R_has_slot, libR), Int, (Ptr{S4Sxp}, Ptr{SymSxp}), s, sym) == 1
        return sexp(ccall((:R_do_slot, libR), Ptr{UnknownSxp}, (Ptr{S4Sxp}, Ptr{SymSxp}), s, sym))
    else
        throw(BoundsError())
    end
end
getindex(s::Ptr{S4Sxp}, sym) = getindex(s,sexp(SymSxp, sym))
getindex(s::RObject{S4Sxp}, sym) = RObject(getindex(s.p,sym))

"extract an element from a S4Sxp by label"
function setindex!(s::Ptr{S4Sxp}, value::Ptr{T}, sym::Ptr{SymSxp}) where T<:Sxp
    protect(value)
    try
        t = rcall_p(findNamespace("methods")[:checkSlotAssignment], s, rcopy(String, sym), value)
        sexp(ccall((:R_do_slot_assign, libR), Ptr{UnknownSxp}, (Ptr{S4Sxp}, Ptr{SymSxp}, Ptr{T}), s, sym, t))
    finally
        unprotect(1)
    end
end
setindex!(s::Ptr{S4Sxp}, value, sym) = setindex!(s, sexp(value), sexp(SymSxp, sym))
# for RObjects
setindex!(s::RObject{S4Sxp}, value, sym) = setindex!(sexp(s), value, sym)

"Return a particular attribute of an RObject"
function getattrib(s::Ptr{S}, sym::Ptr{SymSxp}) where S<:Sxp
    sexp(ccall((:Rf_getAttrib,libR),Ptr{UnknownSxp},(Ptr{S},Ptr{SymSxp}),s,sym))
end
getattrib(s::Ptr{S}, sym) where S<:Sxp = getattrib(s,sexp(SymSxp,sym))
getattrib(r::RObject, sym) = RObject(getattrib(r.p,sym))

"Set a particular attribute of an RObject"
function setattrib!(s::Ptr{S},sym::Ptr{SymSxp},t::Ptr{T}) where {S<:Sxp, T<:Sxp}
    ccall((:Rf_setAttrib,libR),Ptr{Cvoid},(Ptr{S},Ptr{SymSxp},Ptr{T}),s,sym,t)
    return nothing
end
setattrib!(s::Ptr{S}, sym, t) where S<:Sxp = setattrib!(s, sexp(SymSxp,sym), sexp(t))
setattrib!(r::RObject, sym, t) = setattrib!(r.p, sym, t)

attributes(s::SxpHead) = sexp(convert(Ptr{SxpHead}, s.attrib))
attributes(s::Sxp) = attributes(s.head)
attributes(s::Ptr{S}) where S<:Sxp = attributes(unsafe_load(s))
attributes(s::RObject) = RObject(attributes(s.p))

"""
Returns the size of an R object.
"""
function size(s::Ptr{S}) where S<:Sxp
    if isFrame(s)
        (length(getattrib(s, Const.RowNamesSymbol)), length(s))
    elseif isArray(s)
        tuple(convert(Array{Int},unsafe_vec(getattrib(s,Const.DimSymbol)))...)
    else
        (length(s),)
    end
end
size(r::RObject) = size(r.p)

"""
Returns the names of an R vector.
"""
getnames(s::Ptr{S}) where S<:Sxp = getattrib(s,Const.NamesSymbol)
getnames(r::RObject) = RObject(getnames(sexp(r)))

"""
Returns the names of an R vector, the result is converted to a Julia symbol array.
"""
names(r::RObject) = rcopy(Vector{Symbol}, getnames(sexp(r)))

"""
Set the names of an R vector.
"""
setnames!(s::Ptr{S}, n::Ptr{StrSxp}) where S<:Sxp = setattrib!(s,Const.NamesSymbol,n)
setnames!(r::RObject, n) = RObject(setnames!(sexp(r),sexp(StrSxp,n)))

"""
Returns the class of an R object.
"""
function getclass(s::Ptr{S}, singleString::Bool=false) where S<:Sxp
    ccall((:R_data_class,libR),Ptr{StrSxp},(Ptr{S},Cint),s,singleString)
end
getclass(s::Ptr{CharSxp}, singleString::Bool=false) = Const.NilValue
getclass(r::RObject, singleString::Bool=false) = RObject(getclass(sexp(r), singleString))

"""
Set the class of an R object.
"""
setclass!(s::Ptr{S},c::Ptr{StrSxp}) where S<:Sxp = setattrib!(s,Const.ClassSymbol,c)
setclass!(r::RObject,c) = RObject(setclass!(sexp(r)),sexp(StrSxp,c))

allocList(n::Int) = ccall((:Rf_allocList,libR),Ptr{ListSxp},(Cint,),n)
allocArray(::Type{S}, n::Integer) where S<:Sxp =
    ccall((:Rf_allocVector,libR),Ptr{S},(Cint,Cptrdiff_t),sexpnum(S),n)

allocArray(::Type{S}, n1::Integer, n2::Integer) where S<:Sxp =
    ccall((:Rf_allocMatrix,libR),Ptr{S},(Cint,Cint,Cint),sexpnum(S),n1,n2)

allocArray(::Type{S}, n1::Integer, n2::Integer, n3::Integer) where S<:Sxp =
    ccall((:Rf_alloc3DArray,libR),Ptr{S},(Cint,Cint,Cint,Cint),sexpnum(S),n1,n2,n3)

function allocArray(::Type{S}, dims::Integer...) where S<:Sxp
    sdims = sexp(RClass{:integer},[dims...])
    ccall((:Rf_allocArray,libR),Ptr{S},(Cint,Ptr{IntSxp}),sexpnum(S),sdims)
end

isnull(s::Ptr{S}) where S<:Sxp = isNull(s)
"""
Check if values correspond to R's NULL object.
"""
isnull(r::RObject) = isnull(r.p)

"""
NA element for each R base class
"""
naeltype(::Type{LglSxp}) = Const.NaInt
naeltype(::Type{IntSxp}) = Const.NaInt
naeltype(::Type{RealSxp}) = Const.NaReal
naeltype(::Type{CplxSxp}) = complex(Const.NaReal)
naeltype(::Type{StrSxp}) = sexp(Const.NaString)
# naeltype(::Type{S}) where S<:Sxp = sexp(Const.NilValue)

naeltype(::Type{RClass{:logical}}) = Const.NaInt
naeltype(::Type{RClass{:integer}}) = Const.NaInt
naeltype(::Type{RClass{:numeric}}) = Const.NaReal
naeltype(::Type{RClass{:complex}}) = complex(Const.NaReal,Const.NaReal)
naeltype(::Type{RClass{:character}}) = sexp(Const.NaString)

# mirror src/main/arithmetic.c
function is_ieee_na(x::Float64)
    @static if reinterpret(UInt32,UInt64[1])[1] == 1  # little endian
        isnan(x) && reinterpret(UInt32,[x])[1] == 0x7a2
    else
        isnan(x) && reinterpret(UInt32,[x])[2] == 0x7a2
    end
end

"""
Check if a value corresponds to R's sentinel NA values.
These function should not be exported.
"""
isNA(x::ComplexF64) = is_ieee_na(real(x)) && is_ieee_na(imag(x))
isNA(x::Float64) = is_ieee_na(x)
isNA(x::Int32) = x == Const.NaInt
isNA(s::Ptr{CharSxp}) = s === sexp(Const.NaString)
isNA(s::Ptr{S}) where S<:VectorSxp = length(s) == 1 ? isNA(s[1]) : false
# all other values are consided as non-NA.
isNA(s::Any) = false

isna(s::Ptr{S}, i::Integer) where S<:VectorSxp = isNA(s[i])
isna(s::Ptr{S}) where S<:VectorSxp = reshape(BitArray([isNA(a) for a in s]), size(s))
"""
Check if the ith member of s coorespond to R's NA values.
"""
isna(r::RObject, i::Integer) = isna(r.p, i)
"""
Check if the members of a vector are NA values. Always return a BitArray.
"""
isna(r::RObject) = isna(r.p)

function anyna(s::Ptr{S}) where S<:VectorSxp
    for a in s
        if isNA(a)
            return true
        end
    end
    return false
end
"""
Check if there are any NA values in the vector.
"""
anyna(r::RObject{S}) where S<:VectorSxp = anyna(r.p)

# StrSxp
"""
Determines the encoding of the CharSxp. This is determined by the 'gp' part of the sxpinfo (this is the middle 16 bits).
 * 0x00_0002_00 (bit 1): set of bytes (no known encoding)
 * 0x00_0004_00 (bit 2): Latin-1
 * 0x00_0008_00 (bit 3): UTF-8
 * 0x00_0040_00 (bit 6): ASCII

We only support ASCII and UTF-8.
"""
function isascii(s::CharSxp)
    if s.head.info & 0x00_0040_00 != 0
        return true
    elseif s.head.info & 0x00_0008_00 != 0
        return false
    else
        error("Unsupported string type.")
    end
end
isascii(s::Ptr{CharSxp}) = isascii(unsafe_load(s))
isascii(r::RObject{CharSxp}) = isascii(sexp(r))

function isascii(s::Ptr{StrSxp})
    ind = true
    for c in s
        ind &= isNA(c) || isascii(c)
    end
    return ind
end
isascii(r::RObject{StrSxp}) = isascii(sexp(r))

# EnvSxp

function findVarInFrame(e::Ptr{EnvSxp}, s::Ptr{SymSxp})
    ccall((:Rf_findVarInFrame,libR),Ptr{UnknownSxp},(Ptr{EnvSxp},Ptr{SymSxp}),e,s)
end
findVarInFrame(e, s) = findVarInFrame(sexp(e), sexp(SymSxp, s))

function defineVar(s::Ptr{SymSxp}, v::Ptr{S}, e::Ptr{EnvSxp}) where S<:Sxp
    ccall((:Rf_defineVar,libR),Nothing,(Ptr{SymSxp},Ptr{S},Ptr{EnvSxp}),s,v,e)
    nothing
end
defineVar(s, v, p) = defineVar(sexp(SymSxp, s), sexp(v), sexp(p))

"extract the value of symbol s in the environment e"
function getindex(e::Ptr{EnvSxp},s::Ptr{SymSxp})
    v = findVarInFrame(e, s)
    v == sexp(Const.UnboundValue) && throw(BoundsError())
    sexp(v)
end
getindex(e::Ptr{EnvSxp},s) = getindex(e,sexp(SymSxp,s))
getindex(e::RObject{EnvSxp},s) = RObject(getindex(sexp(e),s))

"assign value v to symbol s in the environment e"
function setindex!(e::Ptr{EnvSxp},v::Ptr{S},s::Ptr{StrSxp}) where S<:Sxp
    # `Rf_defineVar` is unsafe to use if the binding is locked.
    # However, `setVarInFrame` is not exported. `base::assign` is
    # an avaliable alternative.
    rcall_p(Const.BaseNamespace["assign"], s, v, envir = e)
end
function setindex!(e::Ptr{EnvSxp},v,s)
    nprotect = 0
    try
        sv = protect(sexp(v))
        nprotect += 1
        ss = protect(sexp(RClass{:character},s))
        nprotect += 1
        setindex!(e,sv,ss)
    finally
        unprotect(nprotect)
    end
end
setindex!(e::RObject{EnvSxp}, v, s) = setindex!(sexp(e), v, s)

"""
    newEnvironment([env])

Create a new environment which extends environment `env` (`globalEnv` by default).
"""
function newEnvironment(env::Ptr{EnvSxp})
    ccall((:Rf_NewEnvironment,libR),Ptr{EnvSxp},
            (Ptr{NilSxp},Ptr{NilSxp},Ptr{EnvSxp}),sexp(Const.NilValue),sexp(Const.Const.NilValue),env)
end
newEnvironment(env::RObject{EnvSxp}) = newEnvironment(sexp(env))
newEnvironment() = newEnvironment(globalEnv)

"find namespace by name of the namespace, it is not error tolerant."
function findNamespace(str::String)
    ccall((:R_FindNamespace,libR),Ptr{EnvSxp}, (Ptr{StrSxp},), sexp(str))
end

"get namespace by name of the namespace. It is safer to be used than findNamespace as it checks bound."
getNamespace(str::String) = reval(rlang(RCall.Const.BaseNamespace["getNamespace"], str))

"Set the variable .Last.value to a given value"
function set_last_value(s::Ptr{S}) where S<:Sxp
    ccall((:SET_SYMVALUE,libR),Nothing,(Ptr{SymSxp},Ptr{UnknownSxp}),sexp(Const.LastvalueSymbol),s)
    nothing
end


text:
using Test, SimpleTensorNetworks.ContractionOrder
using LightGraphs
using LightGraphs: SimpleEdge
using SimpleTensorNetworks.ContractionOrder: abstract_contract_edge!, log2sumexp2

function get_shapes(edges, sizes, neighbors)
    log2_shapes = Vector{Float64}[]
    for i=1:length(neighbors)
        si = Float64[]
        for nb in neighbors[i]
            ie = findfirst(e->e == SimpleEdge(i=>nb) || e == SimpleEdge(nb=>i), edges)
            push!(si, sizes[ie])
        end
        push!(log2_shapes, si)
    end
    return log2_shapes
end

tolog2(shapes) = [log2.(x) for x in shapes]

@testset "abstract contract" begin
    es = [SimpleEdge(1, 2), SimpleEdge(2, 3), SimpleEdge(1, 3), SimpleEdge(3, 4)]
    log2_sizes = log2.([2, 3, 4, 5])
    neighbors = [[2,3], [3,1], [1,2,4], [3]]
    log2_shapes = get_shapes(es, log2_sizes, neighbors)
    strategy = "min_reduce"
    edge_info = edgeinfo(es, log2_shapes, neighbors, strategy)

    edge = SimpleEdge(2, 3)
    tc, sc = abstract_contract_edge!(edge, log2_shapes, edge_info,
        neighbors, strategy)
    @test edge_info.data[SimpleEdge(2, 4)] ≈ 8.0
    @test edge_info.data[SimpleEdge(1, 2)] ≈ 5.0
    @test log2_shapes == [[8], [8, 5], Int[], [5]] |> tolog2
    @test neighbors == [[2], [1, 4], Int[], [2]]
    @test tc ≈ log2(120)

    es = [SimpleEdge(1, 2), SimpleEdge(2, 3), SimpleEdge(1, 3), SimpleEdge(3, 4)]
    log2_sizes = log2.([2, 3, 4, 5])
    neighbors = [[2,3], [3,1], [1,2,4], [3]]
    log2_shapes = get_shapes(es, log2_sizes, neighbors)
    tn = TensorNetworkLayout(SimpleGraph(4, neighbors), log2_shapes)
    order = [SimpleEdge(1, 2),SimpleEdge(1, 3), SimpleEdge(1, 4)]
    tcs, scs = abstract_contract(tn, order)
    @test log2sumexp2(tcs) ≈ log2(89)
    @test maximum(scs) ≈ log2(12*5+12+5)
end

@testset "greedy order" begin
    g = random_regular_graph(10, 3)
    log2_shapes = [[2,2,2] for i=1:10] |> tolog2
    tn = TensorNetworkLayout(g, log2_shapes)
    tcs, scs, order = order_greedy(tn; strategy="min_dim")
    @test length(order) <= length(edges(g))
end

@testset "disconnected graph" begin
    n = 5
    g = SimpleGraph(n)
    add_edge!(g, 2, 3)
    tn = TensorNetworkLayout(g, [ones(degree(g, i)) for i=1:n])
    tc, sc, orders = order_greedy(tn)
    @test orders isa Vector
    tcs, scs = abstract_contract(tn, orders)
    @test maximum(tcs) == 1
end


text:
module PointNshoot
include("StarStuff.jl")

using .StarStuff
using GLMakie

"""
        singleParticleTree(particle::ParticleID,energy::Quantity,r0::Quantity,direction::Vector,gpFrac::Float64)

Build particle tree for a single particle
"""
function singleParticleTree(particle::ParticleID,energy::Quantity,r0::Quantity,direction::Vector,gpFrac::Float64)
    # TODO: allow fields as input
    # define electric and magnetic fields
    efield(x::Array{T}) where {T<:Number} = vec([0.0, 0.0, 0.0])
    bfield = Bfield.model[:jansson2012m]

    # generate initial particle tree
    initialPT = StarStuff.initParticles(
        Vector([particle]),
        Vector([upreferred(energy)]),
        Vector([upreferred(r0)]),
        Vector([direction]),
        gpFrac,
        bfield
    )

    # for construction ParticleTree and input to integrator
    params = Dict(
                  :t => 0.0,
                  :efield => efield,
                  :bfield => bfield
                 )

    # INTEGRATION
    # declare variable name for solution particle tree
    solPT = nothing
    # iterate over Timescales
    for i = 1:length(initialPT.nodes)
        # Build problem using boris algorithm integrator and initial Snapshot
        # of each Timescale
        prob = DiscreteProblem(rintegrator[:boris],
                               initialPT[i,1], (0.0,initialPT.nodes[i].dt*100000), params)
        println("test $i")
        println(typeof(initialPT[i,1]))
        println(typeof(prob))
        # Solve problem using time step specified in Timescale node
        solution = solve(prob, FunctionMap(), dt = initialPT.nodes[i].dt, progress = false)
        println("test $(i) end")
        # Construct a new Timescale out of the solution and add to solution particle tree
        solTS = construct(Timescale,solution.u,Float64[],solution.t,initialPT.nodes[i].dt)
        if solPT === nothing
            solPT = construct(ParticleTree, Vector([solTS]))
        else
            add_node!(solPT,solTS)
        end
    end
    return solPT
end

function best_units(x::Float64)
    # order of magnitude
    oom = mag10(x)
    units = nothing
    # TODO: find method to get strings from units
    if oom < 3
        units = (u"m","m")
    elseif oom < 6
        units = (u"km","km")
    elseif oom < 9
        units = (u"Mm","Mm")
    elseif oom < 14
        units = (u"AU","AU")
    elseif oom < 20
        units = (u"pc","pc")
    else
        units = (u"kpc","kpc")
    end
    return units
end

function show_bfield!(scene::Scene)
    println("Entered show_bfield!")
    # units for all distances
    units_distance = u"kpc"

    # number of points along radial coordinate
    nr = 10
    # innermost radius is 1 kpc IAW Jansson2012
    rmin = 1.0u"kpc"
    # take outermost radius from disk limits in Jansson field
    rmax = Bfield.janssonDict[:rdisklims][end]

    # number of points along z-coordinate
    nz = 20
    # set z-coordinate boundaries to be ±10 kpc for now to bring in some of the halo field
    zmin = -10.0u"kpc"
    zmax = 10.0u"kpc"

    # generate r and z coordinate arrays
    rr = LinRange(rmin,rmax,nr)
    zz = LinRange(zmin,zmax,nz)

    # array to store coordinates
    coords = nothing

    # number of points around innermost radius
    nmin = 5
    # number of points around outermost radius
    nmax = 20
    # generate cartesian coordinates and store in coords
    for (index,r) in enumerate(rr)
        # calculate number of points along this radius
        frac = (nmax-nmin)*index/nr
        n = nmin + Int(round(frac))
        println("Number of points around radius $index = ", n)
        # generate polar angle coordinates
        θ = LinRange(0.0,2*π,n)
        # transform to cartesian
        yy = r*sin.(θ)
        xx = r*cos.(θ)
        # make array for storing coordinates
        nxnynz = length(xx)*length(yy)*length(zz)
        tempcoords = zeros(Float64,(nxnynz,3))*units_distance
        # iterate over all combinations of generated coordinates to form grid
        ii = 1
        for xi in xx
            for yi in yy
                for zi in zz
                    tempcoords[ii,:] .= xi,yi,zi
                    ii+=1
                end
            end
        end
        # concatenate new coordinates
        if index == 1
            coords = tempcoords
        else
            coords = vcat(coords,tempcoords)
        end
    end
    # xx = nothing
    # yy = nothing
    # for (index,r) in enumerate(rr)
    #     # calculate number of points along this radius
    #     frac = (nmax-nmin)*index/nr
    #     n = nmin + Int(round(frac))
    #     println("Number of points around radius $index = ", n)
    #     # generate polar angle coordinates
    #     θ = LinRange(0.0,2*π,n)
    #     # transform to cartesian
    #     y = r*sin.(θ)
    #     x = r*cos.(θ)
    #     if index == 1
    #         yy = y
    #         xx = x
    #     else
    #         append!(yy,y)
    #         append!(xx,x)
    #     end
    # end

    println("Made all spatial coordinates. Calculating magnetic field.")
    # generate array to store magnetic field values
    # bunits = Unitful.unit(Bfield.jansson2012(coords[1,:])[1])
    # bb = zeros(Float64,size(coords))*bunits
    # calculate magnetic field values
    # for i = 1:size(coords)[1]
    #     bb[i,:] .= Bfield.jansson2012(coords[i,:])
    # end
    # println("Finished calculating magnetic field")
    janssonx(x,y,z) = Bfield.jansson2012kpc([x,y,z])[1]*1.0e10
    janssony(x,y,z) = Bfield.jansson2012kpc([x,y,z])[1]*1.0e10
    janssonz(x,y,z) = Bfield.jansson2012kpc([x,y,z])[1]*1.0e10

    coords = ustrip.(units_distance, coords)

    bz = [janssonz(v...) for v ∈ eachrow(coords)]
    println(bz[1:10])
    println(min(bz...))
    println(max(bz...))

    scatter!(scene,coords[:,1],coords[:,2],coords[:,3],color=bz*10,markersize=100)
    # volume!(scene,coords[:,1],coords[:,2],coords[:,3],jansson1,algorithm=:mip)
    # volume!(scene,rand(32,32,32), algorithm=:mip)
end

function pnsScene!(scene::Scene,pt::ParticleTree;filename::Union{String,Nothing}=nothing)
    # Extract coordinates from particle tree
    times, traj = getParticleTrajectory(pt,1,1,1)
    # println("dt = ", times[51]-times[50])
    # println("steps = ",length(times))
    # println("max time = ", times[end])

    xx = [c.x[1][1] for c in traj]
    yy = [c.x[1][2] for c in traj]
    zz = [c.x[1][3] for c in traj]
    coords = [xx,yy,zz]

    #out_df = DataFrame(t=times,x=xx,y=yy,z=zz)
    #CSV.write()

    # Shift coordinates so they're ceneterd at zero
    μcoords = mean.(coords)
    shiftedCoords = [coords[i] .- μcoords[i] for i=1:3]
    localOrigin = μcoords*u"m"
    localOrigin = uconvert.(u"kpc",localOrigin)

    # Calculate bounds for plotting
    xtreme,ystreme,zstreme = extrema.(shiftedCoords)
    mincoord = min(xtreme[1],ystreme[1],zstreme[1])
    maxcoord = max(xtreme[2],ystreme[2],zstreme[2])

    bounds = FRect3D([mincoord,mincoord,mincoord],[maxcoord,maxcoord,maxcoord])

    # 3D SCENE
    makelabels(xs,tickunits) = ["$(@sprintf("%.2e",ustrip(tickunits,x*u"m")))" for x in xs]
    maxticks = collect(LinRange(mincoord,0.0,5))
    append!(maxticks,collect(LinRange(0.0,maxcoord,5)[2:end]))
    tickunits,axesunits = best_units(maxcoord-mincoord)
    ticklabels = makelabels(maxticks,tickunits)
    #	Plot coordinates and labels
    lines!(scene, shiftedCoords..., scale_plot=false,limits = bounds)
    xticks!(scene; xtickrange=maxticks, xticklabels=ticklabels)
    yticks!(scene; ytickrange=maxticks, yticklabels=ticklabels)
    zticks!(scene; ztickrange=maxticks, zticklabels=ticklabels)
    xlabel!(scene, "X ($axesunits)")
    ylabel!(scene, "Y ($axesunits)")
    zlabel!(scene, "Z ($axesunits)")
    AbstractPlotting.update!(scene)

    return scene
end

"""
        run()

Run PointNshoot script
"""
function run()
    GLMakie.activate!()
    scene = Scene()
    show_bfield!(scene)
    return scene
end

export run

end # PointNshoot

import .PointNshoot
PointNshoot.run()


text: Dict{String,Any}("key" => Dict{String,Any}("value" => "value","type" => "string"))

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
using OffsetArrays, Test
@testset "UpperboundedIntegers" begin
  function ubsum(x)
    @assert firstindex(x) == 0
    r = LoopVectorization.CloseOpen(LoopVectorization.UpperBoundedInteger(length(x), LoopVectorization.StaticInt(15)))
    s = zero(eltype(x))
    @turbo for i ∈ r
      s += x[i]
    end
    s
  end
  function ubdouble!(y,x)
    @assert firstindex(x) == 0
    r = LoopVectorization.CloseOpen(LoopVectorization.UpperBoundedInteger(length(x), LoopVectorization.StaticInt(15)))
    @turbo for i ∈ r
      y[i] = 2*x[i]
    end
    y
  end
  for l ∈ 1:15
    x = OffsetVector(rand(l), -1)
    @test ubsum(x) ≈ sum(x)
    @test ubdouble!(similar(x), x) == x .* 2
  end
  for l ∈ 70:75
    # If the actual value is greater than the ubber bound, as is the case here,
    # (upper bound = 15, actual length = 70:75), then the result is undefined.
    # Currently, it'll evaluate to the actual value with heavily unrolled reductions
    x = OffsetVector(rand(l), -1)
    @test ubsum(x) ≈ ubsum(x)
    xs = similar(x)
    # both otherwise, the behavior is undefined. It'll evaluate at most
    # UF * W * cld(UB, W)
    # iterations, where UF is the unroll factor, W the vector width,
    # and N the actual length.
    # It may evaluate less than the upper bound, depending on the mask's value.
    # We check that the first few match
    @test @view(ubdouble!(xs, x)[begin:3]) == @view(x[begin:3]) .* 2;
    LoopVectorization.register_size() > 16 && (@test xs[end] ≠ 2x[end])
  end
end


text:
__precompile__()
module UnityGymWrapper
using PyCall
include("Spaces.jl")

export Environment, step!, reset!, close!

const pymlagents = PyNULL()
const pygym = PyNULL()

function __init__()
    copy!(pymlagents, pyimport("mlagents_envs.environment"))
    copy!(pygym, pyimport("gym_unity.envs"))
end

struct Environment{A}
    pyenv::PyObject
    pyreset::PyObject
    pystep::PyObject
    pyclose::PyObject
    actionspace::ActionSpace
    observationspace::Union{ObservationSpace, Vector{ObservationSpace}}

    function Environment(path;
            nographics=true, usevisual=false, uint8visual=false, multipleobs=false, kwargs...)
        args = collect(Iterators.flatten(Tuple{String,String}[("-$(arg[1])", "$(arg[2])") for arg in kwargs]))
        env = pymlagents.UnityEnvironment(path, no_graphics=nographics, args=args)
        try
            gym = pygym.UnityToGymWrapper(
                env, use_visual=usevisual, uint8_visual=uint8visual, allow_multiple_visual_obs=multipleobs
            )
            as = ActionSpace(gym."action_space")
            os = ObservationSpace(gym."observation_space")
            return new{eltype(as)}(gym, gym."reset", gym."step", gym."close", as, os)
        catch
            pycall(env."close", Nothing)
            rethrow()
        end
    end
end

reset!(env) = pycall(env.pyreset, PyArray)

step!(env) = step!(env, sample(env.actionspace))
step!(env, action::Integer) = step!(env, [action])
function step!(env::Environment{A}, action::AbstractArray{A}) where A
    if A <: Integer
        action -= ones(A, size(action)) # "Python conversion"
    end
    return pycall(env.pystep, Tuple{PyArray, Float32, Bool, PyObject}, PyReverseDims(action))
end

close!(env) = pycall(env.pyclose, Nothing)

end


text:
"""
    ClimaTimeSteppers

Ordinary differential equation solvers

JuliaDiffEq terminology:

* _Function_: the right-hand side function df/dt.
  * by default, a function gets wrapped in an `ODEFunction`
  * define new `IncrementingODEFunction` to support incrementing function calls.

* _Problem_: Function, initial u, time span, parameters and options

  du/dt = f(u,p,t) = fL(u,p,t)  + fR(u,p,t)

  fR(u,p,t) == f(u.p,t) - fL(u,p,t)
  fL(u,_,_) == A*u for some `A` (matrix free)

  SplitODEProlem(fL, fR)

  * `ODEProblem` from OrdinaryDiffEq.jl
    - use `jac` option to `ODEFunction` for linear + full IMEX (https://docs.sciml.ai/latest/features/performance_overloads/#ode_explicit_jac-1)
  * `SplitODEProblem` for linear + remainder IMEX
  * `MultirateODEProblem` for true multirate

* _Algorithm_: small objects (often singleton) which indicate what algorithm + options (e.g. linear solver type)
  * define new abstract `DistributedODEAlgorithm`, algorithms in this pacakge will be subtypes of this
  * define new `Multirate` for multirate solvers

* _Integrator_: contains everything necessary to solve. Used as:

  * define new `DistributedODEIntegrator` for solvers in this package

      init(prob, alg, options...) => integrator
      step!(int) => runs single step
      solve!(int) => runs it to end
      solve(prob, alg, options...) => init + solve!

* _Solution_ (not implemented): contains the "solution" to the ODE.

"""
module ClimaTimeSteppers

using KernelAbstractions
using KernelAbstractions.Extras: @unroll
using LinearAlgebra
using StaticArrays
using CUDA
using MPI

array_device(::Union{Array, SArray, MArray}) = CPU()
array_device(::CuArray) = CUDADevice()
realview(x::Union{Array, SArray, MArray}) = x
realview(x::CuArray) = x

import DiffEqBase, SciMLBase, LinearAlgebra, DiffEqCallbacks

include("functions.jl")
include("operators.jl")

abstract type DistributedODEAlgorithm <: DiffEqBase.AbstractODEAlgorithm
end

include("integrators.jl")

# Include concrete implementations
include("solvers/multirate.jl")
include("solvers/lsrk.jl")
include("solvers/ssprk.jl")
include("solvers/ark.jl")
include("solvers/mis.jl")
include("solvers/wickerskamarock.jl")
include("solvers/rosenbrock.jl")

include("callbacks.jl")

end


text:
module TelegraphWalker

using Reexport

@reexport using Random
@reexport using Distributions
@reexport using DataFrames
@reexport using Chain
@reexport using RandomWalker
@reexport using Telegraph

export
	MeanSquaredDisplacement,
	DisplacementPositions,
	TaxiCabDistance,
	TaxiCabDistanceXY,
	taxicabdistanceXY,
	taxicabdistance,
	periodicdomainshift,
	taxicabdistance,
	taxicabdistanceXY,
	walker_distance,
	random_walk_telegraph,
	random_walk_telegraph_no_df,
	dₜ,
	tele_rw_inf_domain

include("telegraph_walker_structs.jl")
include("taxicab_distances.jl")
include("telegraph_walker_couple.jl")

end


text:
module test_lupin_arguments

using Test
using Lupin

f(x) = 2x
g(; x) = 2x

@test apply(f, Argument(3))     == f(3)
@test apply(g, Argument(; x=3)) == g(; x=3)

end # module test_lupin_arguments


text:
using Gadfly

function main(dataPath)
	lines = parse(readall(pipeline(`cat $[dataPath]ratings.csv`, `wc -l`)))
	@time ratings = readdlm(dataPath * "ratings.csv", ',', header=true, dims=(lines, 4))
	ratingsHeader = ratings[2]
	ratings = ratings[1]

	lines = parse(readall(pipeline(`cat $[dataPath]movies.csv`, `wc -l`)))
	@time movies = readdlm(dataPath * "movies.csv", ',', header=true, dims=(lines, 3))
	moviesHeader = movies[2]
	movies = movies[1]

	lines = parse(readall(pipeline(`cat $[dataPath]links.csv`, `wc -l`)))
	@time links = readdlm(dataPath * "links.csv", ',', header=true, dims=(lines, 3))
	linksHeader = links[2]
	links = links[1]

	lines = parse(readall(pipeline(`cat $[dataPath]tags.csv`, `wc -l`)))
	@time tags = readdlm(dataPath * "tags.csv", ',', header=true, dims=(lines, 4))
	tagsHeader = tags[2]
	tags = tags[1]

	newMovies = deepcopy(movies[:,1])
	newRatings = deepcopy(ratings[:,2])
	newLinks = deepcopy(links[:,1])
	newTags = deepcopy(tags[:,2])

	@printf "Finished loading, going into loop\n"
	@printf "%s/%s" 0 length(movies[:,1])

	for i = 1:length(movies[:, 1])
		if i % 1_000 == 0
			@printf "\n%s/%s" i length(movies[:,1])
		elseif i % 100 == 0
			@printf "."
		end

		if movies[i,1] != i
			newMovies[i] = i
			newRatings[ratings[:,2] .== movies[i,1]] = i
			newLinks[links[:,1] .== movies[i,1]] = i
			newTags[tags[:,2] .== movies[i,1]] = i
		end
	end

	@printf "\nFinished loop, moving changes into main array\n"

	movies[:,1] = newMovies
	ratings[:,2] = newRatings
	links[:,1] = newLinks
	tags[:,2] = newTags

	movies = [moviesHeader; movies]
	ratings = [ratingsHeader; ratings]
	links = [linksHeader; links]
	tags = [tagsHeader; tags]

	@printf "Finished moving, begin saving\n"

	if !isdir(dataPath * "fixed")
		mkdir(dataPath * "fixed")
	end
	writecsv(dataPath * "fixed/movies.csv", movies)
	writecsv(dataPath * "fixed/ratings.csv", ratings)
	writecsv(dataPath * "fixed/links.csv", links)
	writecsv(dataPath * "fixed/tags.csv", tags)

	@printf "Done :)\n"
end

main("ml-latest/fixed/")


text:
export single_point_crossover,
    random_node_crossover,
    aligned_node_crossover,
    proportional_crossover,
    output_graph_crossover,
    subgraph_crossover,
    crossover

"single point crossover, genes from p1 up to a random point, then genes from p2"
function single_point_crossover(cfg::NamedTuple, c1::CGPInd, c2::CGPInd)
    cpoint = rand(2:(min(length(c1.chromosome), length(c2.chromosome))-2))
    if rand(Bool)
        ngenes = deepcopy([c1.chromosome[1:cpoint]; c2.chromosome[(cpoint+1):end]])
        return CGPInd(cfg, ngenes)
    else
        ngenes = deepcopy([c2.chromosome[1:cpoint]; c1.chromosome[(cpoint+1):end]])
        return CGPInd(cfg, ngenes)
    end
end

"returns random input genes from c1 and c2 equally, not a full crossover operator"
function random_inputs(c1::CGPInd, c2::CGPInd)
    child_inputs = zeros(c1.n_in)
    parent = rand(Bool, c1.n_in)
    child_inputs[parent] = c1.chromosome[1:c1.n_in][parent]
    child_inputs[.~parent] = c2.chromosome[1:c1.n_in][.~parent]
    child_inputs
end

"returns random output genes from c1 and c2 equally, not a full crossover operator"
function random_outputs(c1::CGPInd, c2::CGPInd)
    child_outputs = zeros(c1.n_out)
    parent = rand(Bool, c1.n_out)
    s = length(c1.chromosome) - c1.n_out + 1
    child_outputs[parent] = c1.chromosome[s:end][parent]
    child_outputs[.~parent] = c2.chromosome[s:end][.~parent]
    child_outputs
end

"take random nodes from each parent equally, up to the size of the smaller parent"
function random_node_crossover(cfg::NamedTuple, c1::CGPInd, c2::CGPInd)
    l = min(length(c1.chromosome), length(c2.chromosome)) - max(c1.n_out, c2.n_out)
    p = rand(Bool, Int64(round(l / 3)))
    p_nodes = repeat(p, inner=3)
    append!(p_nodes, rand(Bool, c1.n_out))
    genes = rand(length(c1.chromosome))
    genes[p_nodes] = c1.chromosome[p_nodes]
    genes[.~p_nodes] = c2.chromosome[.~p_nodes]
    CGPInd(cfg, genes)
end

"align nodes based on position, then take from each parent equally, only works on PCGP"
function aligned_node_crossover(cfg::NamedTuple, c1::CGPInd, c2::CGPInd)
    p1_pos = get_positions(c1)
    p2_pos = get_positions(c2)
    min_nodes = min(length(c1.nodes)-c1.n_in, length(c2.nodes)-c2.n_in)
    p1_inds = c1.n_in + collect(1:min_nodes)
    p1_nodes = Array{Int64}(undef, 0)
    p2_nodes = Array{Int64}(undef, 0)
    for node in c1.n_in+(1:min_nodes)
        if rand() < 0.5
            i = indmin(abs.(p1_pos[p1_inds] - p2_pos[node]))
            append!(p1_nodes, [p1_inds[i]])
            p1_inds = deleteat!(p1_inds, i)
        else
            append!(p2_nodes, [node])
        end
    end
    genes = [random_inputs(c1, c2); random_outputs(c1, c2); get_genes(
        c1, p1_nodes); get_genes(c2, p2_nodes)]
    CGPInd(cfg, genes)
end

"choose a point between the genes from the two parents for a child individual"
function proportional_crossover(cfg::NamedTuple, c1::CGPInd, c2::CGPInd)
    genes = Array{Float64}(undef, 0)
    if length(c1.chromosome) == length(c2.chromosome)
        r = rand(length(c1.chromosome))
        genes = ((1 .- r) .* c1.chromosome) .+ (r .* c2.chromosome)
    elseif length(c1.chromosome) > length(c2.chromosome)
        r = rand(length(c2.chromosome))
        genes = ((1 .- r) .* c1.chromosome[1:length(c2.chromosome)]) .+ (r .* c2.chromosome)
        genes = [genes; c1.chromosome[(length(c2.chromosome)+1):end]]
    else
        r = rand(length(c1.chromosome))
        genes = ((1 .- r) .* c1.chromosome) .+ (r .* c2.chromosome[1:length(c1.chromosome)])
        genes = [genes; c2.chromosome[(length(c1.chromosome)+1):end]]
    end
    CGPInd(cfg, genes)
end

"crossover with a list of nodes"
function node_crossover(cfg::NamedTuple, c1::CGPInd, c2::CGPInd,
                        c1_nodes::Array{Int16}, c2_nodes::Array{Int16},
                        output_genes::AbstractArray = Array{Float64}(undef, 0))
    s = length(c1.chromosome) - c1.n_out
    if length(output_genes) == 0
        for o in 1:c1.n_out
            g = rand() < 0.5 ? c1.chromosome[s+o] : c2.chromosome[s+o]
            push!(output_genes, g)
        end
    end
    genes = Array{Float64}(undef, 0)
    for i in (c1.n_in+1):length(c1.nodes)
        if i in c1_nodes
            if i in c2_nodes
                g = rand() < 0.5 ? get_genes(c1, i) : get_genes(c2, i)
                append!(genes, g)
            else
                append!(genes, get_genes(c1, i))
            end
        else
            if i in c2_nodes
                append!(genes, get_genes(c2, i))
            else
                append!(genes, rand(3))
            end
        end
    end
    append!(genes, output_genes)
    CGPInd(cfg, genes)
end

"split outputs equally between parents, then construct a child from their corresponding input graphs"
function output_graph_crossover(cfg::NamedTuple, c1::CGPInd, c2::CGPInd)
    p1_nodes = Array{Int16}(undef, 0)
    p2_nodes = Array{Int16}(undef, 0)
    output_genes = Array{Float64}(undef, 0)
    s = length(c1.chromosome) - c1.n_out
    for output in 1:c1.n_out
        if rand() < 0.5
            append!(p1_nodes, get_output_trace(c1, output))
            push!(output_genes, c1.chromosome[s+output])
        else
            append!(p2_nodes, get_output_trace(c2, output))
            push!(output_genes, c2.chromosome[s+output])
        end
    end
    c1_nodes = sort(unique(p1_nodes))
    c2_nodes = sort(unique(p2_nodes))
    node_crossover(cfg, c1, c2, c1_nodes, c2_nodes, output_genes)
end

"Take subgraphs from both parents equally, adding all nodes of the chosen subgraphs to the child"
function subgraph_crossover(cfg::NamedTuple, c1::CGPInd, c2::CGPInd)
    fc1 = forward_connections(c1)
    fc2 = forward_connections(c2)
    c1_nodes = []; c2_nodes = []
    for i in (c1.n_in+1):min(length(c1.nodes), length(c2.nodes))
        if rand(Bool)
            append!(c1_nodes, fc1[i])
        else
            append!(c2_nodes, fc2[i])
        end
    end
    if length(c1.nodes) < length(c2.nodes)
        for i in length(c1.nodes):length(c2.nodes)
            if rand(Bool)
                append!(c2_nodes, fc2[i])
            end
        end
    elseif length(c2.nodes) < length(c1.nodes)
        for i in length(c2.nodes):length(c1.nodes)
            if rand(Bool)
                append!(c1_nodes, fc1[i])
            end
        end
    end
    c1_nodes = Array{Int16}(unique(intersect(collect((c1.n_in+1):length(c1.nodes)), c1_nodes)))
    c2_nodes = Array{Int16}(unique(intersect(collect((c2.n_in+1):length(c2.nodes)), c2_nodes)))
    node_crossover(cfg, c1, c2, c1_nodes, c2_nodes)
end


text:
"""
    treeda_xgb(X, y;
        subsample = 1, colsample_bytree = 1, colsample_bynode = 1,
        max_depth = 6, min_child_weight = 5,
        lambda = 0, verbose = false, kwargs...)
Discrimination tree with XGBoost.
* `X` : X-data (n obs., p variables).
* `y` : Univariate Y-data (n obs.).
* `subsample` : Proportion of rows sampled in `X`
    for building the tree.
* `colsample_bytree` : Proportion of columns sampled in `X`
    for building the tree.
* `colsample_bynode` : Proportion of columns sampled at each node
    in the columns selected for the tree.
* `max_depth` : Maximum depth of the tree.
* `min_child_weight` : Minimum nb. observations that each leaf
    needs to have.
* `lambda` : L2 regularization term on weights.
    Increasing this value will make model more conservative.
* `verbose` : If true, fitting information are printed.
* `kwargs` : Optional named arguments to pass in function `xgboost`
    of `XGBoost.jl` (https://xgboost.readthedocs.io/en/latest/parameter.html).

The function builds a single tree using package `XGboost.jl'.

The sampling of the observations and variables are without replacement.

## References

XGBoost.jl
https://github.com/dmlc/XGBoost.jl

Breiman, L., Friedman, J. H., Olshen, R. A., and Stone, C. J. Classification
And Regression Trees. Chapman & Hall, 1984.

Gey, S., 2002. Bornes de risque, détection de ruptures, boosting :
trois thèmes statistiques autour de CART en régression (These de doctorat).
Paris 11. http://www.theses.fr/2002PA112245
"""
function treeda_xgb(X, y;
        subsample = 1, colsample_bytree = 1, colsample_bynode = 1,
        max_depth = 6, min_child_weight = 5,
        lambda = 0, verbose = false, kwargs...)
    X = ensure_mat(X)
    y = vec(y)
    p = size(X, 2)
    ztab = tab(y)
    y_num = recodcat2num(y; start = 0)
    num_class = length(ztab.keys)
    num_round = 1
    fm = xgboost(X, num_round; label = y_num,
        seed = Int64(round(rand(1)[1] * 1e5)),
        booster = :gbtree,
        tree_method = :auto,
        eta = 1, # learning rate
        subsample = subsample,
        colsample_bytree = colsample_bytree, colsample_bylevel = 1,
        colsample_bynode = colsample_bynode,
        max_depth = max_depth, min_child_weight = min_child_weight,
        lambda = lambda,
        objective = "multi:softmax",
        num_class = num_class,
        silent = !verbose, kwargs...)
    featur = collect(1:p)
    TreedaXgb(fm, featur, ztab.keys, ztab.vals)
end

"""
    rfda_xgb(X, y; rep = 50,
        subsample = .7,
        colsample_bytree = 1, colsample_bynode = 1/3,
        max_depth = 6, min_child_weight = 5,
        lambda = 0, verbose = false, kwargs...)
Random forest discrimination with XGBoost.
* `X` : X-data (n obs., p variables).
* `y` : Univariate Y-data (n obs.).
* `rep` : Nb. trees to build in the forest.
* `subsample` : Proportion of rows sampled in `X`
    for building each tree.
* `colsample_bytree` : Proportion of columns sampled in `X` at each tree.
* `colsample_bynode` : Proportion of columns sampled at each node
    in the columns selected for the tree.
* `max_depth` : Maximum depth of the trees.
* `min_child_weight` : Minimum nb. observations that each leaf
    needs to have.
* `lambda` : L2 regularization term on weights.
    Increasing this value will make model more conservative.
* `verbose` : If true, fitting information are printed.
* `kwargs` : Optional named arguments to pass in function `xgboost`
    of `XGBoost.jl` (https://xgboost.readthedocs.io/en/latest/parameter.html).

The function uses package `XGboost.jl' to build the forest.
See https://xgboost.readthedocs.io/en/latest/tutorials/rf.html.

## References

XGBoost.jl
https://github.com/dmlc/XGBoost.jl

Breiman, L., 1996. Bagging predictors. Mach Learn 24, 123–140.
https://doi.org/10.1007/BF00058655

Breiman, L., 2001. Random Forests. Machine Learning 45, 5–32.
https://doi.org/10.1023/A:1010933404324

Genuer, R., 2010. Forêts aléatoires : aspects théoriques,
sélection de variables et applications. PhD Thesis. Université Paris Sud - Paris XI.

Gey, S., 2002. Bornes de risque, détection de ruptures, boosting :
trois thèmes statistiques autour de CART en régression (These de doctorat).
Paris 11. http://www.theses.fr/2002PA112245
"""
function rfda_xgb(X, y; rep = 50,
        subsample = .7,
        colsample_bytree = 1, colsample_bynode = 1/3,
        max_depth = 6, min_child_weight = 5,
        lambda = 0, verbose = false, kwargs...)
    X = ensure_mat(X)
    y = vec(y)
    p = size(X, 2)
    ztab = tab(y)
    y_num = recodcat2num(y; start = 0)
    num_class = length(ztab.keys)
    num_round = 1
    fm = xgboost(X, num_round; label = y_num,
        seed = Int64(round(rand(1)[1] * 1e5)),
        booster = :gbtree,
        tree_method = :auto,
        num_parallel_tree = rep,
        eta = 1, # learning rate
        subsample = subsample,
        colsample_bytree = colsample_bytree, colsample_bylevel = 1,
        colsample_bynode = colsample_bynode,
        max_depth = max_depth, min_child_weight = min_child_weight,
        lambda = lambda,
        objective = "multi:softmax",
        num_class = num_class,
        silent = !verbose, kwargs...)
    featur = collect(1:p)
    TreedaXgb(fm, featur, ztab.keys, ztab.vals)
end

"""
    xgboostda(X, y; rep = 50, eta = .3,
        subsample = .7, colsample_bytree = 1, colsample_bynode = 1/3,
        max_depth = 6, min_child_weight = 5,
        lambda = 1, verbose = false, kwargs...)
XGBoost discrimination.
* `X` : X-data (n obs., p variables).
* `y` : Univariate Y-data (n obs.).
* `rep` : Nb. trees to build.
* `eta` : Learning rate ([0, 1]).
* `subsample` : Proportion of rows sampled in `X`
    for building each tree.
* `colsample_bytree` : Proportion of columns sampled in `X` at each tree.
* `colsample_bynode` : Proportion of columns sampled at each node
    in the columns selected for the tree.
* `max_depth` : Maximum depth of the trees.
* `min_child_weight` : Minimum nb. observations that each leaf
    needs to have.
* `lambda` : L2 regularization term on weights.
    Increasing this value will make model more conservative.
* `verbose` : If true, fitting information are printed.
* `kwargs` : Optional named arguments to pass in function `xgboost`
    of `XGBoost.jl` (https://xgboost.readthedocs.io/en/latest/parameter.html).

The function uses package `XGboost.jl' to build the forest.
See https://xgboost.readthedocs.io/en/latest/tutorials/rf.html.

## References

XGBoost
https://xgboost.readthedocs.io/en/latest/index.html

XGBoost.jl
https://github.com/dmlc/XGBoost.jl
"""
function xgboostda(X, y; rep = 50, eta = .3,
        subsample = .7, colsample_bytree = 1, colsample_bynode = 1/3,
        max_depth = 6, min_child_weight = 5,
        lambda = 1, verbose = false, kwargs...)
    X = ensure_mat(X)
    y = vec(y)
    p = size(X, 2)
    ztab = tab(y)
    y_num = recodcat2num(y; start = 0)
    num_class = length(ztab.keys)
    num_round = rep
    fm = xgboost(X, num_round; label = y_num,
        seed = Int64(round(rand(1)[1] * 1e5)),
        booster = :gbtree,
        tree_method = :auto,
        num_parallel_tree = 1,
        eta = eta,
        subsample = subsample,
        colsample_bytree = colsample_bytree,
        colsample_bylevel = 1, colsample_bynode = colsample_bynode,
        max_depth = max_depth, min_child_weight = min_child_weight,
        objective = "multi:softmax",
        num_class = num_class,
        lambda = lambda,
        silent = !verbose, kwargs...)
    featur = collect(1:p)
    TreedaXgb(fm, featur, ztab.keys, ztab.vals)
end

"""
    predict(object::TreedaXgb, X)
Compute Y-predictions from a fitted model.
* `object` : The fitted model.
* `X` : X-data for which predictions are computed.
"""
function predict(object::TreedaXgb, X)
    X = ensure_mat(X)
    m = size(X, 1)
    pred = XGBoost.predict(object.fm, X)
    pred = replacebylev(pred, object.lev)
    pred = reshape(pred, m, 1) ;
    (pred = pred,)
end


text:
import MPI
using Test
using Harlequin

include("test_mapmaking.jl")


text:
# Dummy row
struct SourceRow{Source}
    source::Source
    table_name::Symbol
end

# The model row pass will build and modify a model of a row
function model_row(call)
    arguments, keywords = split_call(call)
    model_row_call(arguments...; keywords...)
end

# Very few functions modify the row model, so leave the row unchanged by default
function model_row_call(arbitrary_function, iterator, arguments...)
    model_row(iterator)
end

# Create a model row when getproperty is called on a Database
function get_column(source_row, column_name)
    SourceCode(source_row.source,
        Expr(:call, getproperty, source_row, column_name)
    )
end

function model_row_call(::typeof(getproperty), source_tables::Database, table_name)
    source = get_source(source_tables)
    column_names = get_column_names(source, table_name)
    source_row = SourceRow(source, table_name)
    NamedTuple{column_names}(
        map(column_names) do column_name
        get_column(source_row, column_name)
    end
    )
end

# Map is the only function which directly modifies the model row
function model_row_call(::typeof(QueryOperators.map), iterator, call, call_expression)
    call(model_row(iterator))
end

# Grouped rows have their own dedicated model type
struct GroupRow{Group,Row}
    group::Group
    row::Row
end

function QueryOperators.key(group_of_rows::GroupRow)
    getfield(group_of_rows, :group)
end

function getproperty(group_of_rows::GroupRow, column_name::Symbol)
    getproperty(getfield(group_of_rows, :row), column_name)
end

function length(group_of_rows::GroupRow)
    length(first(getfield(group_of_rows, :row)))
end

function model_row_call(::typeof(QueryOperators.groupby), ungrouped, group_function, group_function_expression, map_selector, map_function_expression)
    model = map_selector(model_row(ungrouped))
    GroupRow(group_function(model), model)
end


text:
module TreeTools

using FastaIO
using JSON
using Dates
## Includes
include("objects.jl")
include("objectsmethods.jl")
include("mutations.jl")
include("prunegraft.jl")
include("datamethods.jl")
include("reading.jl")
include("writing.jl")
include("misc.jl")
include("lbi.jl")

end

## Todo
# the child field of `TreeNode` should be a set and not an array since ordering is not relevant? Howver it makes accessing more difficult. For now, giving up on this idea since I do not benefit from `Set` specific function implemented in Julia: they ultimately fall back to `===` which will consider equal nodes to be different.

max_stars_repo_name: PierreBarrat/FluPredictibility

max_issues_repo_name: PierreBarrat/FluPredictibility

max_forks_repo_name: PierreBarrat/FluPredictibility


text:
include("src/SSPS.jl")

using .SSPS

julia_main()

max_stars_repo_name: gitter-lab/ssps

max_issues_repo_name: gitter-lab/ssps

max_forks_repo_name: gitter-lab/ssps


text:
abstract type Material end

function emit(material::Material, u::AbstractFloat, v::AbstractFloat, point::AbstractVec)::AbstractVec
    Vec(0.0, 0.0, 0.0)
end

max_stars_repo_name: goooxu/JuliaRayTracing

max_issues_repo_name: goooxu/JuliaRayTracing

max_forks_repo_name: goooxu/JuliaRayTracing


text:
using ContinuousPOMDPTreeSearchExperiments
using Base.Test

using POMDPModels
using POMDPToolbox
using LaserTag
using ParticleFilters
using QMDP

srand(4)

sld = SimpleLightDark()
@test isterminal(sld, sld.radius+1)
p = solve(LDHSolver(), sld)
filter = SIRParticleFilter(sld, 1000)
for (s, b, a, r, sp, o) in stepthrough(sld, p, filter, "sbarspo", max_steps=100)
    @show (s, a, r, sp, o)
    @show mean(b)
end

qp = solve(QMDPSolver(), sld, verbose=true)
filter = SIRParticleFilter(sld, 1000)
for (s, b, a, r, sp, o) in stepthrough(sld, qp, filter, "sbarspo", max_steps=100)
    @show (s, a, r, sp, o)
    @show mean(b)
end

m = LQG1D(1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 3)
filter = Kalman1D(m)
policy = LinearFeedback(1/2)
rsum = 0.0
for (s, b, a, r, sp, o) in stepthrough(m, policy, filter, "sbarspo", max_steps=100)
    @show (s, a, r, sp, o)
    @show mean(b)
    rsum += r
end
@show rsum


text:
"Check if formula contains a predicate name."
has_pred(formula::Const, pred_names) = formula.name in pred_names
has_pred(formula::Var, pred_names) = false
has_pred(formula::Compound, pred_names) = formula.name in pred_names ||
    any(has_pred(f, pred_names) for f in formula.args)
has_pred(formula::Term, domain::Domain) =
    has_pred(formula, keys(domain.predicates))

"Check if formula contains a numeric fluent."
has_fluent(formula::Term, state::State) =
    has_pred(formula, keys(state.fluents))
has_fluent(formula::Term, domain::Domain) =
    has_pred(formula, keys(domain.functions))

"Check if formula contains an axiom reference."
has_axiom(formula::Term, domain::Domain) = length(domain.axioms) > 0 &&
    has_pred(formula, (ax.head.name for ax in domain.axioms))

"Check if formula contains a universal or existential quantifier."
has_quantifier(formula::Term) =
    has_pred(formula, Set(Symbol[:forall, :exists]))

"Get subscript from State"
get_value(state::State, key::String) = state[key]

"Check if term is in State"
has_term_in_state(domain::Domain, state::State, term::Term) =
    return state[domain, term]

flatten_goal(problem::Problem) =
    return flatten_conjs(problem.goal)

global g_preconditions = Dict{Symbol,Vector{Vector{Term}}}()

function cache_global_preconditions(domain::Domain)
    for (act_name, act_def) in domain.actions
        conds = filter_negative_preconds(act_def)
        g_preconditions[act_name] = conds
    end
end

"Filter out negative preconditions"
function filter_negative_preconds(action_def::Action)
    conds = get_preconditions(action_def; converter=to_dnf)
    conds = [c.args for c in conds.args]
    for c in conds filter!(t -> t.name != :not, c) end
    return conds
end

"Compute axioms without negative literals"
function compute_hsp_axioms(domain::Domain)
    axioms = regularize_clauses(domain.axioms) # Regularize domain axioms
    axioms = [Clause(ax.head, [t for t in ax.body if t.name != :not])
              for ax in axioms] # Remove negative literals
    domain.axioms = Clause[]
    return domain, axioms
end

"Create state class"
function create_stat(terms::Vector{<:Term}, types::Vector{<:Term}=Term[])
    return State(types, facts, Dict{Symbol, Any}())
end

"Compute costs of one-step derivation of domain axioms"
function compute_costs_one_step_derivation(facts::Set{Term}, fact_costs::Dict{Any,Any}, ax::Clause, heur::String)
    facts_costs = convert_dict_any_to_term(fact_costs)
    _, subst = resolve(ax.body, [Clause(f, []) for f in facts])
    for s in subst
        body = [substitute(t, s) for t in ax.body]
        if heur == "delete_relaxation/h_add"
            cost = sum([0; [get(fact_costs, f, 0) for f in body]])
        else
            cost = maximum([0; [get(fact_costs, f, 0) for f in body]])
        end
        derived = substitute(ax.head, s)
        if cost < get(fact_costs, derived, Inf)
            fact_costs[derived] = cost end
    end
    return fact_costs
end

"Compute costs of all effects of available actions"
function compute_cost_action_effect(fact_costs::Dict{Any,Any}, act::Term, domain::Domain, additions::Dict{Any, Any}, heur::String)
    act_args = domain.actions[act.name].args
    subst = Subst(var => val for (var, val) in zip(act_args, act.args))
    # Look-up preconds and substitute vars
    conds = g_preconditions[act.name]
    conds = [[substitute(t, subst) for t in c] for c in conds]
    # Compute cost of reaching each action
    if heur == "delete_relaxation/h_add"
        cost = minimum([[sum([0; [get(fact_costs, f, 0) for f in conj]]) for conj in conds]; Inf])
    else
        cost = minimum([[maximum([0; [get(fact_costs, f, 0) for f in conj]]) for conj in conds]; Inf])
    end
    # Compute cost of reaching each added fact
    added = [substitute(a, subst) for a in additions[string(act.name)]]
    cost = cost + get_arbitrary_cost(act, domain) # returns 1 if arbitrary cost isn't found (only total-cost key is checked)
    for fact in added
        if cost < get(fact_costs, fact, Inf)
            fact_costs[fact] = cost end
    end
    return fact_costs
end

"Get arbitrary cost of action term"
function get_arbitrary_cost(act::Term, domain::Domain)
    effects = domain.actions[act.name].effect
    for e in effects.args
        if length(e.args) == 0
            continue
        end
        if string(e.args[1]) == "total-cost"
            if string(e.name) == "increase"
                return parse(Float64, string(e.args[2])) * 1
            elseif string(e.name) == "decrease"
                return parse(Float64, string(e.args[2])) * -1
            end
        end
    end
    return 1
end

"Return len of Julia Set"
function length_set(set::Set)
    return length(set)
end

"Initialize facts costs"
function convert_any_to_terms(set_any::Set{Any})
    return [parse(Term, s) for s in set_any]
end

"Initialize fact costs"
function init_facts_costs(facts::Set{Term})
    return Dict{Term,Float64}(f => 0 for f in facts)
end

"Get facts and state"
function get_facts_and_state(facts_costs::Dict{Any,Any}, types::Set{Term})
    real_facts_costs = convert_dict_any_to_term(facts_costs)
    facts = Set(keys(real_facts_costs))
    state = State(types, facts, Dict{Symbol,Any}())
    return facts, state
end

"Convert any any to term float64 dicts"
function convert_dict_any_to_term(dict::Dict{Any, Any})
    return Dict{Term,Float64}(f => c for (f, c) in dict)
end

"Convert any any to symbol vector(vector(term)) dicts"
function convert_dict_any_to_preconds(dict::Dict{Any, Any})
    return Dict{Symbol, Vector{Vector{Term}}}(Symbol(f) => c for (f, c) in dict)
end


text:
using Test
using ColorTypes: RGB, Gray, N0f8, red, green, blue
import ColorVectorSpace
using FileIO, ImageCore, Dates, Statistics, StatsBase

import FFMPEG

import VideoIO

const testdir = dirname(@__FILE__)
const videodir = joinpath(testdir, "..", "videos")
const tempvidname = "testvideo.mp4"
const tempvidpath = joinpath(tempdir(), tempvidname)
const required_accuracy = 0.07

# VideoIO.TestVideos.available()
VideoIO.TestVideos.download_all()

include("utils.jl") # Testing utility functions

include("avptr.jl")
include("reading.jl")
include("writing.jl")
include("accuracy.jl")

GC.gc()
rm(tempvidpath, force = true)

include("bugs.jl")

#VideoIO.TestVideos.remove_all()


text:
using DatasToolbox
using DataFrames
# using ScikitLearn
# fit! = ScikitLearn.fit!
# predict = ScikitLearn.predict
# @sk_import ensemble: RandomForestRegressor
using XGBoost
predict = XGBoost.predict

df_train = DataFrame(ξ=1:2:1000, ζ=2:2:1000, τ=1:500)
df_test = DataFrame(ξ=502:2:520, ζ=501:2:520, τ=501:510)
df = [df_train; df_test]

tsh = TimeSeriesHandler{Float64}(df, :τ, 5, input_cols=[:ξ], output_cols=[:ξ],
                                 n_test_sequences=0,
                                 userange=true)
split!(tsh, 505)
# calling this again to be safe
assign!(tsh)

X_train, y_train = getSquashedTrainData(tsh, flatten=true)

booster = xgboost(X_train, 500, label=y_train, objective="reg:linear")

# yhat = generateTest(X -> predict(booster, X), tsh, on_matrix=true)
yhat = generateSequence(X -> predict(booster, X), tsh, 10, on_matrix=true)

max_stars_repo_name: JuliaTagBot/DatasToolbox

max_issues_repo_name: JuliaTagBot/DatasToolbox

max_forks_repo_name: JuliaTagBot/DatasToolbox