text:
# ============================================================
# Solution Algorithm
# ============================================================

"""
Calculate time step

"""
function timestep(KS::SolverSet, uq::AbstractUQ, sol::AbstractSolution, simTime::Real)

    tmax = 0.0

    if KS.set.nSpecies == 1

        if KS.set.space[1:2] == "1d"
            @inbounds Threads.@threads for i = 1:KS.pSpace.nx
                sos = uq_sound_speed(sol.prim[i], KS.gas.γ, uq)
                vmax = KS.vSpace.u1 + maximum(sos)
                tmax = max(tmax, vmax / KS.pSpace.dx[i])
            end
        elseif KS.set.space[1:2] == "2d"
            @inbounds Threads.@threads for j = 1:KS.pSpace.ny
                for i = 1:KS.pSpace.nx
                    sos = uq_sound_speed(sol.prim[i, j], KS.gas.γ, uq)
                    vmax = max(KS.vSpace.u1, KS.vSpace.v1) + maximum(sos)
                    tmax = max(tmax, vmax / KS.pSpace.dx[i, j] + vmax / KS.pSpace.dy[i, j])
                end
            end
        end

    elseif KS.set.nSpecies == 2

        @inbounds Threads.@threads for i = 1:KS.pSpace.nx
            prim = sol.prim[i]
            sos = uq_sound_speed(prim, KS.gas.γ, uq)
            vmax = max(maximum(KS.vSpace.u1), maximum(abs.(prim[2, :, :]))) + maximum(sos)
            tmax = max(tmax, vmax / KS.pSpace.dx[i])
        end

    end

    dt = KS.set.cfl / tmax
    dt = ifelse(dt < (KS.set.maxTime - simTime), dt, KS.set.maxTime - simTime)

    return dt

end

function timestep(
    KS::T1,
    ctr::T2,
    simTime,
    uq::T3,
) where {
    T1<:AbstractSolverSet,
    T2<:Union{AbstractArray{ControlVolume1D1F,1},AbstractArray{ControlVolume1D2F,1}},
    T3<:AbstractUQ,
}

    tmax = 0.0
    Threads.@threads for i = 1:KS.pSpace.nx
        @inbounds prim = ctr[i].prim
        sos = uq_sound_speed(prim, KS.gas.γ, uq)
        vmax = max(maximum(KS.vSpace.u1), maximum(abs.(prim[2, :]))) + maximum(sos)
        tmax = max(tmax, vmax / ctr[i].dx)
    end

    dt = KS.set.cfl / tmax
    dt = ifelse(dt < (KS.set.maxTime - simTime), dt, KS.set.maxTime - simTime)

    return dt

end

function timestep(
    KS::SolverSet,
    ctr::AbstractArray{ControlVolume1D4F,1},
    simTime::Real,
    uq::AbstractUQ,
)
    tmax = 0.0

    Threads.@threads for i = 1:KS.pSpace.nx
        @inbounds prim = ctr[i].prim
        sos = uq_sound_speed(prim, KS.gas.γ, uq)
        vmax = max(maximum(KS.vSpace.u1), maximum(abs.(prim[2, :, :]))) + maximum(sos)
        tmax = max(tmax, vmax / ctr[i].dx, KS.gas.sol / ctr[i].dx)
    end

    dt = KS.set.cfl / tmax
    dt = ifelse(dt < (KS.set.maxTime - simTime), dt, KS.set.maxTime - simTime)

    return dt
end

function timestep(
    KS::SolverSet,
    ctr::AbstractArray{ControlVolume1D3F,1},
    simTime::Real,
    uq::AbstractUQ,
)
    tmax = 0.0

    @inbounds Threads.@threads for i = 1:KS.pSpace.nx
        prim = ctr[i].prim
        sos = uq_sound_speed(prim, KS.gas.γ, uq)
        vmax = max(maximum(KS.vSpace.u1), maximum(abs.(prim[2, :, :]))) + maximum(sos)
        smax = maximum(abs.(ctr[i].lorenz))
        tmax = max(tmax, vmax / ctr[i].dx, KS.gas.sol / ctr[i].dx, smax / KS.vSpace.du[1])
    end

    dt = KS.set.cfl / tmax
    dt = ifelse(dt < (KS.set.maxTime - simTime), dt, KS.set.maxTime - simTime)

    return dt
end

"""
Reconstruct solution

"""
function reconstruct!(KS::SolverSet, sol::Solution1D1F)

    @inbounds Threads.@threads for i = 1:KS.pSpace.nx
        KitBase.reconstruct3!(
            sol.sw[i],
            sol.w[i-1],
            sol.w[i],
            sol.w[i+1],
            0.5 * (KS.pSpace.dx[i-1] + KS.pSpace.dx[i]),
            0.5 * (KS.pSpace.dx[i] + KS.pSpace.dx[i+1]),
        )

        KitBase.reconstruct3!(
            sol.sf[i],
            sol.f[i-1],
            sol.f[i],
            sol.f[i+1],
            0.5 * (KS.pSpace.dx[i-1] + KS.pSpace.dx[i]),
            0.5 * (KS.pSpace.dx[i] + KS.pSpace.dx[i+1]),
        )
    end

end

function reconstruct!(KS::SolverSet, sol::Solution2D2F)

    #--- x direction ---#
    @inbounds Threads.@threads for j = 1:KS.pSpace.ny
        KitBase.reconstruct2!(
            sol.sw[1, j][:, :, 1],
            sol.w[1, j],
            sol.w[2, j],
            0.5 * (KS.pSpace.dx[1, j] + KS.pSpace.dx[2, j]),
        )
        KitBase.reconstruct2!(
            sol.sh[1, j][:, :, :, 1],
            sol.h[1, j],
            sol.h[2, j],
            0.5 * (KS.pSpace.dx[1, j] + KS.pSpace.dx[2, j]),
        )
        KitBase.reconstruct2!(
            sol.sb[1, j][:, :, :, 1],
            sol.b[1, j],
            sol.b[2, j],
            0.5 * (KS.pSpace.dx[1, j] + KS.pSpace.dx[2, j]),
        )

        KitBase.reconstruct2!(
            sol.sw[KS.pSpace.nx, j][:, :, 1],
            sol.w[KS.pSpace.nx-1, j],
            sol.w[KS.pSpace.nx, j],
            0.5 * (KS.pSpace.dx[KS.pSpace.nx-1, j] + KS.pSpace.dx[KS.pSpace.nx, j]),
        )
        KitBase.reconstruct2!(
            sol.sh[KS.pSpace.nx, j][:, :, :, 1],
            sol.h[KS.pSpace.nx-1, j],
            sol.h[KS.pSpace.nx, j],
            0.5 * (KS.pSpace.dx[KS.pSpace.nx-1, j] + KS.pSpace.dx[KS.pSpace.nx, j]),
        )
        KitBase.reconstruct2!(
            sol.sb[KS.pSpace.nx, j][:, :, :, 1],
            sol.b[KS.pSpace.nx-1, j],
            sol.b[KS.pSpace.nx, j],
            0.5 * (KS.pSpace.dx[KS.pSpace.nx-1, j] + KS.pSpace.dx[KS.pSpace.nx, j]),
        )
    end

    @inbounds Threads.@threads for j = 1:KS.pSpace.ny
        for i = 2:KS.pSpace.nx-1
            KitBase.reconstruct3!(
                sol.sw[i, j][:, :, 1],
                sol.w[i-1, j],
                sol.w[i, j],
                sol.w[i+1, j],
                0.5 * (KS.pSpace.dx[i-1, j] + KS.pSpace.dx[i, j]),
                0.5 * (KS.pSpace.dx[i, j] + KS.pSpace.dx[i+1, j]),
            )

            KitBase.reconstruct3!(
                sol.sh[i, j][:, :, :, 1],
                sol.h[i-1, j],
                sol.h[i, j],
                sol.h[i+1, j],
                0.5 * (KS.pSpace.dx[i-1, j] + KS.pSpace.dx[i, j]),
                0.5 * (KS.pSpace.dx[i, j] + KS.pSpace.dx[i+1, j]),
            )
            KitBase.reconstruct3!(
                sol.sb[i, j][:, :, :, 1],
                sol.b[i-1, j],
                sol.b[i, j],
                sol.b[i+1, j],
                0.5 * (KS.pSpace.dx[i-1, j] + KS.pSpace.dx[i, j]),
                0.5 * (KS.pSpace.dx[i, j] + KS.pSpace.dx[i+1, j]),
            )
        end
    end

    #--- y direction ---#
    @inbounds Threads.@threads for i = 1:KS.pSpace.nx
        KitBase.reconstruct2!(
            sol.sw[i, 1][:, :, 2],
            sol.w[i, 1],
            sol.w[i, 2],
            0.5 * (KS.pSpace.dy[i, 1] + KS.pSpace.dy[i, 2]),
        )
        KitBase.reconstruct2!(
            sol.sh[i, 1][:, :, :, 2],
            sol.h[i, 1],
            sol.h[i, 2],
            0.5 * (KS.pSpace.dy[i, 1] + KS.pSpace.dy[i, 2]),
        )
        KitBase.reconstruct2!(
            sol.sb[i, 1][:, :, :, 2],
            sol.b[i, 1],
            sol.b[i, 2],
            0.5 * (KS.pSpace.dy[i, 1] + KS.pSpace.dy[i, 2]),
        )

        KitBase.reconstruct2!(
            sol.sw[i, KS.pSpace.ny][:, :, 2],
            sol.w[i, KS.pSpace.ny-1],
            sol.w[i, KS.pSpace.ny],
            0.5 * (KS.pSpace.dy[i, KS.pSpace.ny-1] + KS.pSpace.dy[i, KS.pSpace.ny]),
        )
        KitBase.reconstruct2!(
            sol.sh[i, KS.pSpace.ny][:, :, :, 2],
            sol.h[i, KS.pSpace.ny-1],
            sol.h[i, KS.pSpace.ny],
            0.5 * (KS.pSpace.dy[i, KS.pSpace.ny-1] + KS.pSpace.dy[i, KS.pSpace.ny]),
        )
        KitBase.reconstruct2!(
            sol.sb[i, KS.pSpace.ny][:, :, :, 2],
            sol.b[i, KS.pSpace.ny-1],
            sol.b[i, KS.pSpace.ny],
            0.5 * (KS.pSpace.dy[i, KS.pSpace.ny-1] + KS.pSpace.dy[i, KS.pSpace.ny]),
        )
    end

    @inbounds Threads.@threads for j = 2:KS.pSpace.ny-1
        for i = 1:KS.pSpace.nx
            KitBase.reconstruct3!(
                sol.sw[i, j][:, :, 2],
                sol.w[i, j-1],
                sol.w[i, j],
                sol.w[i, j+1],
                0.5 * (KS.pSpace.dy[i, j-1] + KS.pSpace.dy[i, j]),
                0.5 * (KS.pSpace.dy[i, j] + KS.pSpace.dy[i, j+1]),
            )

            KitBase.reconstruct3!(
                sol.sh[i, j][:, :, :, 2],
                sol.h[i, j-1],
                sol.h[i, j],
                sol.h[i, j+1],
                0.5 * (KS.pSpace.dy[i, j-1] + KS.pSpace.dy[i, j]),
                0.5 * (KS.pSpace.dy[i, j] + KS.pSpace.dy[i, j+1]),
            )
            KitBase.reconstruct3!(
                sol.sb[i, j][:, :, :, 2],
                sol.b[i, j-1],
                sol.b[i, j],
                sol.b[i, j+1],
                0.5 * (KS.pSpace.dy[i, j-1] + KS.pSpace.dy[i, j]),
                0.5 * (KS.pSpace.dy[i, j] + KS.pSpace.dy[i, j+1]),
            )
        end
    end

end

"""
Update solution

"""
function update!(
    KS::SolverSet,
    uq::AbstractUQ,
    sol::Solution1D1F,
    flux::Flux1D1F,
    dt::Float64,
    residual::Array{Float64,1},
)

    w_old = deepcopy(sol.w)
    #=
        @. sol.w[2:KS.pSpace.nx-1] +=
            (flux.fw[2:end - 2] - flux.fw[3:end-1]) / KS.pSpace.dx[2:KS.pSpace.nx-1]
        uq_conserve_prim!(sol, KS.gas.γ, uq)
    =#
    @inbounds Threads.@threads for i = 1:KS.pSpace.nx
        @. sol.w[i] += (flux.fw[i] - flux.fw[i+1]) / KS.pSpace.dx[i]
        sol.prim[i] .= uq_conserve_prim(sol.w[i], KS.gas.γ, uq)
    end

    τ = uq_vhs_collision_time(sol, KS.gas.μᵣ, KS.gas.ω, uq)
    M = [uq_maxwellian(KS.vSpace.u, sol.prim[i], uq) for i in eachindex(sol.prim)]

    @inbounds Threads.@threads for i = 1:KS.pSpace.nx
        for j in axes(sol.w[1], 2)
            @. sol.f[i][:, j] =
                (
                    sol.f[i][:, j] +
                    (flux.ff[i][:, j] - flux.ff[i+1][:, j]) / KS.pSpace.dx[i] +
                    dt / τ[i][j] * M[i][:, j]
                ) / (1.0 + dt / τ[i][j])
        end
    end

    # record residuals
    sumRes = zeros(axes(KS.ib.wL, 1))
    sumAvg = zeros(axes(KS.ib.wL, 1))
    for j in axes(sumRes, 1)
        for i = 1:KS.pSpace.nx
            sumRes[j] += sum((sol.w[i][j, :] .- w_old[i][j, :]) .^ 2)
            sumAvg[j] += sum(abs.(sol.w[i][j, :]))
        end
    end
    @. residual = sqrt(sumRes * KS.pSpace.nx) / (sumAvg + 1.e-7)

end

function update!(
    KS::SolverSet,
    uq::AbstractUQ,
    sol::Solution2D2F,
    flux::Flux2D2F,
    dt::Float64,
    residual::Array{Float64,1},
)

    w_old = deepcopy(sol.w)

    @inbounds Threads.@threads for j = 1:KS.pSpace.ny
        for i = 1:KS.pSpace.nx
            @. sol.w[i, j] +=
                (flux.fw1[i, j] - flux.fw1[i+1, j] + flux.fw2[i, j] - flux.fw2[i, j+1]) /
                (KS.pSpace.dx[i, j] * KS.pSpace.dy[i, j])
            sol.prim[i, j] .= uq_conserve_prim(sol.w[i, j], KS.gas.γ, uq)
        end
    end

    τ = uq_vhs_collision_time(sol, KS.gas.μᵣ, KS.gas.ω, uq)
    H = [
        uq_maxwellian(KS.vSpace.u, KS.vSpace.v, sol.prim[i, j], uq) for
        i in axes(sol.prim, 1), j in axes(sol.prim, 2)
    ]
    B = deepcopy(H)
    for i in axes(sol.prim, 1), j in axes(sol.prim, 2)
        for k in axes(B[1, 1], 3)
            B[i, j][:, :, k] .=
                B[i, j][:, :, k] .* KS.gas.K ./ (2.0 * sol.prim[i, j][end, k])
        end
    end

    @inbounds Threads.@threads for i = 1:KS.pSpace.nx
        for j = 1:KS.pSpace.ny
            for k in axes(sol.w[1, 1], 2)
                @. sol.h[i, j][:, :, k] =
                    (
                        sol.h[i, j][:, :, k] +
                        (
                            flux.fh1[i, j][:, :, k] - flux.fh1[i+1, j][:, :, k] +
                            flux.fh2[i, j][:, :, k] - flux.fh2[i, j+1][:, :, k]
                        ) / (KS.pSpace.dx[i, j] * KS.pSpace.dy[i, j]) +
                        dt / τ[i, j][k] * H[i, j][:, :, k]
                    ) / (1.0 + dt / τ[i, j][k])
                @. sol.b[i, j][:, :, k] =
                    (
                        sol.b[i, j][:, :, k] +
                        (
                            flux.fb1[i, j][:, :, k] - flux.fb1[i+1, j][:, :, k] +
                            flux.fb2[i, j][:, :, k] - flux.fb2[i, j+1][:, :, k]
                        ) / (KS.pSpace.dx[i, j] * KS.pSpace.dy[i, j]) +
                        dt / τ[i, j][k] * B[i, j][:, :, k]
                    ) / (1.0 + dt / τ[i, j][k])
            end
        end
    end

    # record residuals
    sumRes = zeros(axes(KS.ib.wL, 1))
    sumAvg = zeros(axes(KS.ib.wL, 1))
    @inbounds for k in axes(sumRes, 1)
        for j = 1:KS.pSpace.ny
            for i = 1:KS.pSpace.nx
                sumRes[k] += sum((sol.w[i, j][k, :] .- w_old[i, j][k, :]) .^ 2)
                sumAvg[k] += sum(abs.(sol.w[i, j][k, :]))
            end
        end
    end
    @. residual = sqrt(sumRes * KS.pSpace.nx * KS.pSpace.ny) / (sumAvg + 1.e-7)

end

function step!(
    KS::SolverSet,
    uq::AbstractUQ,
    sol::Solution2D2F,
    flux::Flux2D2F,
    dt::Float64,
    residual::Array{Float64,1},
)

    sumRes = zeros(axes(KS.ib.wL, 1))
    sumAvg = zeros(axes(KS.ib.wL, 1))

    Threads.@threads for j = 1:KS.pSpace.ny
        for i = 1:KS.pSpace.nx
            step!(
                KS,
                uq,
                sol.w[i, j],
                sol.prim[i, j],
                sol.h[i, j],
                sol.b[i, j],
                flux.fw1[i, j],
                flux.fh1[i, j],
                flux.fb1[i, j],
                flux.fw1[i+1, j],
                flux.fh1[i+1, j],
                flux.fb1[i+1, j],
                flux.fw2[i, j],
                flux.fh2[i, j],
                flux.fb2[i, j],
                flux.fw2[i, j+1],
                flux.fh2[i, j+1],
                flux.fb2[i, j+1],
                dt,
                KS.pSpace.dx[i, j] * KS.pSpace.dy[i, j],
                sumRes,
                sumAvg,
            )
        end
    end

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

end

function step!(
    KS::SolverSet,
    uq::AbstractUQ,
    w::Array{Float64,2},
    prim::Array{Float64,2},
    h::AbstractArray{Float64,3},
    b::AbstractArray{Float64,3},
    fwL::Array{Float64,2},
    fhL::AbstractArray{Float64,3},
    fbL::AbstractArray{Float64,3},
    fwR::Array{Float64,2},
    fhR::AbstractArray{Float64,3},
    fbR::AbstractArray{Float64,3},
    fwU::Array{Float64,2},
    fhU::AbstractArray{Float64,3},
    fbU::AbstractArray{Float64,3},
    fwD::Array{Float64,2},
    fhD::AbstractArray{Float64,3},
    fbD::AbstractArray{Float64,3},
    dt::Float64,
    area::Float64,
    sumRes::Array{Float64,1},
    sumAvg::Array{Float64,1},
)

    w_old = deepcopy(w)

    @. w += (fwL - fwR + fwD - fwU) / area
    prim .= uq_conserve_prim(w, KS.gas.γ, uq)

    τ = uq_vhs_collision_time(prim, KS.gas.μᵣ, KS.gas.ω, uq)
    H = uq_maxwellian(KS.vSpace.u, KS.vSpace.v, prim, uq)
    B = similar(H)
    for k in axes(H, 3)
        B[:, :, k] .= H[:, :, k] .* KS.gas.K ./ (2.0 * prim[end, k])
    end

    for k in axes(h, 3)
        @. h[:, :, k] =
            (
                h[:, :, k] +
                (fhL[:, :, k] - fhR[:, :, k] + fhD[:, :, k] - fhU[:, :, k]) / area +
                dt / τ[k] * H[:, :, k]
            ) / (1.0 + dt / τ[k])
        @. b[:, :, k] =
            (
                b[:, :, k] +
                (fbL[:, :, k] - fbR[:, :, k] + fbD[:, :, k] - fbU[:, :, k]) / area +
                dt / τ[k] * B[:, :, k]
            ) / (1.0 + dt / τ[k])
    end

    for i = 1:4
        sumRes[i] += sum((w[i, :] .- w_old[i, :]) .^ 2)
        sumAvg[i] += sum(abs.(w[i, :]))
    end

end

function update!(
    KS::SolverSet,
    uq::AbstractUQ,
    ctr::AbstractArray{ControlVolume1D1F,1},
    face::AbstractArray{Interface1D1F,1},
    dt,
    residual::Array{<:Real,1};
    coll = :bgk::Symbol,
    bc = :fix::Symbol,
) where {T1<:AbstractSolverSet}

    sumRes = zeros(3)
    sumAvg = zeros(3)

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

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

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

end

function update!(
    KS::SolverSet,
    uq::AbstractUQ,
    ctr::AbstractArray{ControlVolume1D3F,1},
    face::AbstractArray{Interface1D3F,1},
    dt::Real,
    residual::Array{<:Real,2};
    coll = :bgk::Symbol,
    bc = :extra::Symbol,
    isMHD = true::Bool,
)

    sumRes = zeros(5, 2)
    sumAvg = zeros(5, 2)

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

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

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

end

function update!(
    KS::SolverSet,
    uq::AbstractUQ,
    ctr::AbstractArray{ControlVolume1D4F,1},
    face::AbstractArray{Interface1D4F,1},
    dt::Real,
    residual::Array{<:Real,2};
    coll = :bgk::Symbol,
    bc = :extra::Symbol,
    isMHD = true::Bool,
)

    sumRes = zeros(5, 2)
    sumAvg = zeros(5, 2)

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

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

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

end

function update_boundary!(
    KS::SolverSet,
    uq::AbstractUQ,
    ctr::AbstractArray{<:AbstractControlVolume1D,1},
    face::Array{<:AbstractInterface1D,1},
    dt::Real,
    residual::Array{<:AbstractFloat};
    coll = :bgk::Symbol,
    bc = :extra::Symbol,
    isMHD = true::Bool,
)

    if bc != :fix
        resL = zeros(size(KS.ib.wL, 1), KS.set.nSpecies)
        avgL = zeros(size(KS.ib.wL, 1), KS.set.nSpecies)
        resR = zeros(size(KS.ib.wL, 1), KS.set.nSpecies)
        avgR = zeros(size(KS.ib.wL, 1), KS.set.nSpecies)

        i = 1
        j = KS.pSpace.nx
        if KS.set.space[3:4] == "0f"
            step!(KS, uq, face[i], ctr[i], face[i+1], dt, resL, avgL)
            step!(KS, uq, face[j], ctr[j], face[j+1], dt, resR, avgR)
        elseif KS.set.space[3:4] in ["3f", "4f"]
            step!(KS, uq, face[i], ctr[i], face[i+1], dt, resL, avgL, coll, isMHD)
            step!(KS, uq, face[j], ctr[j], face[j+1], dt, resR, avgR, coll, isMHD)
        else
            step!(KS, uq, face[i], ctr[i], face[i+1], dt, resL, avgL, coll)
            step!(KS, uq, face[j], ctr[j], face[j+1], dt, resR, avgR, coll)
        end

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

    ng = 1 - first(eachindex(KS.pSpace.x))
    if bc == :extra

        for i = 1:ng
            ctr[1-i].w .= ctr[1].w
            ctr[1-i].prim .= ctr[1].prim
            ctr[KS.pSpace.nx+i].w .= ctr[KS.pSpace.nx].w
            ctr[KS.pSpace.nx+i].prim .= ctr[KS.pSpace.nx].prim

            if KS.set.space[3:4] == "1f"
                ctr[1-i].f .= ctr[1].f
                ctr[KS.pSpace.nx+i].f .= ctr[KS.pSpace.nx].f
            elseif KS.set.space[3:4] == "2f"
                ctr[1-i].h .= ctr[1].h
                ctr[1-i].b .= ctr[1].b
                ctr[KS.pSpace.nx+i].h .= ctr[KS.pSpace.nx].h
                ctr[KS.pSpace.nx+i].b .= ctr[KS.pSpace.nx].b
            elseif KS.set.space[3:4] == "3f"
                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].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
            elseif KS.set.space[3:4] == "4f"
                ctr[1-i].h0 .= ctr[1].h0
                ctr[1-i].h1 .= ctr[1].h1
                ctr[1-i].h2 .= ctr[1].h2
                ctr[1-i].h3 .= ctr[1].h3
                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].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].h3 .= ctr[KS.pSpace.nx].h3
                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
            else
                throw("incorrect amount of distribution functions")
            end
        end

    elseif bc == :period

        for i = 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[KS.pSpace.nx+i].w .= ctr[i].w
            ctr[KS.pSpace.nx+i].prim .= ctr[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].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

            if KS.set.space[3:4] == "1f"
                ctr[1-i].f .= ctr[KS.pSpace.nx+1-i].f
                ctr[KS.pSpace.nx+i].f .= ctr[i].f
            elseif KS.set.space[3:4] == "2f"
                ctr[1-i].h .= ctr[KS.pSpace.nx+1-i].h
                ctr[1-i].b .= ctr[KS.pSpace.nx+1-i].b
                ctr[KS.pSpace.nx+i].h .= ctr[i].h
                ctr[KS.pSpace.nx+i].b .= ctr[i].b
            elseif KS.set.space[3:4] == "3f"
                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].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
            elseif KS.set.space[3:4] == "4f"
                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].h3 .= ctr[KS.pSpace.nx+1-i].h3
                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].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].h3 .= ctr[i].h3
                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
            else
                throw("incorrect amount of distribution functions")
            end
        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)

        if KS.set.space[3:4] == "1f"
            @. ctr[0].f = 0.5 * (ctr[-1].f + ctr[1].f)
            @. ctr[KS.pSpace.nx+1].f = 0.5 * (ctr[KS.pSpace.nx].f + ctr[KS.pSpace.nx+2].f)
        elseif KS.set.space[3:4] == "2f"
            @. ctr[0].h = 0.5 * (ctr[-1].h + ctr[1].h)
            @. ctr[0].b = 0.5 * (ctr[-1].b + ctr[1].b)
            @. ctr[KS.pSpace.nx+1].h = 0.5 * (ctr[KS.pSpace.nx].h + ctr[KS.pSpace.nx+2].h)
            @. ctr[KS.pSpace.nx+1].b = 0.5 * (ctr[KS.pSpace.nx].b + ctr[KS.pSpace.nx+2].b)
        elseif KS.set.space[3:4] == "3f"
            @. 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].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)
        elseif KS.set.space[3:4] == "4f"
            @. 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].h3 = 0.5 * (ctr[-1].h3 + ctr[1].h3)
            @. 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].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].h3 =
                0.5 * (ctr[KS.pSpace.nx].h3 + ctr[KS.pSpace.nx+2].h3)
            @. 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
            throw("incorrect amount of distribution functions")
        end

    else

    end

end


text:
"""
    struct InfiniteMPS{A<:GenericMPSTensor,B<:MPSBondTensor}

Represents an infinite matrix product state
The state is stored in the centergauge where
    state.AL[i]*state.CR[i] = state.AC[i] = state.CR[i-1]*state.AR[i]
"""
struct InfiniteMPS{A<:GenericMPSTensor,B<:MPSBondTensor}
    AL::PeriodicArray{A,1}
    AR::PeriodicArray{A,1}
    CR::PeriodicArray{B,1}
    AC::PeriodicArray{A,1}
end

Base.size(arr::InfiniteMPS,i) = size(arr.AL,i)
Base.length(arr::InfiniteMPS) = size(arr,1)
Base.eltype(arr::InfiniteMPS) = eltype(arr.AL)
Base.copy(m::InfiniteMPS) = InfiniteMPS(copy(m.AL),copy(m.AR),copy(m.CR),copy(m.AC));
Base.repeat(m::InfiniteMPS,i::Int) = InfiniteMPS(repeat(m.AL,i),repeat(m.AR,i),repeat(m.CR,i),repeat(m.AC,i));
Base.similar(st::InfiniteMPS) = InfiniteMPS(similar(st.AL),similar(st.AR),similar(st.CR),similar(st.AC))
TensorKit.norm(st::InfiniteMPS) = norm(st.AC[1]);
virtualspace(psi::InfiniteMPS, n::Integer) = _firstspace(psi.AL[n+1])

function InfiniteMPS(pspaces::AbstractArray{S,1},Dspaces::AbstractArray{S,1};kwargs...) where S
    InfiniteMPS([TensorMap(rand,Defaults.eltype,Dspaces[mod1(i-1,length(Dspaces))]*pspaces[i],Dspaces[i]) for i in 1:length(pspaces)];kwargs...)
end

#allow users to pass in simple arrays
function InfiniteMPS(A::AbstractArray{T,1}; kwargs...) where T<:GenericMPSTensor

    #we make a copy, and are therfore garantueeing no side effects for the user
    AR = PeriodicArray(A[:]);

    #initial guess for CR
    CR = PeriodicArray([isomorphism(Matrix{eltype(AR[i])},space(AR[i+1],1),space(AR[i+1],1)) for i in 1:length(A)]);

    AL = similar(AR);

    uniform_leftorth!(AL,CR,AR;kwargs...);
    uniform_rightorth!(AR,CR,AL;kwargs...);

    AC = similar(AR)
    for loc = 1:length(A)
        AC[loc] = AL[loc]*CR[loc]
    end

    CRtype = tensormaptype(spacetype(T),1,1,eltype(T));
    return InfiniteMPS{T,CRtype}(AL,AR,CR,AC)
end

function InfiniteMPS(A::AbstractArray{T,1},cinit::B;kwargs...) where {T<:GenericMPSTensor, B<:MPSBondTensor}
    CR = PeriodicArray(fill(copy(cinit),length(A)));
    AL = PeriodicArray(A[:])
    AR = similar(AL);
    uniform_rightorth!(AR,CR,AL;kwargs...);

    AC = similar(AR)
    for loc = 1:length(A)
        AC[loc] = AL[loc]*CR[loc]
    end
    return InfiniteMPS{T,B}(AL,AR,CR,AC)
end

function TensorKit.normalize!(st::InfiniteMPS)
    normalize!.(st.CR)
    normalize!.(st.AC)
    st
end

"
    l_RR(state,location)
    Left dominant eigenvector of the AR-AR transfermatrix
"
l_RR(state::InfiniteMPS,loc::Int=1) = @tensor toret[-1;-2]:=state.CR[loc-1][1,-2]*conj(state.CR[loc-1][1,-1])

"
    l_RL(state,location)
    Left dominant eigenvector of the AR-AL transfermatrix
"
l_RL(state::InfiniteMPS,loc::Int=1) = state.CR[loc-1]

"
    l_LR(state,location)
    Left dominant eigenvector of the AL-AR transfermatrix
"
l_LR(state::InfiniteMPS,loc::Int=1) = state.CR[loc-1]'

"
    l_LL(state,location)
    Left dominant eigenvector of the AL-AL transfermatrix
"
l_LL(state::InfiniteMPS{A},loc::Int=1) where A= isomorphism(Matrix{eltype(A)}, space(state.AL[loc],1),space(state.AL[loc],1))

"
    r_RR(state,location)
    Right dominant eigenvector of the AR-AR transfermatrix
"
r_RR(state::InfiniteMPS{A},loc::Int=length(state)) where A= isomorphism(Matrix{eltype(A)},domain(state.AR[loc]),domain(state.AR[loc]))

"
    r_RL(state,location)
    Right dominant eigenvector of the AR-AL transfermatrix
"
r_RL(state::InfiniteMPS,loc::Int=length(state)) = state.CR[loc]'

"
    r_LR(state,location)
    Right dominant eigenvector of the AL-AR transfermatrix
"
r_LR(state::InfiniteMPS,loc::Int=length(state)) = state.CR[loc]

"
    r_LL(state,location)
    Right dominant eigenvector of the AL-AL transfermatrix
"
r_LL(state::InfiniteMPS,loc::Int=length(state))= @tensor toret[-1;-2]:=state.CR[loc][-1,1]*conj(state.CR[loc][-2,1])

function TensorKit.dot(a::InfiniteMPS,b::InfiniteMPS;krylovdim = 30)
    init = TensorMap(rand,ComplexF64,space(a.AL[1],1),space(b.AL[1],1))
    num = lcm(length(a),length(b))
    (vals,vecs,convhist) = eigsolve(x->transfer_left(x,b.AL[1:num],a.AL[1:num]),init,1,:LM,Arnoldi(krylovdim=krylovdim))
    convhist.converged == 0 && @info "dot mps not converged"
    return vals[1]
end

Base.circshift(st::InfiniteMPS,n) = InfiniteMPS(circshift(st.AL,n),circshift(st.AR,n),circshift(st.CR,n),circshift(st.AC,n))

ext: jl


text:
module SlurmSweeps
using Logging
using DelimitedFiles
using Random
using Sockets

include("server.jl")
include("client.jl")
include("instance.jl")
include("setup.jl")

end # module


text:
#######################################################################
function test_run_ppca(σ, seed)
#######################################################################

    x, = simulate3Ddata(σ, seed)

    W, μ, σ, proj, fwd = ppca(x')

    figure(1)

    cla()

    plot3D(x[:,1], x[:,2], x[:,3], "bo", alpha=0.3)

    z = sort([proj(x[i,:]) for i=1:size(x,1)])

    zgrid = collect(LinRange(minimum(z), maximum(z), 100))

    X = reduce(hcat, fwd.(zgrid))'

    plot3D(X[:,1], X[:,2], X[:,3], "-")

    return W, μ

end


text:
#!/usr/bin/env juila
"""
Run bgzip on a file. The script will either decompress or
compress a file depending on whether infile ends in ".gz".
Returns the name of the outfile.

bgzip_file(infile)

"""
function bgzip_file(infile::String)
    if occursin(r".gz$", infile)
        run(`bash -c "bgzip -d -@8 $infile"`)
        filename = splitext(infile)
        outfile = filename[1]
    elseif ~occursin(r".gz$", infile)
        run(`bash -c "bgzip -@8 $infile"`)
        filename = join([infile, "gz"], ".")
        outfile = filename
    else
        println("ERROR: Could not bgzip or unzip file $filein.\n\n")
        exit()
    end
    return outfile
end

max_stars_repo_name: ScottWatkins/STRdiff

max_issues_repo_name: ScottWatkins/STRdiff

max_forks_repo_name: ScottWatkins/STRdiff


text:
"""
    penalty objective
"""
struct PenaltyObjective{T} <: Objective
    α::T
    idx::Int
end

function objective(Z, obj::PenaltyObjective, idx, T)
    J = 0.0
    for t = 1:T-1
        s = Z[idx.u[t][obj.idx]]
        J += s
    end
    return obj.α * J
end

function objective_gradient!(∇J, Z, obj::PenaltyObjective, idx, T)
    for t = 1:T-1
        ∇J[idx.u[t][obj.idx]] += obj.α
    end
    return nothing
end

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
using DelimitedFiles
using FFTW

function create_sine_wave(freq, res, filename)
    a = [sin(freq*2*pi*i/res) for i = 1:res]
    writedlm(filename, a)
end

function create_sine_fft(freq, res, filename)
    a = [sin(freq*2*pi*i/res) for i = 1:res]
    writedlm(filename, abs.(fft(a)))
end

function create_3_sine_waves(freq1, freq2, freq3, res, filename)
    a = [sin(freq1*2*pi*i/res) for i = 1:res]
    a += [sin(freq2*2*pi*i/res) for i = 1:res]
    a += [sin(freq3*2*pi*i/res) for i = 1:res]
    writedlm(filename, a)
end

function create_3_sine_fft(freq1, freq2, freq3, res, filename)
    a = [sin(freq1*2*pi*i/res) for i = 1:res]
    a += [sin(freq2*2*pi*i/res) for i = 1:res]
    a += [sin(freq3*2*pi*i/res) for i = 1:res]
    writedlm(filename, abs.(fft(a)))
end

max_stars_repo_name: leios/thesis

max_issues_repo_name: leios/thesis

max_forks_repo_name: leios/thesis


text:
#!/usr/bin/env julia

using Luxor, Colors

include("julia-logo.jl")

function spiral_logo()
    gsave()
    scale(.3, .3)
    r = 200
    setcolor("gray")
    for i in 0:pi/8:2pi
        gsave()
        translate(r * cos(i), r * sin(i))
        rotate(i)
        julialogo()
        grestore()
    end
    grestore()
end

function expandingspiral()
    gsave()
    scale(.3, .3)
    r = 200
    for i in pi:pi/12:6pi
        gsave()
        translate(i/3 * r * cos(i), i/3 * r * sin(i))
        scale(0.8, 0.8)
        rotate(i)
        julialogo()
        grestore()
    end
    grestore()
end

function dropshadow()
    steps=20
    # white-gray ramp
    gramp = linspace(colorant"white", colorant"gray60", steps)
    gsave()
    r = 200
    setopacity(0.1)
    for i in 1:steps
        sethue(gramp[i])
        translate(-0.6,-0.5)
        julialogo(false)
    end
    julialogo()
    grestore()
end

function colorgrid()
    #cols = colormap("RdBu", 5; mid=0.5, logscale=false)
    #cols = sequential_palette(rand(10:360), 5, b=0.1)
    cols = distinguishable_colors(25)
    gsave()
    c = 0
    for row in 100:100:500
        for column in 100:100:500
            gsave()
            setcolor(color(cols[c+=1]))
            translate(row, column)
            scale(0.3, 0.3)
            julialogo(false)
            grestore()
        end
    end
    grestore()
end

function draw_logo(fname)
    Drawing(1600,1600, fname)
    origin()
    background("white")

    translate(-500,-200)
    spiral_logo()

    translate(750,0)
    expandingspiral()

    translate(-1000,500)
    dropshadow()

    translate(700, -100)
    colorgrid()

    finish()
    println("finished test: output in $(fname)")
end

draw_logo("/tmp/julia-logo-draw-eps.eps")


text:
"""
    build a single histogram containing all grads and weight information
"""
function hist_kernel_grad!(h::CuDeviceArray{T,3}, δ𝑤::CuDeviceMatrix{T}, xid::CuDeviceMatrix{UInt8},
    𝑖::CuDeviceVector{S}, 𝑗::CuDeviceVector{S}) where {T,S}

    nbins = size(h, 2)

    it = threadIdx().x
    id, jd = blockDim().x, blockDim().y
    ib, j = blockIdx().x, blockIdx().y
    ig, jg = gridDim().x, gridDim().y

    shared = @cuDynamicSharedMem(T, (size(h, 1), size(h, 2)))
    fill!(shared, 0)
    sync_threads()

    i_tot = length(𝑖)
    iter = 0
    @inbounds while iter * id * ig < i_tot
        i = it + id * (ib - 1) + iter * id * ig
        @inbounds if i <= length(𝑖) && j <= length(𝑗)
            i_idx = 𝑖[i]
            CUDA.atomic_add!(pointer(shared, (xid[i_idx, 𝑗[j]] - 1) * 3 + 1), δ𝑤[1, i_idx])
            CUDA.atomic_add!(pointer(shared, (xid[i_idx, 𝑗[j]] - 1) * 3 + 2), δ𝑤[2, i_idx])
            CUDA.atomic_add!(pointer(shared, (xid[i_idx, 𝑗[j]] - 1) * 3 + 3), δ𝑤[3, i_idx])
        end
        iter += 1
    end
    sync_threads()
    # loop over i blocks
    if it <= nbins
        @inbounds hid = Base._to_linear_index(h, 1, it, 𝑗[j])
        @inbounds CUDA.atomic_add!(pointer(h, hid), shared[1, it])

        @inbounds hid = Base._to_linear_index(h, 2, it, 𝑗[j])
        @inbounds CUDA.atomic_add!(pointer(h, hid), shared[2, it])

        @inbounds hid = Base._to_linear_index(h, 3, it, 𝑗[j])
        @inbounds CUDA.atomic_add!(pointer(h, hid), shared[3, it])
    end
    sync_threads()
    return nothing
end

# base approach - block built along the cols first, the rows (limit collisions)
function update_hist_gpu!(
    ::L,
    h::CuArray{T,3},
    δ𝑤::CuMatrix{T},
    X_bin::CuMatrix{UInt8},
    𝑖::CuVector{S},
    𝑗::CuVector{S}, K;
    MAX_THREADS=256) where {L <: GradientRegression,T,S}

    nbins = size(h, 2)
    thread_i = max(nbins, min(MAX_THREADS, length(𝑖)))
    threads = (thread_i, 1)
    blocks = (8, length(𝑗))
    # println("threads: ", threads, " | blocks: ", blocks)
    @cuda blocks = blocks threads = threads shmem = sizeof(T) * size(h, 1) * nbins hist_kernel_grad!(h, δ𝑤, X_bin, 𝑖, 𝑗)
    CUDA.synchronize()
    return nothing
end

"""
    build a single histogram containing all grads and weight information
"""
function hist_kernel_gauss!(h::CuDeviceArray{T,3}, δ𝑤::CuDeviceMatrix{T}, xid::CuDeviceMatrix{UInt8},
    𝑖::CuDeviceVector{S}, 𝑗::CuDeviceVector{S}) where {T,S}

    nbins = size(h, 2)

    it, k = threadIdx().x, threadIdx().y
    id, jd = blockDim().x, blockDim().y
    ib, j = blockIdx().x, blockIdx().y
    ig, jg = gridDim().x, gridDim().y

    shared = @cuDynamicSharedMem(T, (size(h, 1), size(h, 2)))
    fill!(shared, 0)
    sync_threads()

    i_tot = length(𝑖)
    iter = 0
    @inbounds while iter * id * ig < i_tot
        i = it + id * (ib - 1) + iter * id * ig
        @inbounds if i <= length(𝑖) && j <= length(𝑗)
            i_idx = 𝑖[i]
            CUDA.atomic_add!(pointer(shared, (xid[i_idx, 𝑗[j]] - 1) * 5 + k), δ𝑤[k, i_idx])
        end
        iter += 1
    end
    sync_threads()
    # loop over i blocks
    if it <= nbins
        @inbounds hid = Base._to_linear_index(h, k, it, 𝑗[j])
        @inbounds CUDA.atomic_add!(pointer(h, hid), shared[k, it])
    end
    sync_threads()
    return nothing
end

# base approach - block built along the cols first, the rows (limit collisions)
function update_hist_gpu!(
    ::L,
    h::CuArray{T,3},
    δ𝑤::CuMatrix{T},
    X_bin::CuMatrix{UInt8},
    𝑖::CuVector{S},
    𝑗::CuVector{S}, K;
    MAX_THREADS=128) where {L <: GaussianRegression,T,S}

    nbins = size(h, 2)
    thread_i = max(nbins, min(MAX_THREADS, length(𝑖)))
    threads = (thread_i, 5)
    blocks = (8, length(𝑗))
    @cuda blocks = blocks threads = threads shmem = sizeof(T) * size(h, 1) * nbins hist_kernel_gauss!(h, δ𝑤, X_bin, 𝑖, 𝑗)
    CUDA.synchronize()
    return nothing
end

"""
    Multi-threads split_set!
        Take a view into left and right placeholders. Right ids are assigned at the end of the length of the current node set.
"""
function split_chunk_kernel!(left::CuDeviceVector{S}, right::CuDeviceVector{S}, 𝑖::CuDeviceVector{S}, X_bin, feat, cond_bin, offset, chunk_size, lefts, rights) where {S}

    it = threadIdx().x
    bid = blockIdx().x
    gdim = gridDim().x

    left_count = 0
    right_count = 0

    i = chunk_size * (bid - 1) + 1
    bid == gdim ? bsize = length(𝑖) - chunk_size * (bid - 1) : bsize = chunk_size
    i_max = i + bsize - 1

    @inbounds while i <= i_max
        @inbounds if X_bin[𝑖[i], feat] <= cond_bin
            left_count += 1
            left[offset + chunk_size * (bid - 1) + left_count] = 𝑖[i]
        else
            right_count += 1
            right[offset + chunk_size * (bid - 1) + right_count] = 𝑖[i]
        end
        i += 1
    end
    @inbounds lefts[bid] = left_count
    @inbounds rights[bid] = right_count
    sync_threads()
    return nothing
end

function split_views_kernel!(out::CuDeviceVector{S}, left::CuDeviceVector{S}, right::CuDeviceVector{S}, offset, chunk_size, lefts, rights, sum_lefts, cumsum_lefts, cumsum_rights) where {S}

    it = threadIdx().x
    bid = blockIdx().x
    gdim = gridDim().x

    # bsize = lefts[bid] + rights[bid]
    bid == 1 ? cumsum_left = 0 : cumsum_left = cumsum_lefts[bid-1]
    bid == 1 ? cumsum_right = 0 : cumsum_right = cumsum_rights[bid-1]

    iter = 1
    i_max = lefts[bid]
    @inbounds while iter <= i_max
        out[offset + cumsum_left + iter] = left[offset + chunk_size * (bid - 1) + iter]
        iter += 1
    end

    iter = 1
    i_max = rights[bid]
    @inbounds while iter <= i_max
        out[offset + sum_lefts + cumsum_right + iter] = right[offset + chunk_size * (bid - 1) + iter]
        iter += 1
    end
    sync_threads()
    return nothing
end

function split_set_threads_gpu!(out, left, right, 𝑖, X_bin, feat, cond_bin, offset)
    𝑖_size = length(𝑖)

    nblocks = ceil(Int, min(length(𝑖) / 128, 2^10))
    chunk_size = floor(Int, length(𝑖) / nblocks)

    lefts = CUDA.zeros(Int, nblocks)
    rights = CUDA.zeros(Int, nblocks)

    # threads = 1
    @cuda blocks = nblocks threads = 1 split_chunk_kernel!(left, right, 𝑖, X_bin, feat, cond_bin, offset, chunk_size, lefts, rights)
    CUDA.synchronize()

    sum_lefts = sum(lefts)
    cumsum_lefts = cumsum(lefts)
    cumsum_rights = cumsum(rights)
    @cuda blocks = nblocks threads = 1 split_views_kernel!(out, left, right, offset, chunk_size, lefts, rights, sum_lefts, cumsum_lefts, cumsum_rights)

    CUDA.synchronize()
    return (view(out, offset + 1:offset + sum_lefts), view(out, offset + sum_lefts + 1:offset + length(𝑖)))
end

"""
    update_gains!
        GradientRegression
"""
function update_gains_gpu!(
    loss::L,
    node::TrainNodeGPU{T},
    𝑗::AbstractVector{S},
    params::EvoTypes, K;
    MAX_THREADS=512) where {L<:GradientRegression,T,S}

    cumsum!(node.hL, node.h, dims=2)
    node.hR .= view(node.hL, :, params.nbins:params.nbins, :) .- node.hL
    # cumsum!(view(node.hL, :, :, :), view(node.h, :, :, :), dims=2)
    # cumsum!(view(histR, :, :, :, nid), reverse!(view(hist, :, :, :, nid), dims=2), dims=2)
    # view(node.hR, :, :, :) .= view(node.hL, :, params.nbins:params.nbins, :) .- view(node.hL, :, :, :)

    threads = min(params.nbins, MAX_THREADS)
    blocks = length(𝑗)
    @cuda blocks = blocks threads = threads hist_gains_gpu_kernel!(node.gains, node.hL, node.hR, 𝑗, params.nbins, params.λ, params.min_weight)
    # hist_gains_gpu!(loss, node.gains, node.hL, node.hR, 𝑗, params.nbins, params.λ)
    CUDA.synchronize()
    return nothing
end

function hist_gains_gpu_kernel!(gains::CuDeviceMatrix{T}, hL::CuDeviceArray{T,3}, hR::CuDeviceArray{T,3}, 𝑗::CuDeviceVector{S}, nbins, λ::T, min_𝑤::T) where {T,S}

    i = threadIdx().x
    @inbounds j = 𝑗[blockIdx().x]

    if i == nbins
        gains[i, j] = hL[1, i, j]^2 / (hL[2, i, j] + λ * hL[3, i, j]) / 2
    elseif hL[3, i, j] > min_𝑤 && hR[3, i, j] > min_𝑤
        gains[i, j] = (hL[1, i, j]^2 / (hL[2, i, j] + λ * hL[3, i, j]) +
        hR[1, i, j]^2 / (hR[2, i, j,] + λ * hR[3, i, j])) / 2
    end
    sync_threads()
    return nothing
end

"""
    update_gains!
        GaussianRegression
"""
function update_gains_gpu!(
    loss::L,
    node::TrainNodeGPU{T},
    𝑗::AbstractVector{S},
    params::EvoTypes, K;
    MAX_THREADS=512) where {L<:GaussianRegression,T,S}

    cumsum!(node.hL, node.h, dims=2)
    node.hR .= view(node.hL, :, params.nbins:params.nbins, :) .- node.hL

    thread_i = min(params.nbins, MAX_THREADS)
    threads = thread_i
    blocks = length(𝑗)
    @cuda blocks = blocks threads = threads hist_gains_gpu_kernel_gauss!(node.gains, node.hL, node.hR, 𝑗, params.nbins, params.λ, params.min_weight)
    # hist_gains_gpu!(loss, node.gains, node.hL, node.hR, 𝑗, params.nbins, params.λ)
    CUDA.synchronize()
    return nothing
end

function hist_gains_gpu_kernel_gauss!(gains::CuDeviceMatrix{T}, hL::CuDeviceArray{T,3}, hR::CuDeviceArray{T,3}, 𝑗::CuDeviceVector{S}, nbins, λ::T, min_𝑤::T) where {T,S}

    i = threadIdx().x
    @inbounds j = 𝑗[blockIdx().x]

    if i == nbins
        gains[i, j] = (hL[1, i, j]^2 / (hL[3, i, j] + λ * hL[5, i, j]) + hL[2, i, j]^2 / (hL[4, i, j] + λ * hL[5, i, j])) / 2
    elseif hL[5, i, j] > min_𝑤 && hR[5, i, j] > min_𝑤
        gains[i, j] = (hL[1, i, j]^2 / (hL[3, i, j] + λ * hL[5, i, j]) +
        hR[1, i, j]^2 / (hR[3, i, j,] + λ * hR[5, i, j])) / 2 +
        (hL[2, i, j]^2 / (hL[4, i, j] + λ * hL[5, i, j]) +
        hR[2, i, j]^2 / (hR[4, i, j,] + λ * hR[5, i, j])) / 2
    end
    sync_threads()
    return nothing
end


text:
# precompile.jl contains precompilation directives for all methods whose compilation
# is triggered by running DataFrames tests and takes more than a given threshold
# of 0.1.
# `precompile` calls are generated using SnoopCompile via:
#     using SnoopCompileCore
#     inf_timing = @snoopi tmin=0.1 include("test/runtests.jl")
#     using SnoopCompile
#     pc = SnoopCompile.parcel(inf_timing)
#     SnoopCompile.write("src/other/precompile_tmp.jl", pc[:DataFrames], always=true)
# and then hand edited. The editing steps are:
# * removing signatures with anonymous functions
# * removing signatures with specific NamedTuples passed as arguments except `x1` as name or kwargs
#   (which is generated internally)
# * changing Int64 to Int (to handle 32-bit architectures correctly)
# * disabling precompilation on Julia older than 1.5
# * run @warnpcfail check for all=true and all=false both on Julia stable and nightly

function precompile(all=false)
    VERSION >= v"1.5" || return nothing

        all || ccall(:jl_generating_output, Cint, ()) == 1 || return nothing

        Base.precompile(Tuple{Aggregate{typeof(std), Nothing},Vector{Union{Missing, Int}},GroupedDataFrame{DataFrame}})
        Base.precompile(Tuple{Core.kwftype(typeof(describe)),NamedTuple{(:cols,), Tuple{Vector{Pair{Symbol, Pair{Function, Symbol}}}}},typeof(describe),DataFrame})
        Base.precompile(Tuple{Core.kwftype(typeof(innerjoin)),NamedTuple{(:on, :makeunique, :validate, :renamecols), Tuple{Symbol, Bool, Pair{Bool, Bool}, Pair{String, String}}},typeof(innerjoin),DataFrame,DataFrame})
        Base.precompile(Tuple{Core.kwftype(typeof(innerjoin)),NamedTuple{(:on, :makeunique, :validate, :renamecols), Tuple{Symbol, Bool, Pair{Bool, Bool}, Pair{String, Symbol}}},typeof(innerjoin),DataFrame,DataFrame})
        Base.precompile(Tuple{Core.kwftype(typeof(innerjoin)),NamedTuple{(:on, :makeunique, :validate, :renamecols), Tuple{Symbol, Bool, Pair{Bool, Bool}, Pair{Symbol, Symbol}}},typeof(innerjoin),DataFrame,DataFrame})
        Base.precompile(Tuple{Core.kwftype(typeof(innerjoin)),NamedTuple{(:on, :makeunique, :validate, :renamecols), Tuple{Vector{Any}, Bool, Pair{Bool, Bool}, Pair{String, Symbol}}},typeof(innerjoin),DataFrame,DataFrame})
        Base.precompile(Tuple{Core.kwftype(typeof(innerjoin)),NamedTuple{(:on, :renamecols), Tuple{Vector{Any}, Pair{String, String}}},typeof(innerjoin),DataFrame,DataFrame})
        Base.precompile(Tuple{Core.kwftype(typeof(innerjoin)),NamedTuple{(:on,), Tuple{Symbol}},typeof(innerjoin),DataFrame,DataFrame})
        Base.precompile(Tuple{Core.kwftype(typeof(leftjoin!)),NamedTuple{(:on, :makeunique, :source), Tuple{Pair{Symbol, Symbol}, Bool, Symbol}},typeof(leftjoin!),DataFrame,DataFrame})
        Base.precompile(Tuple{Core.kwftype(typeof(leftjoin)),NamedTuple{(:on, :makeunique, :validate, :renamecols), Tuple{Vector{Any}, Bool, Pair{Bool, Bool}, Pair{Symbol, String}}},typeof(leftjoin),DataFrame,DataFrame})
        Base.precompile(Tuple{Core.kwftype(typeof(leftjoin)),NamedTuple{(:on, :renamecols, :source), Tuple{Vector{Any}, Pair{String, String}, Symbol}},typeof(leftjoin),DataFrame,DataFrame})
        Base.precompile(Tuple{Core.kwftype(typeof(outerjoin)),NamedTuple{(:on, :makeunique, :validate, :renamecols), Tuple{Vector{Any}, Bool, Pair{Bool, Bool}, Pair{Symbol, Symbol}}},typeof(outerjoin),DataFrame,DataFrame})
        Base.precompile(Tuple{Core.kwftype(typeof(outerjoin)),NamedTuple{(:on, :source, :makeunique), Tuple{Symbol, String, Bool}},typeof(outerjoin),DataFrame,DataFrame})
        Base.precompile(Tuple{Core.kwftype(typeof(outerjoin)),NamedTuple{(:on, :source, :makeunique), Tuple{Symbol, Symbol, Bool}},typeof(outerjoin),DataFrame,DataFrame})
        Base.precompile(Tuple{Core.kwftype(typeof(reduce)),NamedTuple{(:cols, :source), Tuple{Symbol, Nothing}},typeof(reduce),typeof(vcat),NTuple{5, DataFrame}})
        Base.precompile(Tuple{Core.kwftype(typeof(select!)),NamedTuple{(:renamecols,), Tuple{Bool}},typeof(select!),GroupedDataFrame{SubDataFrame{DataFrame, SubIndex{Index, Vector{Int}, Vector{Int}}, UnitRange{Int}}},Union{Regex, AbstractString, Function, Signed, Symbol, Unsigned, Pair, Type, All, Between, Cols, InvertedIndex, AbstractVecOrMat{T} where T}})
        Base.precompile(Tuple{Core.kwftype(typeof(show)),NamedTuple{(:allrows, :allcols, :allgroups, :rowlabel, :summary, :truncate), Tuple{Bool, Bool, Bool, Symbol, Bool, Int}},typeof(show),Base.PipeEndpoint,GroupedDataFrame{DataFrame}})
        Base.precompile(Tuple{Core.kwftype(typeof(show)),NamedTuple{(:show_row_number,), Tuple{Bool}},typeof(show),IOContext{IOBuffer},DataFrame})
        Base.precompile(Tuple{Core.kwftype(typeof(vcat)),NamedTuple{(:cols, :source), Tuple{Symbol, Symbol}},typeof(vcat),DataFrame,Vararg{DataFrame}})
        Base.precompile(Tuple{Core.kwftype(typeof(vcat)),NamedTuple{(:cols,), Tuple{Vector{Symbol}}},typeof(vcat),DataFrame,Vararg{DataFrame}})
        Base.precompile(Tuple{Type{DataFrame},Vector{Tuple{Int, Int}}})
        Base.precompile(Tuple{typeof(_innerjoin_sorted),OnCol{Tuple{Vector{Int}, Vector{String}}, 2},OnCol{Tuple{SubArray{Union{Missing, Int}, 1, Vector{Union{Missing, Int}}, Tuple{Vector{Int}}, false}, SubArray{Union{Missing, String}, 1, Vector{Union{Missing, String}}, Tuple{Vector{Int}}, false}}, 2}})
        Base.precompile(Tuple{typeof(_innerjoin_sorted),OnCol{Tuple{Vector{UInt32}, Vector{UInt32}, Vector{UInt32}}, 3},OnCol{Tuple{Vector{UInt32}, Vector{UInt32}, Vector{UInt32}}, 3}})
        Base.precompile(Tuple{typeof(_innerjoin_unsorted_int),SubArray{Union{Missing, Int}, 1, Vector{Union{Missing, Int}}, Tuple{Vector{Int}}, false},Vector{Missing}})
        Base.precompile(Tuple{typeof(_innerjoin_unsorted_int),Vector{Int32},Vector{Int32}})
        Base.precompile(Tuple{typeof(_innerjoin_unsorted),OnCol{Tuple{Vector{Union{Missing, String}}, Vector{Union{Missing, String}}, Vector{Union{Missing, String}}}, 3},OnCol{Tuple{Vector{String}, Vector{String}, Vector{String}}, 3}})
        Base.precompile(Tuple{typeof(_innerjoin_unsorted),OnCol{Tuple{Vector{Union{Missing, Symbol}}, Vector{Union{Missing, Symbol}}}, 2},OnCol{Tuple{Vector{Union{Missing, Symbol}}, Vector{Union{Missing, Symbol}}}, 2}})
        Base.precompile(Tuple{typeof(_mean_fast),Vector{Vector{Union{Missing, Int}}}})
        Base.precompile(Tuple{typeof(_semijoin_unsorted_int),Vector{Int},Vector{Union{Missing, Int}},BitVector,Bool})
        Base.precompile(Tuple{typeof(_semijoin_unsorted),OnCol{Tuple{Vector{UInt32}, Vector{UInt32}}, 2},OnCol{Tuple{Vector{Union{Missing, UInt32}}, Vector{Union{Missing, UInt32}}}, 2},BitVector,Bool})
        Base.precompile(Tuple{typeof(Base.Broadcast.dotview),SubDataFrame{DataFrame, Index, Vector{Int}},Colon,Symbol})
        Base.precompile(Tuple{typeof(describe),DataFrame})
        Base.precompile(Tuple{typeof(do_call),ComposedFunction{typeof(sum), typeof(skipmissing)},Vector{Int},Vector{Int},Vector{Int},GroupedDataFrame{DataFrame},Tuple{Vector{Union{Missing, String}}},Int})
        Base.precompile(Tuple{typeof(do_call),typeof(cor),Vector{Int},Vector{Int},Vector{Int},GroupedDataFrame{DataFrame},Tuple{Vector{Int}, Vector{Int}},Int})
        Base.precompile(Tuple{typeof(do_call),typeof(sum),Vector{Int},Vector{Int},Vector{Int},GroupedDataFrame{DataFrame},Tuple{Vector{DataFrame}},Int})
        Base.precompile(Tuple{typeof(groupby),DataFrame,Symbol})
        Base.precompile(Tuple{typeof(names),DataFrameRow{DataFrame, SubIndex{Index, Vector{Int}, Vector{Int}}},Cols{Tuple{String}}})
        Base.precompile(Tuple{typeof(permutedims),DataFrame,Symbol})
        Base.precompile(Tuple{typeof(prepare_on_col),Vector{Int},Vector{Union{}}})
        Base.precompile(Tuple{typeof(repeat),DataFrame,Int})
        Base.precompile(Tuple{typeof(row_group_slots),Tuple{BitVector},Val{false},Vector{Int},Bool,Nothing})
        Base.precompile(Tuple{typeof(row_group_slots),Tuple{PooledVector{Int, UInt32, Vector{UInt32}}, PooledVector{String, UInt32, Vector{UInt32}}},Val{false},Vector{Int},Bool,Nothing})
        Base.precompile(Tuple{typeof(row_group_slots),Tuple{PooledVector{String, UInt32, Vector{UInt32}}, Vector{Int}},Val{false},Vector{Int},Bool,Nothing})
        Base.precompile(Tuple{typeof(row_group_slots),Tuple{PooledVector{String, UInt32, Vector{UInt32}}, Vector{Int32}},Val{false},Vector{Int},Bool,Nothing})
        Base.precompile(Tuple{typeof(row_group_slots),Tuple{PooledVector{Union{Missing, Int}, UInt32, Vector{UInt32}}, PooledVector{String, UInt32, Vector{UInt32}}},Val{false},Vector{Int},Bool,Nothing})
        Base.precompile(Tuple{typeof(row_group_slots),Tuple{PooledVector{Union{Missing, String}, UInt32, Vector{UInt32}}, PooledVector{String, UInt32, Vector{UInt32}}},Val{false},Vector{Int},Bool,Nothing})
        Base.precompile(Tuple{typeof(row_group_slots),Tuple{RepeatedVector{Union{Missing, Int}}, RepeatedVector{Union{Missing, Int}}, RepeatedVector{Union{Missing, String}}, Vector{Union{Missing, Int}}},Val{false},Vector{Int},Bool,Nothing})
        Base.precompile(Tuple{typeof(row_group_slots),Tuple{SubArray{Int, 1, Vector{Int}, Tuple{Base.OneTo{Int}}, true}},Tuple{IntegerRefpool{Int}},Tuple{IntegerRefarray{SubArray{Int, 1, Vector{Int}, Tuple{Base.OneTo{Int}}, true}}},Val{false},Vector{Int},Bool,Bool})
        Base.precompile(Tuple{typeof(row_group_slots),Tuple{SubArray{Union{Missing, Int}, 1, Vector{Union{Missing, Int}}, Tuple{Base.OneTo{Int}}, true}},Tuple{IntegerRefpool{Union{Missing, Int}}},Tuple{IntegerRefarray{SubArray{Union{Missing, Int}, 1, Vector{Union{Missing, Int}}, Tuple{Base.OneTo{Int}}, true}}},Val{false},Vector{Int},Bool,Bool})
        Base.precompile(Tuple{typeof(row_group_slots),Tuple{Vector{Bool}},Val{false},Vector{Int},Bool,Nothing})
        Base.precompile(Tuple{typeof(row_group_slots),Tuple{Vector{Float64}, Vector{Int}},Val{false},Vector{Int},Bool,Nothing})
        Base.precompile(Tuple{typeof(row_group_slots),Tuple{Vector{Int}, Vector{Int}, Vector{Int}},Val{false},Vector{Int},Bool,Nothing})
        Base.precompile(Tuple{typeof(row_group_slots),Tuple{Vector{Int}, Vector{Union{Missing, Int}}},Val{false},Vector{Int},Bool,Nothing})
        Base.precompile(Tuple{typeof(row_group_slots),Tuple{Vector{Int32}, Vector{Int32}},Val{false},Vector{Int},Bool,Nothing})
        Base.precompile(Tuple{typeof(row_group_slots),Tuple{Vector{Int32}},Tuple{IntegerRefpool{Int}},Tuple{IntegerRefarray{Vector{Int32}}},Val{false},Vector{Int},Bool,Bool})
        Base.precompile(Tuple{typeof(row_group_slots),Tuple{Vector{Int32}},Val{false},Vector{Int},Bool,Nothing})
        Base.precompile(Tuple{typeof(row_group_slots),Tuple{Vector{String}, Vector{Int32}},Val{false},Vector{Int},Bool,Nothing})
        Base.precompile(Tuple{typeof(row_group_slots),Tuple{Vector{Union{Missing, BigInt}}},Val{false},Vector{Int},Bool,Nothing})
        Base.precompile(Tuple{typeof(row_group_slots),Tuple{Vector{Union{Missing, Bool}}},Tuple{IntegerRefpool{Union{Missing, Int}}},Tuple{IntegerRefarray{Vector{Union{Missing, Bool}}}},Val{false},Vector{Int},Bool,Bool})
        Base.precompile(Tuple{typeof(row_group_slots),Tuple{Vector{Union{Missing, Bool}}},Val{false},Vector{Int},Bool,Nothing})
        Base.precompile(Tuple{typeof(row_group_slots),Tuple{Vector{Union{Missing, Float64}}, Vector{Float64}},Val{false},Vector{Int},Bool,Nothing})
        Base.precompile(Tuple{typeof(row_group_slots),Tuple{Vector{Union{Missing, Float64}}, Vector{Int}},Val{false},Vector{Int},Bool,Nothing})
        Base.precompile(Tuple{typeof(row_group_slots),Tuple{Vector{Union{Missing, Int}}, Vector{Union{Missing, Float64}}},Tuple{IntegerRefpool{Union{Missing, Int}}, IntegerRefpool{Union{Missing, Int}}},Tuple{IntegerRefarray{Vector{Union{Missing, Int}}}, IntegerRefarray{Vector{Union{Missing, Float64}}}},Val{false},Vector{Int},Bool,Bool})
        Base.precompile(Tuple{typeof(row_group_slots),Tuple{Vector{Union{Missing, Int}}},Val{false},Vector{Int},Bool,Nothing})
        Base.precompile(Tuple{typeof(select),SubDataFrame{DataFrame, Index, UnitRange{Int}},Any})
        Base.precompile(Tuple{typeof(select),SubDataFrame{DataFrame, SubIndex{Index, UnitRange{Int}, UnitRange{Int}}, Vector{Int}},Any,Any,Vararg{Any}})
        Base.precompile(Tuple{typeof(setindex!),DataFrame,Vector{Int},Vector{Bool},Int})
        Base.precompile(Tuple{typeof(show),IOBuffer,DataFrameRow{DataFrame, SubIndex{Index, UnitRange{Int}, UnitRange{Int}}}})
        Base.precompile(Tuple{typeof(sort),SubDataFrame{DataFrame, SubIndex{Index, UnitRange{Int}, UnitRange{Int}}, UnitRange{Int}}})
        Base.precompile(Tuple{typeof(subset),GroupedDataFrame{SubDataFrame{DataFrame, Index, UnitRange{Int}}},Any})
        Base.precompile(Tuple{typeof(transform),GroupedDataFrame{DataFrame},Union{Regex, AbstractString, Function, Signed, Symbol, Unsigned, Pair, Type, All, Between, Cols, InvertedIndex, AbstractVecOrMat{T} where T}})
        Base.precompile(Tuple{typeof(view),SubDataFrame{DataFrame, Index, UnitRange{Int}},Int,Int})

    if all
        for v in ([1, 2], [2, 1], [2, 2, 1], Int32[1, 2], Int32[2, 1], Int32[2, 2, 1]),
            op in (identity, x -> string.(x), x -> PooledArrays.PooledArray(string.(x))),
            on in (:v1, [:v1, :v2])
            df = DataFrame(v1=op(v), v2=v)
            combine(groupby(df, on), identity, :v1 => identity,
                    :v2 => ByRow(identity), :v2 => sum)
            innerjoin(df, select(df, on), on=on)
            outerjoin(df, select(df, on), on=on)
        end
    end
    return nothing
end


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    @assert precompile(Scene, ())
    @assert precompile(update_limits!, (Scene,))
    @assert precompile(update_limits!, (Scene, Automatic))
    @assert precompile(update_limits!, (Scene, FRect3D))
end


text:
let
    p = FunctionPolicy(x::Bool->!x)
    @test action(p, true) == false

    s = FunctionSolver(x::Int->2*x)
    p = solve(s, LegacyGridWorld())
    @test action(p, 10) == 20
    @test action(p, 100) == 200
    updater(p) # just to make sure this doesn't error
end

max_stars_repo_name: JuliaPOMDP/POMDPPolicies

max_issues_repo_name: JuliaPOMDP/POMDPPolicies

max_forks_repo_name: JuliaPOMDP/POMDPPolicies


text:
module Transformation

export pearsonresiduals

    include("pearsonresiduals.jl")

end


text:
module Init_Test

using GLFW, FactCheck

facts() do

	@fact GLFW.Init() --> true

end

facts() do

	window = GLFW.CreateWindow(640, 480, "GLFW.jl")
	@fact typeof(window) --> GLFW.Window
	GLFW.DestroyWindow(window)
end

#Initialize debugging window
facts() do

	GLFW.WindowHint(GLFW.OPENGL_DEBUG_CONTEXT, Cint(true))
	window = GLFW.CreateWindow(640, 480, "GLFW.jl")
	@fact typeof(window) --> GLFW.Window
	GLFW.DestroyWindow(window)

end

#Initialize with regular hints
facts() do

	GLFW.DefaultWindowHints()
	window = GLFW.CreateWindow(640, 480, "GLFW.jl")
	@fact typeof(window) --> GLFW.Window
	GLFW.DestroyWindow(window)

end

#Check OpenGL Versions
facts() do

	for i in [1,2]
		for j in [0,1]

			GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, i)
			GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, j)
			GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, convert(Cuint,0))
			GLFW.WindowHint(GLFW.OPENGL_PROFILE,GLFW.OPENGL_ANY_PROFILE)
			window = GLFW.CreateWindow(640, 480, "GLFW.jl")
			@fact typeof(window) --> GLFW.Window
			GLFW.DestroyWindow(window)

		end
	end
end

#Version 3.0
facts() do

	GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
	GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 0)
	GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, convert(Cuint,0))
	GLFW.WindowHint(GLFW.OPENGL_PROFILE,GLFW.OPENGL_ANY_PROFILE)
	window = GLFW.CreateWindow(640, 480, "GLFW.jl")
	@fact typeof(window) --> GLFW.Window
	GLFW.DestroyWindow(window)
end

#Version 3.1
#=
facts() do

	GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
	GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 1)
	GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, convert(Cuint,0))
	GLFW.WindowHint(GLFW.OPENGL_PROFILE,GLFW.OPENGL_ANY_PROFILE)
	window = GLFW.CreateWindow(640, 480, "GLFW.jl")
	@fact typeof(window) --> GLFW.Window
	GLFW.DestroyWindow(window)
end
=#

#Version 3.3 FC
facts() do

	GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
	GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 3)
	GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, convert(Cuint,1))
	GLFW.WindowHint(GLFW.OPENGL_PROFILE,GLFW.OPENGL_CORE_PROFILE)
	window = GLFW.CreateWindow(640, 480, "GLFW.jl")
	@fact typeof(window) --> GLFW.Window
	GLFW.DestroyWindow(window)
end

end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct TGTypeII <: TurbineGov
        R::Float64
        T1::Float64
        T2::Float64
        τ_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        P_ref::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int
        internal::InfrastructureSystemsInternal
    end

Parameters of a Turbine Governor Type II.

# Arguments
- `R::Float64`: Droop parameter, validation range: `(0, nothing)`
- `T1::Float64`: Transient gain time constant, validation range: `(0, nothing)`
- `T2::Float64`: Power fraction time constant, validation range: `(0, nothing)`
- `τ_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Power into the governor limits
- `P_ref::Float64`: Reference Power Set-point, validation range: `(0, nothing)`
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`: The states of the TGTypeI model are:
	x_g1: lead-lag state
- `n_states::Int`: TGTypeII has 1 state
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct TGTypeII <: TurbineGov
    "Droop parameter"
    R::Float64
    "Transient gain time constant"
    T1::Float64
    "Power fraction time constant"
    T2::Float64
    "Power into the governor limits"
    τ_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Reference Power Set-point"
    P_ref::Float64
    ext::Dict{String, Any}
    "The states of the TGTypeI model are:
	x_g1: lead-lag state"
    states::Vector{Symbol}
    "TGTypeII has 1 state"
    n_states::Int
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function TGTypeII(R, T1, T2, τ_limits, P_ref=1.0, ext=Dict{String, Any}(), )
    TGTypeII(R, T1, T2, τ_limits, P_ref, ext, [:xg], 1, InfrastructureSystemsInternal(), )
end

function TGTypeII(; R, T1, T2, τ_limits, P_ref=1.0, ext=Dict{String, Any}(), states=[:xg], n_states=1, internal=InfrastructureSystemsInternal(), )
    TGTypeII(R, T1, T2, τ_limits, P_ref, ext, states, n_states, internal, )
end

# Constructor for demo purposes; non-functional.
function TGTypeII(::Nothing)
    TGTypeII(;
        R=0,
        T1=0,
        T2=0,
        τ_limits=(min=0.0, max=0.0),
        P_ref=0,
        ext=Dict{String, Any}(),
    )
end

"""Get [`TGTypeII`](@ref) `R`."""
get_R(value::TGTypeII) = value.R
"""Get [`TGTypeII`](@ref) `T1`."""
get_T1(value::TGTypeII) = value.T1
"""Get [`TGTypeII`](@ref) `T2`."""
get_T2(value::TGTypeII) = value.T2
"""Get [`TGTypeII`](@ref) `τ_limits`."""
get_τ_limits(value::TGTypeII) = value.τ_limits
"""Get [`TGTypeII`](@ref) `P_ref`."""
get_P_ref(value::TGTypeII) = value.P_ref
"""Get [`TGTypeII`](@ref) `ext`."""
get_ext(value::TGTypeII) = value.ext
"""Get [`TGTypeII`](@ref) `states`."""
get_states(value::TGTypeII) = value.states
"""Get [`TGTypeII`](@ref) `n_states`."""
get_n_states(value::TGTypeII) = value.n_states
"""Get [`TGTypeII`](@ref) `internal`."""
get_internal(value::TGTypeII) = value.internal

"""Set [`TGTypeII`](@ref) `R`."""
set_R!(value::TGTypeII, val) = value.R = val
"""Set [`TGTypeII`](@ref) `T1`."""
set_T1!(value::TGTypeII, val) = value.T1 = val
"""Set [`TGTypeII`](@ref) `T2`."""
set_T2!(value::TGTypeII, val) = value.T2 = val
"""Set [`TGTypeII`](@ref) `τ_limits`."""
set_τ_limits!(value::TGTypeII, val) = value.τ_limits = val
"""Set [`TGTypeII`](@ref) `P_ref`."""
set_P_ref!(value::TGTypeII, val) = value.P_ref = val
"""Set [`TGTypeII`](@ref) `ext`."""
set_ext!(value::TGTypeII, val) = value.ext = val


text:
export findComponentsWithSignifFreq, getMotionFreq, findPeaks,
       bins2hertz, hertz2bins, generateSpectrum,
       parabelFit, parabelMax

function findComponentsWithSignifFreq(u)
  uFD = zeros(size(u,1),Int64(size(u,2)/2))
  for i = 1:size(u,1)
    uFD[i,:] = generateSpectrum(squeeze(u[i,:]))
  end
  value = zeros(size(uFD,1))
  for i = 1:size(uFD,1)
    value[i] = maximum(squeeze(uFD[i,:]))
  end
  value[1] = 0
  return findmax(value)
end

"""
    getMotionFreq(bSF::MPIFile, bMeas::MPIFile, choosePeak::Int)

Determines frequency [Hz] of signal modulation for each patch and frame

...
- `bSF::MPIFile`: Systemfunction ==> required because the motion frequency is only determined from components with high SNR
- `bMeas::MPIFile`: measurement
- `choosePeak::Int`: Within the measurement signal different modulations can be found and also higher harmonics.
 Choose the number of the peak you want
...
"""
function getMotionFreq(bSF::MPIFile, bMeas::MPIFile, choosePeak::Int,
		       SNRThreshMF = 10, minFreqMF = 80e3, recChannelsMF = [2]) # frequency selection

    freqs = filterFrequencies(bSF,SNRThresh=SNRThreshMF,minFreq=minFreqMF,recChannels=recChannelsMF)
    repFreq = 1 / dfCycle(bMeas);
    peaksInHz = 0
    selectedPeakInHz = zeros(acqNumFrames(bMeas),acqNumPatches(bMeas))
    uspec = getMeasurementsFD(bMeas, frequencies=freqs, numAverages=1, spectralLeakageCorrection=false)
    norm, ind = findComponentsWithSignifFreq(squeeze(uspec[:,1:acqNumPeriodsPerPatch(bMeas),1]))
    @debug "Maximum freq" ind

	for j=1:acqNumFrames(bMeas)
    	for i = 1:acqNumPatches(bMeas)
        	uFT = generateSpectrum(squeeze(uspec[ind,(i-1)*acqNumPeriodsPerPatch(bMeas)+1:i*acqNumPeriodsPerPatch(bMeas),j]))
        	frind = collect(1:acqNumPeriodsPerPatch(bMeas)/2)
       		frHz = bins2hertz(frind,repFreq,acqNumPeriodsPerPatch(bMeas))

            peaksInInt = findPeaks(uFT,stdTh=1)
            peaksInHz = bins2hertz(peaksInInt, 1/dfCycle(bMeas), acqNumPeriodsPerPatch(bMeas))
            @debug "Found peaks: " peaksInHz

			selectedPeak = parabelMax(parabelFit(peaksInInt[choosePeak]-1:peaksInInt[choosePeak]+1,log.(uFT[peaksInInt[choosePeak]-1:peaksInInt[choosePeak]+1])))
			selectedPeakInHz[j,i] = bins2hertz(selectedPeak,repFreq,acqNumPeriodsPerPatch(bMeas))
		end
    end
    return selectedPeakInHz
end

function parabelFit(x,y)
    p = [x.^2 x ones(x)]\y
    return p
end

function parabelMax(par)
    return -par[2]/2/par[1]
end

function generateSpectrum(data)
  data = applyHannWindow(broadcast(-, data, mean(data)))
  ftc = fft(data)

  ft = ftc
  ft = abs.(ftc[2:div(end,2)+1])
  helpv = size(ft,1)
  ft[1:end] += abs.(ftc[end:-1:end-helpv+1])

  return ft
end

function applyHannWindow(data)

    N = size(data,1)
    n = (0:N-1)

    return data .* (0.5*( 1 .-cos.(2*pi*n/(N-1)) ))
end

function reduceHighSpectralValuesToPeaks(data, isHigh)
  diracpeaks = []
  idx = 1
  while !isempty(findall(isHigh[idx:end]))
      from = findfirst(isHigh[idx:end]) + idx - 1
      till = findfirst(isHigh[from+1:end].==0) + from - 1
      push!(diracpeaks, findmax(data[from:till])[2]+from-1)
      idx = till + 1
  end
  return diracpeaks
end

function enforceMinDistanceOfPeaks(diracpeaks, minBinDist, data)
  keep = false
  while sum(.!keep) > 0
      keep = ( diracpeaks[2:end]-diracpeaks[1:end-1] ) .> minBinDist
      push!(keep,true) # append last index
      # keep the greater ones
      rem = findall(.!keep)
      @debug rem
      for i in rem
	  @debug diracpeaks[i]
          if data[diracpeaks[i]] > data[diracpeaks[i+1]]
              keep[i] = true
              keep[i+1] = false
          end
      end
      diracpeaks = diracpeaks[keep .== true]

  end
  return diracpeaks
end

"""
	findPeaks(data; stdTh=1.5, minBinDist=0, minBin=1)
Find peaks in the data.

- data:         column vector of data
- stdTh:        peaks are at least stdTh*standarddeviation above the mean
- minBinDist:   min number of bins between two peaks
- minBin:       no peaks in bins below minBin
- return:       peak bins
"""
function findPeaks(data; stdTh=1.5, minBinDist=0, minBin=1)
    # smooth data??

    if length(findall(data[1:end-1].<data[2:end])) == 0 || length(findall(data[1:end-1].>data[2:end])) == 0
	warn("This should not happen... No variation of data")
	return []
    end

    # exclude beginning and end of data
    frst = maximum( [findfirst( data[1:end-1].<data[2:end] ), ceil(Int32, minBin)] )
    lst = minimum( [findlast( data[1:end-1].>data[2:end] ), floor(Int32, length(data)-minBin+1)] ) # last entry of data never used
    data = data[frst:lst]

    # detrending
    #x = 1:length(data)
    #linPar = [x ones(x)]\data
    #data = data - (linPar[1]*x+linPar[2])

    # mean,std
    mu = mean(data)
    sigma = std(data)

    # threshold on diff to mean
    peaks = data .> (mu + stdTh*sigma)
    if isempty(peaks)
        @error "No peaks were found!"
    end

    # no maximum at last position
    peaks[end] = false

    diracpeaks = reduceHighSpectralValuesToPeaks(data,peaks)
    diracpeaks = enforceMinDistanceOfPeaks(diracpeaks,minBinDist,data)

    diracpeaks = diracpeaks .+ frst.-1
    return diracpeaks
end

"""
	bins2hertz(bin, fs, nSamples)
Convert FFT bin index to Hertz.

- bin:      bin index
- fs:       sampling frequency
- nSamples: number of samples FFT was applied to
- return:   Hertz
"""
function bins2hertz(bin, fs, nSamples)
  return fs/nSamples*bin
end

"""
	hertz2bins(hz, fs, nSamples)
Convert Hertz to FFT bin index.

- hz:       Hertz
- fs:       sampling frequency
- nSamples: number of samples FFT was applied to
- return:   bins
"""
function hertz2bins(hz, fs, nSamples)
  return nSamples/fs*hz
end


text:
not_implemented_for(x) = error("Not implemented for $(x). You might want to put:  `using Makie` into your code!")

#TODO only have one?
const Theme = Attributes

default_theme(scene, T) = Attributes()

function default_theme(scene)
    light = Vec3f0[Vec3f0(1.0,1.0,1.0), Vec3f0(0.1,0.1,0.1), Vec3f0(0.9,0.9,0.9), Vec3f0(20,20,20)]
    Theme(
        color = theme(scene, :color),
        visible = theme(scene, :visible),
        linewidth = 1,
        light = light,
        transformation = automatic,
        model = automatic,
        alpha = 1.0,
    )
end

"""
    `image(x, y, image)` / `image(image)`

Plots an image on range `x, y` (defaults to dimensions).
"""
@atomic(Image) do scene
    Theme(;
        default_theme(scene)...,
        colormap = [RGBAf0(0,0,0,1), RGBAf0(1,1,1,1)],
        colorrange = automatic,
        fxaa = false,
    )
end

# could be implemented via image, but might be optimized specifically by the backend
"""
    `heatmap(x, y, values)` or `heatmap(values)`

Plots a heatmap as an image on `x, y` (defaults to interpretation as dimensions).
"""
@atomic(Heatmap) do scene
    Theme(;
        default_theme(scene)...,
        colormap = theme(scene, :colormap),
        colorrange = automatic,
        linewidth = 0.0,
        levels = 1,
        fxaa = true,
        interpolate = false
    )
end

"""
    `volume(volume_data)`

Plots a volume. Available algorithms are:
* `:iso` => IsoValue
* `:absorption` => Absorption
* `:mip` => MaximumIntensityProjection
* `:absorptionrgba` => AbsorptionRGBA
* `:indexedabsorption` => IndexedAbsorptionRGBA
"""
@atomic(Volume) do scene
    Theme(;
        default_theme(scene)...,
        fxaa = true,
        algorithm = :iso,
        absorption = 1f0,
        isovalue = 0.5f0,
        isorange = 0.05f0,
        colormap = theme(scene, :colormap),
        colorrange = (0, 1)
    )
end

"""
    `surface(x, y, z)`

Plots a surface, where `(x, y, z)` are supposed to lie on a grid.
"""
@atomic(Surface) do scene
    Theme(;
        default_theme(scene)...,
        colormap = theme(scene, :colormap),
        colorrange = automatic,
        fxaa = true,
    )
end

"""
    `lines(x, y, z)` / `lines(x, y)` / or `lines(positions)`

Creates a connected line plot for each element in `(x, y, z)`, `(x, y)` or `positions`.
"""
@atomic(Lines) do scene
    Theme(;
        default_theme(scene)...,
        linewidth = 1.0,
        colormap = theme(scene, :colormap),
        colorrange = automatic,
        linestyle = nothing,
        fxaa = false
    )
end

"""
    `linesegments(x, y, z)` / `linesegments(x, y)` / `linesegments(positions)`

Plots a line for each pair of points in `(x, y, z)`, `(x, y)`, or `positions`.

**Attributes**:
The same as for [`lines`](@ref)
"""
@atomic(LineSegments) do scene
    default_theme(scene, Lines)
end

# alternatively, mesh3d? Or having only mesh instead of poly + mesh and figure out 2d/3d via dispatch
"""
    `mesh(x, y, z)`, `mesh(mesh_object)`, `mesh(x, y, z, faces)`, or `mesh(xyz, faces)`

Plots a 3D mesh.
"""
@atomic(Mesh) do scene
    Theme(;
        default_theme(scene)...,
        fxaa = true,
        interpolate = false,
        shading = true,
        colormap = theme(scene, :colormap),
        colorrange = automatic,
    )
end

"""
    `scatter(x, y, z)` / `scatter(x, y)` / `scatter(positions)`

Plots a marker for each element in `(x, y, z)`, `(x, y)`, or `positions`.
"""
@atomic(Scatter) do scene
    Theme(;
        default_theme(scene)...,
        marker = Circle,
        markersize = 0.1,
        strokecolor = RGBA(0, 0, 0, 0),
        strokewidth = 0.0,
        glowcolor = RGBA(0, 0, 0, 0),
        glowwidth = 0.0,
        rotations = Billboard(),
        colormap = theme(scene, :colormap),
        colorrange = automatic,
        marker_offset = automatic,
        fxaa = false,
        transform_marker = false, # Applies the plots transformation to marker
        uv_offset_width = Vec4f0(0),
        distancefield = nothing
    )
end

"""
    `meshscatter(x, y, z)` / `meshscatter(x, y)` / `meshscatter(positions)`

Plots a mesh for each element in `(x, y, z)`, `(x, y)`, or `positions` (similar to `scatter`).
`markersize` is a scaling applied to the primitive passed as `marker`
"""
@atomic(MeshScatter) do scene
    Theme(;
        default_theme(scene)...,
        marker = Sphere(Point3f0(0), 1f0),
        markersize = 0.1,
        rotations = Quaternionf0(0, 0, 0, 1),
        colormap = theme(scene, :colormap),
        colorrange = automatic,
        fxaa = true
    )
end

"""
    `text(string)`

Plots a text.
"""
@atomic(Text) do scene
    Theme(;
        default_theme(scene)...,
        strokecolor = (:black, 0.0),
        strokewidth = 0,
        font = theme(scene, :font),
        align = (:left, :bottom),
        rotation = 0.0,
        textsize = 20,
        position = Point2f0(0),
    )
end

const atomic_function_symbols = (
        :text, :meshscatter, :scatter, :mesh, :linesegments,
        :lines, :surface, :volume, :heatmap, :image
)

const atomic_functions = getfield.(Ref(AbstractPlotting), atomic_function_symbols)

function color_and_colormap!(plot, intensity = plot[:color])
    if isa(intensity[], AbstractArray{<: Number})
        haskey(plot, :colormap) || error("Plot $T needs to have a colormap to allow the attribute color to be an array of numbers")
        replace_automatic!(plot, :colorrange) do
            lift(extrema_nan, intensity)
        end
        true
    else
        delete!(plot, :colorrange)
        false
    end
end

"""
    `calculated_attributes!(plot::AbstractPlot)`

Fill in values that can only be calculated when we have all other attributes filled
"""
calculated_attributes!(plot::T) where T = calculated_attributes!(T, plot)

"""
    `calculated_attributes!(trait::Type{<: AbstractPlot}, plot)`
trait version of calculated_attributes
"""
calculated_attributes!(trait, plot) = nothing

function calculated_attributes!(::Type{<: Mesh}, plot)
    need_cmap = color_and_colormap!(plot)
    need_cmap || delete!(plot, :colormap)
    return
end

function calculated_attributes!(::Type{<: Union{Heatmap, Image}}, plot)
    plot[:color] = plot[3]
    color_and_colormap!(plot)
end
function calculated_attributes!(::Type{<: Surface}, plot)
    color_and_colormap!(plot, plot[3])
end
function calculated_attributes!(::Type{<: MeshScatter}, plot)
    color_and_colormap!(plot)
end

function calculated_attributes!(::Type{<: Scatter}, plot)
    # calculate base case
    color_and_colormap!(plot)
    replace_automatic!(plot, :marker_offset) do
        # default to middle
        lift(x-> Vec2f0.((x .* (-0.5f0))), plot[:markersize])
    end
end

function calculated_attributes!(::Type{<: Union{Lines, LineSegments}}, plot)
    color_and_colormap!(plot)
end

# # to allow one color per edge
# function calculated_attributes!(plot::LineSegments)
#     plot[:color] = lift(plot[:color], plot[1]) do c, p
#         if (length(p) ÷ 2) == length(c)
#             [c[k] for k in 1:length(c), l in 1:2]
#         else
#             c
#         end
#     end
# end

function (PT::Type{<: Combined})(parent, transformation, attributes, input_args, converted)
    PT(parent, transformation, attributes, input_args, converted, AbstractPlot[])
end
plotsym(::Type{<:AbstractPlot{F}}) where F = Symbol(typeof(F).name.mt.name)

function (PlotType::Type{<: AbstractPlot{Typ}})(scene::SceneLike, attributes::Attributes, args::Tuple) where Typ
    # make sure all arguments are a node
    # with a sensible name
    arg_nodes = node.(ntuple(i-> Symbol("input $i"), length(args)), args)
    args_converted = map(arg_nodes...) do args...
        # do the argument conversion inside a lift
        args = convert_arguments(PlotType, args...)
        PlotType2 = plottype(args...)
        args
    end
    # now get a signal node/signal for each argument
    N = length(value(args_converted))
    names = argument_names(PlotType, N)
    node_args_seperated = ntuple(N) do i
        map(args_converted, name = string(names[i])) do x
            if i <= length(x)
                x[i]
            else
                error("You changed the number of arguments for $PlotType. This isn't allowed!")
            end
        end
    end
    plot_attributes, rest = merged_get!(()-> default_theme(scene, PlotType), plotsym(PlotType), scene, attributes)
    trans = get(plot_attributes, :transformation, automatic)
    transformation = if to_value(trans) == automatic
        Transformation(scene)
    elseif isa(to_value(trans), Transformation)
        to_value(trans)
    else
        t = Transformation(scene)
        transform!(t, to_value(trans))
        t
    end

    replace_automatic!(plot_attributes, :model) do
        transformation.model
    end
    # The argument type of the final plot object is the assumened to stay constant after
    # argument conversion. This might not always hold, but it simplifies
    # things quite a bit
    ArgTyp = typeof(value(args_converted))
    # construct the fully qualified plot type, from the possible incomplete (abstract)
    # PlotType
    FinalType = basetype(PlotType){Typ, ArgTyp}
    # create the plot, with the full attributes, the input signals, and the final signal nodes.
    plot_obj = FinalType(scene, transformation, plot_attributes, arg_nodes, node_args_seperated)
    calculated_attributes!(plot_obj)
    plot_obj, rest
end

"""
    `plot_type(plot_args...)`

The default plot type for any argument is `lines`.
Any custom argument combination that has only one meaningful way to be plotted should overload this.
e.g.:
```example
    # make plot(rand(5, 5, 5)) plot as a volume
    plottype(x::Array{<: AbstractFlot, 3}) = Volume
```
"""
plottype(plot_args...) = Combined{Any, Tuple{typeof.(to_value.(plot_args))...}}

function Plot(args::Vararg{Any, N}) where N
    Combined{Any, <: Tuple{args...}}
end
Base.@pure function Plot(::Type{T}) where T
    Combined{Any, <: Tuple{T}}
end
Base.@pure function Plot(::Type{T1}, ::Type{T2}) where {T1, T2}
    Combined{Any, <: Tuple{T1, T2}}
end
# creates a new scene
plot(args...; kw_args...) = plot!(Scene(), Any, args...; kw_args...)
plot(P::Type, args...; kw_args...) = plot!(Scene(), P, args...; kw_args...)

# creates a new childscene
plot(scene::SceneLike, args...; kw_args...) = plot!(Scene(scene), Any, args...; kw_args...)
plot(scene::SceneLike, P::Type, args...; kw_args...) = plot!(Scene(scene), P, args...; kw_args...)

# plots to global current scene
plot!(args...; kw_args...) = plot!(current_scene(), Any, args...; kw_args...)
plot!(P::Type, args...; kw_args...) = plot!(current_scene(), P, Attributes(kw_args), args...)
plot!(P::Type, attributes::Attributes, args...) = plot!(current_scene(), P, attributes, args...)

# plots to scene
plot!(scene::SceneLike, args...; kw_args...) = plot!(scene, Any, args...; kw_args...)
plot!(scene::SceneLike, P::Type, args...; kw_args...) = plot!(scene, P, Attributes(kw_args), args...)

# function plot!(scene::SceneLike, plot::Combined{F, Arg}) where {F, Arg}
#     plot!(plot, arguments(plot)...; attributes(plot)...)
# end
# function plot!(scene::SceneLike, plot::Atomic{F, Arg}) where {F, Arg}
#     plot!(plot, arguments(plot)...; attributes(plot)...)
# end

function plot!(scene::SceneLike, ::Type{Any}, attributes::Attributes, args...)
    PlotType = plottype(value.(args)...)
    plot!(scene, PlotType, attributes, args...)
end

plot!(p::Atomic) = p

function plot!(scene::SceneLike, ::Type{PlotType}, attributes::Attributes, args...) where PlotType
    # create "empty" plot type - empty meaning containing no plots, just attributes + arguments
    plot_object, non_plot_kwargs = PlotType(scene, attributes, args)
    # call user defined overload to fill the plot type
    plot!(plot_object)
    # call the assembly recipe, that also adds this to the scene
    # kw_args not consumed by PlotType will be passed forward to plot! as non_plot_kwargs
    plot!(scene, plot_object, non_plot_kwargs)
end

function plot!(scene::Combined, ::Type{PlotType}, attributes::Attributes, args...) where PlotType
    # create "empty" plot type - empty meaning containing no plots, just attributes + arguments
    plot_object, non_plot_kwargs = PlotType(scene, attributes, args)
    # call user defined overload to fill the plot type
    plot!(plot_object)
    push!(scene.plots, plot_object)
    scene
end


text:
module Latexify
using Requires
using LaTeXStrings
using InteractiveUtils
using Markdown
using MacroTools: postwalk
import MacroTools
using Printf
using Formatting

export latexify, md, copy_to_clipboard, auto_display, set_default, get_default,
    reset_default, @latexrecipe, render, @latexify, @latexrun, @latexdefine

## Allow some backwards compatibility until its time to deprecate.
export latexequation, latexarray, latexalign, latexraw, latexinline, latextabular, mdtable

export StyledNumberFormatter, FancyNumberFormatter

COPY_TO_CLIPBOARD = false
function copy_to_clipboard(x::Bool)
    global COPY_TO_CLIPBOARD = x
end

AUTO_DISPLAY = false
function auto_display(x::Bool)
    global AUTO_DISPLAY = x
end

include("unicode2latex.jl")
include("function2latex.jl")
include("latexraw.jl")
include("latexoperation.jl")
include("latexarray.jl")
include("latexalign.jl")
include("latexbracket.jl")
include("latexinline.jl")
include("latexequation.jl")
include("latextabular.jl")
include("default_kwargs.jl")
include("recipes.jl")
include("macros.jl")

include("mdtable.jl")
include("mdtext.jl")
include("md.jl")

include("utils.jl")

include("numberformatters.jl")

include("latexify_function.jl")

### Add support for additional packages without adding them as dependencies.
function __init__()
    @require DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" begin
        include("plugins/ParameterizedFunctions.jl")
    end
    @require DiffEqBiological = "eb300fae-53e8-50a0-950c-e21f52c2b7e0" begin
        include("plugins/DiffEqBiological.jl")
    end
    @require SymEngine = "123dc426-2d89-5057-bbad-38513e3affd8" begin
        include("plugins/SymEngine.jl")
    end
    @require DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" begin
        include("plugins/DataFrames.jl")
    end
end

macro generate_test(expr)
    return :(clipboard("@test $($(string(expr))) == replace(\nraw\"$($(esc(expr)))\", \"\\r\\n\"=>\"\\n\")\n"))
end

"""
    @append_latexify_test!(fname, expr)

Generate a Latexify test and append it to the file `fname`.

The expression `expr` should return a string when evaluated.

Example use:
```
Latexify.@append_latexify_test!("./tests/latexify_tests.jl", latexify(:(x/y)))
```

The macro returns the output of the expression and can often be rendered
for a visual check that the test itself is ok.
```
Latexify.@append_latexify_test!("./tests/latexify_tests.jl", latexify(:(x/y))) |> render
```
"""
macro append_latexify_test!(fname, expr)
    fname = esc(fname)
    return :(
    str = "@test $($(string(expr))) == replace(\nraw\"$($(esc(expr)))\", \"\\r\\n\"=>\"\\n\")\n\n";
    open($fname, "a") do f
        write(f,str)
    end;
    $(esc(expr))
    )
end

"""
    @append_test!(fname, expr)

Both execute and append code to a test file.

The code can be either a normal expression or a string.
Example use:
```
Latexify.@append_test A = [1 2; 3 4]
```

Useful for adding code that generates objects to be used in latexify tests.
"""
macro append_test!(fname, str)
    fname = esc(fname)
    returnobj = str isa String ? Meta.parse(str) : str
    printobj = str isa String ? str : string(MacroTools.striplines(str))
    return :(
    open($fname, "a") do f
        write(f, $(esc(printobj)))
        write(f, "\n\n")
    end;
    $(esc(returnobj))
    )
end

end


text:
PropertyCardType = CardType{:property}()
const PropertyCard = NastranCard{PropertyCardType}

type PSHELL <: PropertyCard id::Int64
    mat1_id::Int64
    thickness::Nullable{Float64}
    nonstructural_mass::Float64
end
function convert(::Type{PSHELL},card)
    id = card[2]::Int64
    mat1_id = card[3]::Int64
    thickness = card[4] == "" ? Void : card[4]::Float64
    nonstructural_mass = card[9] == "" ? 0.0 : card[9]::Float64
    PSHELL(id,mat1_id,thickness,nonstructural_mass)
end

type PROD <: PropertyCard
    id::Int64
    mat_id::Int64
    A::Float64
    J::Float64
    C::Float64
    nonstructural_mass::Float64
end
function convert(::Type{PROD},card)
    id = card[2]::Int64
    mat_id = card[3]::Int64
    A = card[4]::Float64
    J = card[5]::Float64
    C = card[6] == "" ? 0.0 : card[6]::Float64
    nonstructural_mass = card[7] == "" ? 0.0 : card[7]::Float64
    PROD(id,mat_id,A,J,C,nonstructural_mass)
end

type PBUSH <: PropertyCard
    id::Int64
end
function convert(::Type{PBUSH},card)
    id = card[2]::Int64
    PBUSH(id)
end

type PBEAM <: PropertyCard
    id::Int64
    mat_id::Int64
    nonstructural_mass::Float64
    nonstructural_inertiaA::Float64
    nonstructural_inertiaB::Float64
end

is_continuation(field) = field == "YES"
is_double_continuation(field) = field == "YESA" || field == "NO"
function convert(::Type{PBEAM},card)
    id = card[2]::Int64
    mat_id = card[3]::Int64
    nonstructural_mass = card[9] == "" ? 0.0 : card[9]::Float64
    nonstructural_inertiaA = 0.0
    nonstructural_inertiaB = 0.0
    i = 10
    while i < length(card)
        if is_continuation(card[i])
        elseif is_double_continuation(card[i])
            i += 8
        else
            if card[i+4] != ""
                nonstructural_inertiaA = card[i+4]::Float64
            end
            if card[i+5] != ""
                nonstructural_inertiaB = card[i+5]::Float64
            end
        end
        i += 8
    end
    PBEAM(id,mat_id,nonstructural_mass,nonstructural_inertiaA,nonstructural_inertiaB)

end

type PELAS <: PropertyCard
    id::Int64
end
function convert(::Type{PELAS},card)
    id = card[2]::Int64
    PELAS(id)
end

type PAERO1 <: PropertyCard
    id::Int64
end
function convert(::Type{PAERO1},card)
    id = card[2]::Int64
    PAERO1(id)
end

type PSOLID <: PropertyCard
    id::Int64
end
function convert(::Type{PSOLID},card)
    id = card[2]::Int64
    PSOLID(id)
end


text:
import accrue.cryptoerase.runtime.Condition;

class C {
  static final Condition{L} c = new Condition();
  void m() {
      int{H} x = 37;
      if (x > 0) {
	  c.set();
      }
  }

  public static void main(String[] args) {
    new C().m();
  }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
using PyPlot
Tmax = 5
npoints = 10
g(n) = range(1/Tmax^n,1., length=npoints)
h(n) = range(0, n*log10(Tmax), length=npoints)
j(n) = range(0, 1-exp(-n*(Tmax-1)), length=npoints)

##
figure()
for n in range(0, 5, length=11)
    #plot(g(n).^(-1/n), g(n), label="T^$(-n)", ".-")
    plot(1 .- 1/n*log.(1 .- j(n)), j(n), "x-", label="1-exp(-$n(T -1))")
end
n = 1
plot(10 .^(1/n*h(n)), h(n), label="log T", "v-k")
xlabel("T")
ylabel("equally spaced f(T)")
legend(bbox_to_anchor=(1.04,1), borderaxespad=0)
plt.tight_layout()


text:
using Statistics, StatsBase
using PyCall
xr = pyimport("xarray")
using Elefridge
using ZfpCompression
using JLD

path = "/Users/milan/cams"
allvars = ["go3","so2","aermr04","aermr05","aermr06","ch4","co","co2","no"]
# allvars = ["no2"]

for vari in allvars
    println("--------")
    println("Compressing $vari")

    filelist = filter(x->endswith(x,vari*".grib"),readdir(path))
    X = xr.open_dataarray(joinpath(path,filelist[1]),engine="cfgrib").data

    # volume and error norms (columns):
    # vol, mean, median, 90% and max  of decimal error
    # mean, median, 90% and max of abs error
    # total 9

    # compression methods (rows):
    # ZFP lossless, ZFP tol = 1e-7,1e-5,1e-11, precision=14,17,20
    # ZFP lossless with RoundNearest16,RoundNearest24
    # total 9

    # storing array
    EV = Array{Float32,2}(undef,9,9)

    ## LOSSLESS
    Xc = zfp_compress(X)
    EV[1,1] = sizeof(X)/sizeof(Xc)     # compression factor
    EV[1,2:end] .= 0

    println("Lossless, $(EV[1,1])")

    # ZFP tol
    for (i,tol) in enumerate([1e-7,1e-9,1e-11])
        Xc = zfp_compress(X,tol=tol)
        EV[1+i,1] = sizeof(X)/sizeof(Xc)     # compression factor

        X2 = similar(X)
        zfp_decompress!(X2,Xc,tol=tol)

        if any(X2 .< 0)
            err = fill(-1,size(X))
        else
            err = abs.(log10.(X ./ X2))         # decimal error
        end
        EV[1+i,2] = mean(err)
        EV[1+i,3] = median(err)
        EV[1+i,4] = percentile(vec(err),90)
        EV[1+i,5] = maximum(err)

        err = abs.(X.-X2)   # absolute error
        EV[1+i,6] = mean(err)
        EV[1+i,7] = median(err)
        EV[1+i,8] = percentile(vec(err),90)
        EV[1+i,9] = maximum(err)
        println("tol=$tol, cf=$(EV[1+i,1]), max dec=$(EV[1+i,5]), max abs=$(EV[1+i,9])")
    end

    # ZFP precision
    for (i,prec) in enumerate([14,17,20])
        Xc = zfp_compress(X,precision=prec)
        EV[4+i,1] = sizeof(X)/sizeof(Xc)     # compression factor

        X2 = similar(X)
        zfp_decompress!(X2,Xc,precision=prec)

        if any(X2 .< 0)
            err = fill(-1,size(X))
        else
            err = abs.(log10.(X ./ X2))         # decimal error
        end
        EV[4+i,2] = mean(err)
        EV[4+i,3] = median(err)
        EV[4+i,4] = percentile(vec(err),90)
        EV[4+i,5] = maximum(err)

        err = abs.(X.-X2)   # absolute error
        EV[4+i,6] = mean(err)
        EV[4+i,7] = median(err)
        EV[4+i,8] = percentile(vec(err),90)
        EV[4+i,9] = maximum(err)
        println("precision=$prec, cf=$(EV[4+i,1]), max dec=$(EV[4+i,5]), max abs=$(EV[4+i,9])")
    end

    # RoundNearest16/24 + zfp lossless
    for (i,bits) in enumerate([7,15])
        X2 = round(X,bits)

        Xc = zfp_compress(X2)
        EV[7+i,1] = sizeof(X)/sizeof(Xc)     # compression factor

        err = abs.(log10.(X ./ X2))         # decimal error
        EV[7+i,2] = mean(err)
        EV[7+i,3] = median(err)
        EV[7+i,4] = percentile(vec(err),90)
        EV[7+i,5] = maximum(err)

        err = abs.(X.-X2)   # absolute error
        EV[7+i,6] = mean(err)
        EV[7+i,7] = median(err)
        EV[7+i,8] = percentile(vec(err),90)
        EV[7+i,9] = maximum(err)
        println("RN=$bits, cf=$(EV[7+i,1]), max dec=$(EV[7+i,5]), max abs=$(EV[7+i,9])")
    end

    save(joinpath(path,"error/$(vari)_zfp.jld"),"EV",EV)
    println("$vari stored.")
end


text:
"""
    function updateRanking(;
        players::Vector{Player},
        ranking_scheme_to_apply::String,
        ranking_scheme_to_update::String,
        last_tournament_id::String,
        verbose::Bool
    )

This function is used to update the player's rankings after a tournament.

Arguments
=========
- `players`:
    an array of players in the simulation.
- `ranking_scheme_to_apply`:
    the ranking scheme to apply.
    options:
    1) "ATP"
    2) "Elo"
    3) "Dwz"
- `ranking_scheme_to_update`:
    the ranking_scheme to update for the player.
- `last_tournament_id`:
    the identifier for the last tournament that was played.
- `verbose`:
    whether to print the results.
"""

function updateRanking(;
    players::Vector{Player},
    ranking_scheme_to_apply::String,
    ranking_scheme_to_update::String,
    last_tournament_id::String,
    verbose::Bool
)
    if ranking_scheme_to_apply == "ATP"
        atp_ranking(players, ranking_scheme_to_update)
    elseif ranking_scheme_to_apply == "Elo"
        elo_ranking(players, ranking_scheme_to_update, last_tournament_id)
    elseif ranking_scheme_to_apply == "Dwz"
        dwz_ranking(players, ranking_scheme_to_update, last_tournament_id)
    else
        throw(DomainError(ranking_scheme, "unsupported ranking scheme"))
    end
    if verbose
        println("########")
        println("$ranking_scheme_to_update updated with $ranking_scheme_to_apply")
        println("########")
        println()
    end
    return nothing
end

max_stars_repo_name: Danial-Hussain/Tennis-Ranking

max_issues_repo_name: Danial-Hussain/Tennis-Ranking

max_forks_repo_name: Danial-Hussain/Tennis-Ranking