text:
using HiQGA.transD_GP, PyPlot
## model fixed parts, i.e., water
zfixed   = [0]
ρfixed   = [0.22]
nmax = 200
# z grid spec starts, the first z and first ρ will be unused in MT
zstart = 35.
extendfrac, dz = 1.028, 10
zall, znall, zboundaries = transD_GP.setupz(zstart, extendfrac, dz=dz, n=50, showplot=true)
z, ρ, nfixed = transD_GP.makezρ(zboundaries; zfixed=zfixed, ρfixed=ρfixed)
## fill in detail in ohm-m
ρ[(z.>=zstart) .& (z.<100)] .= 4.
ρ[(z.>=100) .& (z.<250)] .= 14
ρ[(z.>=250) .&(z.<600)] .= 0.8
ρ[z.>=600] .= 100
##
T = 10 .^range(-2, 0, length=10)
F = transD_GP.MT1DInversion.create_synthetic(ρ =  ρ[2:end], zboundaries = zboundaries,
                                            freqs = 1 ./T, rseed=125, showplot=true, noisefrac=0.05, logscaledepth=false, irxlayer=1)

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA

max_forks_repo_name: GeoscienceAustralia/HiQGA


text:
module StructsOfArrays

export StructOfArrays, ScalarRepeat

struct StructOfArrays{T,N,U<:Tuple} <: AbstractArray{T,N}
    arrays::U
end

mutable struct ScalarRepeat{T}
    scalar::T
end
Base.ndims(::ScalarRepeat) = 1
Base.getindex(s::ScalarRepeat, i...) = s.scalar
#should setindex! really be allowed? It will set the index for the whole row...
Base.setindex!(s::ScalarRepeat{T}, value, i...) where {T} = (s.scalar = T(value))
Base.eltype(::ScalarRepeat{T}) where {T} = T

Base.start(::ScalarRepeat) = 1
Base.next(sr::ScalarRepeat, i) = sr.scalar, i+1
Base.done(sr::ScalarRepeat, i) = false

# since this is used in hot loops, and T.types[.] doesn't play well with compiler
# this needs to be a generated function
@generated function is_tuple_struct(::Type{T}) where T
    is_ts = length(T.types) == 1 && T.types[1] <: Tuple
    :($is_ts)
end
struct_eltypes(::T) where {T} = struct_eltypes(T)
function struct_eltypes(::Type{T}) where T
    if is_tuple_struct(T) #special case tuple types (E.g. FixedSizeVectors)
        return eltypes = T.types[1].parameters
    else
        return eltypes = T.types
    end
end

make_iterable(x::AbstractArray) = x
make_iterable(x) = ScalarRepeat(x)

@generated function StructOfArrays(::Type{T}, dim1::Integer, rest::Integer...) where T
    (!isleaftype(T) || T.mutable) && return :(throw(ArgumentError("can only create an StructOfArrays of leaf type immutables")))
    isempty(T.types) && return :(throw(ArgumentError("cannot create an StructOfArrays of an empty or bitstype")))
    dims = (dim1, rest...)
    N = length(dims)
    eltypes  = struct_eltypes(T)
    arrtuple = Tuple{[Array{eltypes[i],N} for i = 1:length(eltypes)]...}

    :(StructOfArrays{T,$N,$arrtuple}(
        ($([:(Array($(eltypes[i]), (dim1, rest...))) for i = 1:length(eltypes)]...),)
    ))
end
StructOfArrays(T::Type, dims::Tuple{Vararg{Integer}}) = StructOfArrays(T, dims...)

function StructOfArrays(T::Type, a, rest...)
    arrays = map(make_iterable, (a, rest...))
    N = ndims(arrays[1])
    eltypes = map(eltype, arrays)
    s_eltypes = struct_eltypes(T)
    any(ix->ix[1]!=ix[2], zip(eltypes,s_eltypes)) && throw(ArgumentError(
        "fieldtypes of $T must be equal to eltypes of arrays: $eltypes"
    ))
    any(x->ndims(x)!=N, arrays) && throw(ArgumentError(
        "cannot create an StructOfArrays from arrays with different ndims"
    ))
    arrtuple = Tuple{map(typeof, arrays)...}
    StructOfArrays{T, N, arrtuple}(arrays)
end

Base.IndexStyle(::Type{<:StructOfArrays}) = IndexLinear()

@generated function Base.similar(A::StructOfArrays, ::Type{T}, dims::Dims) where T
    if isbits(T) && length(T.types) > 1
        :(StructOfArrays(T, dims))
    else
        :(Array(T, dims))
    end
end

Base.convert(::Type{StructOfArrays{T,N}}, A::AbstractArray{S,N}) where {T,S,N} =
    copy!(StructOfArrays(T, size(A)), A)
Base.convert(::Type{StructOfArrays{T}}, A::AbstractArray{S,N}) where {T,S,N} =
    convert(StructOfArrays{T,N}, A)
Base.convert(::Type{StructOfArrays}, A::AbstractArray{T,N}) where {T,N} =
    convert(StructOfArrays{T,N}, A)

function Base.size(A::StructOfArrays{T, N, U}) where {T, N, U}
    for elem in A.arrays
        if isa(elem, AbstractArray)
            return size(elem)
        end
    end
    # if none has a size, size is inf!
    ntuple(Val{N}) do i
        typemax(Int)
    end
end

@generated function Base.getindex(A::StructOfArrays{T}, i::Integer...) where T
    n = length(struct_eltypes(T))
    Expr(:block, Expr(:meta, :inline),
         :($T($([:(A.arrays[$j][i...]) for j = 1:n]...)))
    )
end

function _getindex(x::T, i) where T
    is_tuple_struct(T) ? x[i] : getfield(x, i)
end
@generated function Base.setindex!(A::StructOfArrays{T}, x, i::Integer...) where T
    n = length(struct_eltypes(T))
    quote
        $(Expr(:meta, :inline))
        v = convert(T, x)
        $([:(A.arrays[$j][i...] = _getindex(v, $j)) for j = 1:n]...)
        x
    end
end
end # module

ext: jl


text:
### Bad data processing
function baddata(settings, numsys, x, z, v, W, G, J, Nmeasure, idx, badsave, rbelow)
    ########## Sparse inverse ##########
    S, L, U, d, p, q, Rs = constructfact(G)
    parent = etree(S)
    R = symbfact(S, parent)
    Gi = sparseinv(L, U, d, p, q, Rs, R)

    ########## Diagonal entries of residual matrix ##########
    JGi = J * Gi
    c = fill(0.0, size(J, 1))
    for i in idx
        c[i[1]] += JGi[i] * J[i]
    end

    ########## Critical data ##########
    xcr = (J' * J) \ (J' * z)
    hcr = J * xcr

    ########## Largest normalized residual ##########
    h = J * x
    rmax = 0.0; idxr = 0
    @inbounds for i = 1:Nmeasure
        if abs(z[i] - hcr[i]) > settings.bad[:critical]
            rnor = abs(z[i] - h[i]) / sqrt(abs(v[i] - c[i]))
            if rnor > rmax
                idxr = i
                rmax = rnor
            end
        end
    end

    if rmax < settings.bad[:treshold]
        rbelow = false
    end

    if idxr == 0
        println("All measurements are marked as critical, the bad data processing is terminated.")
        rbelow = false
    end

    ########## Remove measurement ##########
    if rbelow
        @inbounds for i in idx
            if i[1] == idxr
                J[i] = 0.0
            end
        end
        z[idxr] = 0.0
        @inbounds for i in W.rowval
            if i == idxr
                W.nzval[i] = 0.0
            end
        end
        badsave = [badsave; [idxr rmax]]
    end

    return J, z, W, badsave, rbelow
end # baddata

### Observability restoration based on the Gram matrix associated to a reduced network and measurement set
function restorationgram(system, settings, numsys, num, measurements, islands, Nisland, J, mean, weight,
    Ybus, slack, Tslack, branch, from, to, fromPij, toPij, busPi, branchPij, busTi, onPi, onPij, onTi, Nvol, Nflow,
    meanPij, meanPi, meanTi, variPi, variTi, variPij, admitance, Gshunt, transShift, Pshift)
    ######### Bus-islands ##########
    islandsWhere = fill(0, numsys.Nbus)
    @inbounds for (k, island) in enumerate(islands)
        for i in island
            islandsWhere[i] = k
        end
    end

    ########## Tie buses and branches ##########
    tie_branch = copy(branch)
    @inbounds for i = 1:numsys.Nbranch
        if islandsWhere[from[i]] == islandsWhere[to[i]]
            tie_branch[i] = 0
        end
    end

    tie_bus = fill(0, numsys.Nbus)
    @inbounds for i in tie_branch
        if i != 0
            tie_bus[from[i]] = from[i]
            tie_bus[to[i]] = to[i]
        end
    end

    ########## Reduced Jacobian Wb ##########
    Nreduced = Nvol
    @inbounds for i = 1:num.legacyNi
        if onPi[i] == 1 && tie_bus[busPi[i]] != 0
            Nreduced += 1
        end
    end
    idx = fill(0, Nreduced)
    N = Nflow; cnt = 1
    @inbounds for i = 1:num.legacyNi
        if onPi[i] == 1
            N += 1
            if tie_bus[busPi[i]] != 0
                idx[cnt] = N
                cnt += 1
            end
        end
    end
    @inbounds for i = 1:num.pmuNv
        if onTi[i] == 1
            N += 1
            idx[cnt] = N
            cnt += 1
        end
    end
    Jt = copy(transpose(J))
    col = fill(0, numsys.Nbus)
    @inbounds for (k, island) in enumerate(islands)
        for j in island
            col[j] = k
        end
    end
    N = 0
    @inbounds for (k, i) in enumerate(idx)
        for j in Jt.colptr[i]:(Jt.colptr[i + 1] - 1)
            N += 1
        end
    end
    rown = fill(0, N); coln = fill(0, N); valn = fill(0.0, N)
    cnt = 1
    @inbounds for (k, i) in enumerate(idx)
        for j in Jt.colptr[i]:(Jt.colptr[i + 1] - 1)
            rown[cnt] = k
            coln[cnt] = col[Jt.rowval[j]]
            valn[cnt] = Jt.nzval[j]
            cnt += 1
        end
    end
    Wb = sparse(rown, coln, valn, length(idx), Nisland)

    ########## Pseudo-measurements Jacobian Wp ##########
    Nelement = 0; Nmeasuren = 0
    Nflown = 0; Ninjn = 0; Nvoln = 0;
    if settings.observe[:Pij] != 0
        @inbounds for i = 1:num.legacyNf
            if onPij[i] == 0 && tie_branch[branchPij[i]] != 0
               Nelement += 2; Nmeasuren += 1; Nflown += 1
            end
        end
    end
    if settings.observe[:Pi] != 0
        @inbounds for i = 1:num.legacyNi
            if onPi[i] == 0 && tie_bus[busPi[i]] != 0
                Nmeasuren += 1; Ninjn += 1
                for j in Ybus.colptr[busPi[i]]:(Ybus.colptr[busPi[i] + 1] - 1)
                    Nelement += 1
                end
            end
        end
    end
    if settings.observe[:Ti] != 0
        @inbounds for i = 1:num.pmuNv
            if onTi[i] == 0 && tie_bus[busTi[i]] != 0
                Nelement += 1; Nmeasuren += 1; Nvoln += 1
            end
        end
    end

    rown = fill(0, Nelement); coln = similar(rown); jacn = fill(0.0, Nelement)
    meann = fill(0.0, Nmeasuren); weightn = similar(meann)
    idx = fill(0, Nmeasuren); index = 1; rowindex = 1
    Xflown = fill(0, Nflown, 2); Xinjn = fill(0, Ninjn, 2); Xvoln = fill(0, Nvoln, 2)
    if settings.observe[:Pij] != 0
        @inbounds for i = 1:num.legacyNf
            if onPij[i] == 0 && tie_branch[branchPij[i]] != 0
                rown[index] = rowindex
                coln[index] = fromPij[i]
                jacn[index] = admitance[branchPij[i]]
                index += 1
                rown[index] = rowindex
                coln[index] = toPij[i]
                jacn[index] = -admitance[branchPij[i]]

                if fromPij[i] == from[branchPij[i]]
                    meann[rowindex] = meanPij[i] + transShift[branchPij[i]] * (pi / 180)  * admitance[branchPij[i]]
                else
                    meann[rowindex] = meanPij[i] - transShift[branchPij[i]] * (pi / 180) * admitance[branchPij[i]]
                end

                weightn[rowindex] = 1 / settings.observe[:Pij]
                idx[rowindex] = rowindex

                Xflown[rowindex, 1] = i; Xflown[rowindex, 2] = rowindex
                index += 1; rowindex += 1
            end
        end
    end

    cnt = 1
    if settings.observe[:Pi] != 0
        @inbounds for i = 1:num.legacyNi
            if onPi[i] == 0 && tie_bus[busPi[i]] != 0
                for j in Ybus.colptr[busPi[i]]:(Ybus.colptr[busPi[i] + 1] - 1)
                    rown[index] = rowindex
                    coln[index] = Ybus.rowval[j]
                    jacn[index] = Ybus.nzval[j]
                    index += 1
                end

                meann[rowindex] = meanPi[i] - Pshift[i] - Gshunt[i] / system.basePower
                weightn[rowindex] = 1 / settings.observe[:Pi]
                idx[rowindex] = rowindex
                Xinjn[cnt, 1] = i; Xinjn[cnt, 2] = rowindex
                rowindex += 1; cnt += 1
            end
        end
    end

    cnt = 1
    if settings.observe[:Ti] != 0
        @inbounds for i = 1:num.pmuNv
            if onTi[i] == 0 && tie_bus[busTi[i]] != 0
                rown[index] = rowindex
                coln[index] = busTi[i]
                jacn[index] = 1.0

                meann[rowindex] = meanTi[i] - Tslack * (pi / 180)
                weightn[rowindex] = 1 / settings.observe[:Ti]
                idx[rowindex] = rowindex

                Xvoln[cnt, 1] = i; Xvoln[cnt, 2] = rowindex
                index += 1; rowindex += 1; cnt += 1
            end
        end
    end
    JpsT = sparse(coln, rown, jacn, numsys.Nbus, Nmeasuren)

    ########## Reduced Jacobian Wc ##########
    N = 0
    @inbounds for (k, i) in enumerate(idx)
        for j in JpsT.colptr[i]:(JpsT.colptr[i + 1] - 1)
            N += 1
        end
    end
    rown = fill(0, N); coln = fill(0, N); valn = fill(0.0, N)
    cnt = 1
    @inbounds for (k, i) in enumerate(idx)
        for j in JpsT.colptr[i]:(JpsT.colptr[i + 1] - 1)
            rown[cnt] = k
            coln[cnt] = col[JpsT.rowval[j]]
            valn[cnt] = JpsT.nzval[j]
            cnt += 1
        end
    end
    Wp = sparse(rown, coln, valn, length(idx), Nisland)

    ########## Restore observability ##########
    M = [Wb; Wp]
    M = M * M'
    F = qr(Matrix(M))
    R = F.R

    pseudo = Int64[]
    @inbounds for i = (Nreduced + 1):(Nreduced + Nmeasuren)
        if abs(R[i, i]) > settings.observe[:pivot]
            push!(pseudo, i - Nreduced)
        end
    end
    @inbounds for i in pseudo
        for j in Xflown[:, 2]
            if i == j
                onPij[Xflown[j, 1]] = 2
                variPij[Xflown[j, 1]] = settings.observe[:Pij]
            end
        end
        for j in Xinjn[:, 2]
            if i == j
                onPi[Xinjn[j, 1]] = 2
                variPi[Xinjn[j, 1]] = settings.observe[:Pi]
            end
        end
        for j in Xvoln[:, 2]
            if i == j
                onTi[Xvoln[j, 1]] = 2
                variTi[Xvoln[j, 1]] = settings.observe[:Ti]
            end
        end
    end

    Npseudo = length(pseudo)
    J = [J; transpose(JpsT[:, pseudo])]
    mean = [mean; meann[pseudo]]
    weight = [weight; weightn[pseudo]]

    return J, mean, weight, Npseudo
end # restorationgram

### Gaussian belief propagation observability restoration
function restorationbp(system, settings, numsys, num, measurements, islands, Nisland, J, mean, weight,
    Ybus, branch, from, to, busPi, onPi, meanPi, variPi, Pshift, Gshunt)
    ######### Bus-islands ##########
    islandsWhere = fill(0, numsys.Nbus)
    @inbounds for (k, island) in enumerate(islands)
        for i in island
            islandsWhere[i] = k
        end
    end

    ######### Tie buses and branches ##########
    tie_branch = copy(branch)
    @inbounds for i = 1:numsys.Nbranch
        if islandsWhere[from[i]] == islandsWhere[to[i]]
            tie_branch[i] = 0
        end
    end
    filter!(x->x!=0, tie_branch)

    tie_bus = fill(0, numsys.Nbus)
    @inbounds for i in tie_branch
        tie_bus[from[i]] = from[i]
        tie_bus[to[i]] = to[i]
    end

    ######### Tie real-time measurements and pseudo-measurements ##########
    Ntie = 0; Neletie = 0
    @inbounds for i = 1:num.legacyNi
        if tie_bus[busPi[i]] != 0
            Ntie += 1
            for j in Ybus.colptr[busPi[i]]:(Ybus.colptr[busPi[i] + 1] - 1)
                Neletie += 1
            end
        end
    end

    rowtie = fill(0, Neletie); coltie = similar(rowtie); jactie = fill(0.0, Neletie); idx = fill(0, Ntie);
    meantie = fill(0.0, Ntie); weighttie = similar(meantie)
    injbus = fill(0, Ntie); incidence = fill(false, numsys.Nbus); pseudo = Int64[]
    Nreal = 0; index = 1; rowindex = 1
    @inbounds for i = 1:num.legacyNi
        if tie_bus[busPi[i]] != 0
            for j in Ybus.colptr[busPi[i]]:(Ybus.colptr[busPi[i] + 1] - 1)
                rowtie[index] = rowindex
                coltie[index] = Ybus.rowval[j]
                jactie[index] = Ybus.nzval[j]
                incidence[Ybus.rowval[j]] = true
                index += 1
            end
            meantie[rowindex] = meanPi[i] - Pshift[i] - Gshunt[i] / system.basePower
            weighttie[rowindex] = 1 / variPi[i]
            idx[rowindex] = rowindex
            injbus[rowindex] = tie_bus[busPi[i]]

            if onPi[i] == 0
                push!(pseudo, rowindex)
            else
                Nreal += 1
            end
            rowindex += 1
        end
    end
    JtieT = sparse(coltie, rowtie, jactie, numsys.Nbus, Ntie)

    ######### Reduce matrix ##########
    N = 0
    @inbounds for (k, i) in enumerate(idx)
        for j in JtieT.colptr[i]:(JtieT.colptr[i + 1] - 1)
            N += 1
        end
    end

    col = fill(0, numsys.Nbus)
    @inbounds for (k, island) in enumerate(islands)
        for j in island
            col[j] = k
        end
    end

    rown = fill(0, N); coln = fill(0, N); valn = fill(0.0, N); cnt = 1
    @inbounds for (k, i) in enumerate(idx)
        for j in JtieT.colptr[i]:(JtieT.colptr[i + 1] - 1)
            rown[cnt] = k
            coln[cnt] = col[JtieT.rowval[j]]
            valn[cnt] = JtieT.nzval[j]
            cnt += 1
        end
    end
    Mt = sparse(coln, rown, valn, Nisland, length(idx))

    ########### Factor graph data ##########
    Nvar, Nfac = size(Mt); Nlink = nnz(Mt)
    Vvar_fac, Wfac_var, to_fac, to_var, colptr, rowptr, idxT, row, col, links = factorgraph(Mt, Nvar, Nfac, Nlink, Nfac)
    Wdir = fill(1e-120, Nvar); Vind = fill(1.0, Nfac); Vind[pseudo] .= 1e120

    ######### Pass messages from singly-connected factor nodes to all indirect links ##########
    @inbounds for i = 1:Nvar
        for j in colptr[i]:(colptr[i + 1] - 1)
            k = to_fac[j]
            Vvar_fac[k] = 1 / Wdir[i]
        end
    end

     ########## Inference ##########
    need = Nisland - 1 - Nreal; pass = 1; bppseudo = Int64[]; residual = 0.0; pre = 0.0; pseudom = Int64[]
    while need > 0
        on = pseudo[pass]
        Wdir[links[on][1]] = 1e120
        Vind[on] = 1

        for iter = 1:settings.observe[:restoreMax]
            @inbounds for i = 1:Nfac
                Vrow = 0.0
                for j in rowptr[i]:(rowptr[i + 1] - 1)
                    Vrow += Vvar_fac[j]
                end
                for j in rowptr[i]:(rowptr[i + 1] - 1)
                    Wfac_var[to_var[j]] = (1) / (Vind[i] + (Vrow -  Vvar_fac[j]))
                end
            end

            @inbounds for i = 1:Nvar
                Wcol = 0.0
                for j in colptr[i]:(colptr[i + 1] - 1)
                    Wcol += Wfac_var[j]
                end
                for j in colptr[i]:(colptr[i + 1] - 1)
                    Vvar_fac[to_fac[j]] = 1 / (Wdir[i] + (Wcol  - Wfac_var[j]))
                end
            end
        end

        residual = Vind[on]
        @inbounds for j in rowptr[on]:(rowptr[on + 1] - 1)
            residual += Vvar_fac[j]
        end

        Wdir[links[on][1]] = 1e-120
        if residual > 1e2
            need -= 1
            push!(bppseudo, injbus[on])
            push!(pseudom, on)
        else
            Vind[on] = 1e120
        end
        pass += 1
    end

    onPi[bppseudo] .= 2
    variPi[bppseudo] .= settings.observe[:Pi]

    J = [J; transpose(JtieT[:, pseudom])]
    mean = [mean; meantie[pseudom]]
    weight = [weight; weighttie[pseudom]]

    return J, mean, weight, Nisland - 1 - Nreal
end # restorationbp

### Flow islands defined by active power flow measurements
function islandsflow(Jflow)
    M = SimpleWeightedGraph(Jflow' * Jflow)
    islandsFlow = connected_components(M)
    NislandsFlow = size(islandsFlow, 1)

    return islandsFlow, NislandsFlow
end # islandsflow

### Gaussian belief propagation island detection
function islandsbp(settings, numsys, J, bus, islandsFlow, NislandsFlow)
    ########## Independent system of equations ##########
    Jt = copy(J')
    Ht = independentset(Jt)
    bus_indecies = copy(bus)

    ########## Factor graph data ##########
    Nvar, Nfac = size(Ht); Nlink = 0; Nind = 0
    Wvirtual = 0.0; Wdirect = 1e240
    Wdir = fill(Wvirtual, Nvar)

    ########## Find islands with isolated buses ##########
    variable_isolate = fill(0, Nvar)
    @inbounds for i = 1:Nfac
        for j in Ht.colptr[i]:(Ht.colptr[i + 1] - 1)
            variable_isolate[Ht.rowval[j]] += 1
        end
        var_in_column = Ht.colptr[i + 1] - Ht.colptr[i]
        if var_in_column == 1
            variable_isolate[Ht.rowval[Ht.colptr[i]]] -= 1
        else
            Nlink += var_in_column
            Nind += 1
        end
    end

    ########## Set probe nodes ##########
    islandsWhere = fill(0, numsys.Nbus)
    numvar = fill(0, NislandsFlow)
    @inbounds for (k, island) in enumerate(islandsFlow)
        for i in island
            islandsWhere[i] = k
            numvar[k] += 1
        end
    end
    numvar_max = maximum(numvar)

    if numvar_max > 1
        setprobe = findfirst(a -> a == numvar_max, numvar)
        @inbounds for i in islandsFlow[setprobe]
            Wdir[i] = Wdirect
        end
    else
        @inbounds for i = 1:Nvar
            if variable_isolate[i] != 0
                Wdir[i] = Wdirect
                break
            end
        end
    end

    ########## Find Islands ##########
    detection = true; previous = 0.0; islands = [Int64[]]
    iter_max = settings.observe[:islandMax]; breaks = settings.observe[:islandBreak]
    stop = settings.observe[:islandStop]; treshold = settings.observe[:islandTreshold];
    while detection
        ########## Factor graph data ##########
        Vvar_fac, Wfac_var, to_fac, to_var, colptr, rowptr, idxT, row, col = factorgraph(Ht, Nvar, Nfac, Nlink, Nind)

        ########## Inference ##########
        Xbp = fill(0.0, Nvar)
        @inbounds for iter = 1:iter_max
            for i = 1:Nind
                Vrow = 0.0
                for j in rowptr[i]:(rowptr[i + 1] - 1)
                    Vrow += Vvar_fac[j]
                end
                for j in rowptr[i]:(rowptr[i + 1] - 1)
                    Wfac_var[to_var[j]] = 1 / (Vrow - Vvar_fac[j] + eps(Float64))
                end
            end

            previous = 0.0
            for i = 1:Nvar
                Wcol = 0.0
                for j in colptr[i]:(colptr[i + 1] - 1)
                    Wcol += Wfac_var[j]
                end
                for j in colptr[i]:(colptr[i + 1] - 1)
                    Vvar_fac[to_fac[j]] = 1 / (Wcol + Wdir[i] - Wfac_var[j] + eps(Float64))
                end
                temp = 1 / (Wcol + Wdir[i])
                if temp < 10e30
                    previous = maximum([previous; abs(temp - Xbp[i])])
                end
                Xbp[i] = temp

            end
            if previous < stop && iter > breaks
                break
            end
        end

        ########## Reduce Jacobian ##########
        row = fill(true, Nvar); col = fill(true, Nfac); bip = fill(false, Nvar)
        @inbounds for i in idxT
            if Xbp[i[1]] < treshold
                row[i[1]] = false
                col[i[2]] = false
                bip[i[1]] = true
            end
        end
        Ht = Ht[row, col]
        Nvar, Nfac = size(Ht); Nlink = 0; Nind = 0

        ########## Find island ##########
        observable = bus_indecies[bip]
        unobservable = bus_indecies[.!bip]
        islands = [islands; [observable]]

        deleteat!(bus_indecies, bip)
        if length(bus_indecies) == 0
            break
        end

        ########## Find islands with isolated buses ##########
        variable_isolate = fill(0, Nvar)
        @inbounds for i = 1:Nfac
            for j in Ht.colptr[i]:(Ht.colptr[i + 1] - 1)
                variable_isolate[Ht.rowval[j]] += 1
            end
            var_in_column = Ht.colptr[i + 1] - Ht.colptr[i]
            if var_in_column == 1
                variable_isolate[Ht.rowval[Ht.colptr[i]]] -= 1
            else
                Nlink += var_in_column
                Nind += 1
            end
        end

        ########## Set probe nodes ##########
        unobservable_island = islandsWhere[unobservable]
        setprobe_idx = findfirst(a -> a != 0, variable_isolate)
        if setprobe_idx == nothing
            @inbounds for i in bus_indecies
                islands = [islands; [[i]]]
            end
            break
        end
        setprobe = unobservable_island[setprobe_idx]
        Wdir = fill(Wvirtual, Nvar)
        @inbounds for i = 1:Nvar
            if unobservable_island[i] == setprobe
                Wdir[i] = Wdirect
            end
        end
    end
    deleteat!(islands, 1)
    Nisland = size(islands, 1)

    return islands, Nisland
end # islandsbp

### Topological island detection
function islandstopological(settings, numsys, J, Jflow, Ybus, branch, from, to, busPi, onPi, islands)
    islandsWhere = fill(0, numsys.Nbus)
    @inbounds for (k, island) in enumerate(islands)
        for i in island
            islandsWhere[i] = k
        end
    end

    ########## Tie buses and branches ##########
    tie_branch = copy(branch)
    @inbounds for i = 1:numsys.Nbranch
        if islandsWhere[from[i]] == islandsWhere[to[i]]
            tie_branch[i] = 0
        end
    end
    tie_bus = fill(0, numsys.Nbus)
    @inbounds for i in tie_branch
        if i != 0
            tie_bus[from[i]] = from[i]
            tie_bus[to[i]] = to[i]
        end
    end

    ########## Boundary Injections ##########
    inj = copy(busPi); inj = inj[tie_bus[busPi] .!= 0]
    @inbounds for (k, i) in enumerate(inj)
        if onPi[i] == 0
            inj[k] = 0
        end
    end
    filter!(x->x!=0, inj); filter!(x->x!=0, tie_branch)

    ########## Stage 1: Processing of Individual Boundary Injections ##########
    islands, islandsWhere, inj, tie_branch = topological_merge_pairs(numsys, Ybus, busPi, islands, inj, islandsWhere, tie_branch, from, to)

    ########## Stage K: Processing of a Collection of K Boundary Injection Measurements ##########
    Ninj = length(inj); merge = 1; con = fill(false, numsys.Nbus); island_inj = fill(Int[], Ninj, 1)
    while merge != 0
        for (k, i) in enumerate(inj)
            conection = Ybus.colptr[busPi[i]]:(Ybus.colptr[busPi[i] + 1] - 1)
            con[Ybus.rowval[conection]] .= true
            island_inj[k] = sort(unique(islandsWhere[con]))
            con[Ybus.rowval[conection]] .= false
        end
        merge_index = decision_tree(island_inj)

        if merge_index != false
            merge_islands = []
            for i in merge_index
                for j in island_inj[i]
                    push!(merge_islands, j)
                end
            end
            merge_islands = unique(merge_islands)
            start = merge_islands[1]
            for i = 2:length(merge_islands)
                next = merge_islands[i]
                islands[start] = [islands[start]; islands[next]]
                islands[next] = []
            end
            empty = isempty.(islands)
            deleteat!(islands, empty)
        else
            break
        end

        islandsWhere = fill(0, numsys.Nbus)
        @inbounds for (k, island) in enumerate(islands)
            for i in island
                islandsWhere[i] = k
            end
        end
        @inbounds for (k, i) in enumerate(tie_branch)
            if islandsWhere[from[i]] == islandsWhere[to[i]]
                tie_branch[k] = 0
            end
        end
        tie_bus = fill(0, numsys.Nbus)
        @inbounds for i in tie_branch
            if i != 0
                tie_bus[from[i]] = from[i]
                tie_bus[to[i]] = to[i]
            end
        end
        inj = tie_bus[inj]
        filter!(x->x!=0, inj); filter!(x->x!=0, tie_branch)

        islands, islandsWhere, inj, tie_branch = topological_merge_pairs(numsys, Ybus, busPi, islands, inj, islandsWhere, tie_branch, from, to)

        Ninj = length(inj)
        con = fill(false, numsys.Nbus)
        island_inj = fill(Int[], Ninj, 1)
    end
    empty = isempty.(islands)
    deleteat!(islands, empty)
    Nisland = size(islands, 1)

    return islands, Nisland
end # islandstopological

function topological_merge_pairs(numsys, Ybus, busPi, islands, inj, islandsWhere, tie_branch, from, to)
    merge = 1; current = 0
    con = fill(false, numsys.Nbus)
    flag = false
    @inbounds while merge != 0
        merge = 1
        for (k, i) in enumerate(inj)
            current = islandsWhere[i]
            conection = Ybus.colptr[busPi[i]]:(Ybus.colptr[busPi[i] + 1] - 1)
            con[Ybus.rowval[conection]] .= true
            con[islands[current]] .= false
            toOther = islandsWhere[con]
            con[Ybus.rowval[conection]] .= false
            if !(isempty(toOther)) && all(y -> y == toOther[1], toOther)
                merge = 2
                islands[current] = [islands[current]; islands[toOther[1]]]
                islandsWhere[islands[toOther[1]]] .= current
                islands[toOther[1]] = []
                deleteat!(inj, k)
                flag = true
            end
        end
        if merge == 1 || size(islands, 1) == 1
            merge = 0
        end
    end
    empty = isempty.(islands)
    deleteat!(islands, empty)

    if flag
        islandsWhere = fill(0, numsys.Nbus)
        @inbounds for (k, island) in enumerate(islands)
            for i in island
                islandsWhere[i] = k
            end
        end

        @inbounds for (k, i) in enumerate(tie_branch)
            if islandsWhere[from[i]] == islandsWhere[to[i]]
                    tie_branch[k] = 0
            end
        end
        tie_bus = fill(0, numsys.Nbus)
        @inbounds for i in tie_branch
            if i != 0
                tie_bus[from[i]] = from[i]
                tie_bus[to[i]] = to[i]
            end
        end
        inj = tie_bus[inj]
        filter!(x->x!=0, inj)
        filter!(x->x!=0, tie_branch)
    end

    return islands, islandsWhere, inj, tie_branch
end

function combinations_recursive(position, value, result, max_n, k, accumulator)
    for i = value:max_n
        result[position] = position + i
        if position < k
            combinations_recursive(position + 1, i, result, max_n, k, accumulator)
        else
            push!(accumulator, copy(result))
        end
    end
end

function combinations(n, k)
    max_n = n - k
    accumulator = []
    result = zeros(Int, k)
    combinations_recursive(1, 0, result, max_n, k, accumulator)

    return accumulator
end

function check(measurments, indicies, total, required)
    appeared = zeros(Bool, total)
    for index in indicies
        for island in measurments[index]
            appeared[island] = true
        end
    end

    return sum(appeared) == required
end

function decision_tree(measurments)
    total_islands = 0
    for measurment in measurments
        for island in measurment
            total_islands = max(total_islands, island)
        end
    end

    for t = 2:length(measurments)
        total_combinations = combinations(length(measurments), t)
        for combination in total_combinations
            if check(measurments, combination, total_islands , t + 1)
                return combination
            end
        end
    end

    return false
end

### Observability analysis using the nodal variable formulation without the restore routine (on hold util Julia 1.6 releases)
function islandsnodal(numsys, num, settings, J, bus, branch, onPij, branchPij, onPi, busPi, from, to, onTi, busTi)
    ########## Find irrelevant branches ##########
    branch_non = fill(0, numsys.Nbranch)

    @inbounds for i = 1:num.legacyNf
        if onPij[i] == 1
            branch_non[branchPij[i]] = 1
        end
    end
    @inbounds for i = 1:num.legacyNi
        if onPi[i] == 1
            for j = 1:numsys.Nbranch
                if busPi[i] == from[j] || busPi[i] == to[j]
                    branch_non[j] = 1
                end
            end
        end
    end
    @inbounds for i = 1:num.pmuNv
        if onTi[i] == 1
            for j = 1:numsys.Nbranch
                if busTi[i] == from[j] || busTi[i] == to[j]
                    branch_non[j] = 1
                end
            end
        end
    end
    branch_irrelevant = copy(branch_non)

    ########## Branch-bus incidence matrix ##########
    row = collect(1:numsys.Nbranch)
    A = sparse([row; row], [from; to], [branch_non; -branch_non], numsys.Nbranch, numsys.Nbus)

    ########## Initialize position for angle pseudo-measurements ##########
    pseudo = sparse(bus, bus, zeros(numsys.Nbus), numsys.Nbus, numsys.Nbus)
    H = [J; pseudo]
    Ht = copy(transpose(H))
    rowsdel = collect(1:Nmeasure)

    ########## Branches and incident measurements ##########
    branch_measure = [Int64[] for i in 1:numsys.Nbranch]
    cnt = 0
    @inbounds for i = 1:num.legacyNf
        if onPij[i] == 1
            cnt += 1
            for j = 1:numsys.Nbranch
                if branchPij[i] == branch[j]
                    push!(branch_measure[j], cnt)
                end
            end
        end
    end
    @inbounds for i = 1:num.legacyNi
        if onPi[i] == 1
            cnt += 1
            for j = 1:numsys.Nbranch
                if from[j] == busPi[i] || to[j] == busPi[i]
                    push!(branch_measure[j], cnt)
                end
            end
        end
    end

    Pf = 1
    @inbounds while sum(Pf) != 0
        ########## Find zero pivots and replace them with angle pseudo-measurements ##########
        G = transpose(H) * H
        F = qr(Matrix(G))
        D = F.R

        whereIs = Int64[]
        for i = 1:numsys.Nbus
            if abs(D[i, i]) <= 1e-5
                push!(whereIs, i)
            end
        end

        if isempty(whereIs)
            break
        end

        t = fill(0.0, numsys.Nbus)
        for (k, i) in enumerate(whereIs)
            H[Nmeasure + i, i] = 1.0
            t[i] = k
        end

        ########## Branch flow estimates ##########
        T = (transpose(H) * H) \ t
        Pf = abs.(A * T) .>= 1e-5
        for i = 1:numsys.Nbus
            for j in A.colptr[i]:(A.colptr[i + 1] - 1)
                if Pf[A.rowval[j]] == 1
                    A[A.rowval[j], i] = 0
                    branch_non[A.rowval[j]] = 0
                end
            end
        end

        ########## Remove all measurements incidence to unobservable branches ##########
        for i = 1:numsys.Nbranch
            if Pf[i] == 1
                for k in branch_measure[i]
                    if rowsdel[k] !=0
                        rowsdel[k] = 0
                        for j in Ht.colptr[k]:(Ht.colptr[k + 1] - 1)
                            H[k, Ht.rowval[j]] = 0.0
                        end
                    end
                end
            end
        end
        for (k, i) in enumerate(whereIs)
            H[Nmeasure + i, i] = 0.0
        end
    end

    islands = connected_components(SimpleGraph(A' * A))

    return islands
end # islandsnodal

max_stars_repo_name: mcosovic/JuliaGrid

max_issues_repo_name: mcosovic/JuliaGrid

max_forks_repo_name: mcosovic/JuliaGrid


text:
@doc raw"""
    ChambollePock(M, N, cost, x0, ξ0, m, n, prox_F, prox_G_dual, forward_operator, adjoint_DΛ)

Perform the Riemannian Chambolle–Pock algorithm.

Given a `cost` function $\mathcal E\colon\mathcal M \to ℝ$ of the form
```math
\mathcal E(x) = F(x) + G( Λ(x) ),
```
where $F\colon\mathcal M \to ℝ$, $G\colon\mathcal N \to ℝ$,
and $\Lambda\colon\mathcal M \to \mathcal N$. The remaining input parameters are

* `x,ξ` primal and dual start points $x\in\mathcal M$ and $\xi\in T_n\mathcal N$
* `m,n` base points on $\mathcal M$ and $\mathcal N$, respectively.
* `forward_operator` the operator $Λ(⋅)$ or its linearization $DΛ(⋅)[⋅]$, depending on whether `:exact` or `:linearized` is chosen.
* `adjoint_linearized_operator` the adjoint $DΛ^*$ of the linearized operator $DΛ(m)\colon T_{m}\mathcal M \to T_{Λ(m)}\mathcal N$
* `prox_F, prox_G_Dual` the proximal maps of $F$ and $G^\ast_n$

By default, this performs the exact Riemannian Chambolle Pock algorithm, see the opional parameter
`DΛ` for ther linearized variant.

For more details on the algorithm, see[^BergmannHerzogSilvaLouzeiroTenbrinckVidalNunez2020].

# Optional Parameters

* `acceleration` – (`0.05`)
* `dual_stepsize` – (`1/sqrt(8)`) proximnal parameter of the primal prox
* `Λ` (`missing`) the exact operator, that is required if the forward operator is linearized;
  `missing` indicates, that the forward operator is exact.
* `primal_stepsize` – (`1/sqrt(8)`) proximnal parameter of the dual prox
* `relaxation` – (`1.`)
* `relax` – (`:primal`) whether to relax the primal or dual
* `variant` - (`:exact` if `Λ` is missing, otherwise `:linearized`) variant to use.
  Note that this changes the arguments the `forward_operator` will be called.
* `stopping_criterion` – (`stopAtIteration(100)`) a [`StoppingCriterion`](@ref)
* `update_primal_base` – (`missing`) function to update `m` (identity by default/missing)
* `update_dual_base` – (`missing`) function to update `n` (identity by default/missing)
* `retraction_method` – (`ExponentialRetraction()`) the rectraction to use
* `inverse_retraction_method` - (`LogarithmicInverseRetraction()`) an inverse retraction to use.
* `vector_transport_method` - (`ParallelTransport()`) a vector transport to use

[^BergmannHerzogSilvaLouzeiroTenbrinckVidalNunez2020]:
    > R. Bergmann, R. Herzog, M. Silva Louzeiro, D. Tenbrinck, J. Vidal-Núñez:
    > _Fenchel Duality Theory and a Primal-Dual Algorithm on Riemannian Manifolds_,
    > arXiv: [1908.02022](http://arxiv.org/abs/1908.02022)
    > accepted for publication in Foundations of Computational Mathematics
"""
function ChambollePock(
    M::Manifold,
    N::Manifold,
    cost::Function,
    x::P,
    ξ::T,
    m::P,
    n::Q,
    prox_F::Function,
    prox_G_dual::Function,
    forward_operator::Function,
    adjoint_linear_operator::Function;
    kwargs...,
) where {P,T,Q}
    x_res = allocate(x)
    copyto!(x_res, x)
    ξ_res = allocate(ξ)
    copyto!(ξ_res, ξ)
    m_res = allocate(m)
    copyto!(m_res, m)
    n_res = allocate(n)
    copyto!(n_res, n)
    return ChambollePock!(
        M,
        N,
        cost,
        x_res,
        ξ_res,
        m_res,
        n_res,
        prox_F,
        prox_G_dual,
        forward_operator,
        adjoint_linear_operator;
        kwargs...,
    )
end
@doc raw"""
    ChambollePock(M, N, cost, x, ξ, m, n, prox_F, prox_G_dual, forward_operator, adjoint_DΛ)

Perform the Riemannian Chambolle–Pock algorithm in place of `x`, `ξ`, and potenitally `m`,
`n` if they are not fixed. See [`ChambollePock`](@ref) for details and optional parameters.
"""
function ChambollePock!(
    M::Manifold,
    N::Manifold,
    cost::Function,
    x::P,
    ξ::T,
    m::P,
    n::Q,
    prox_F::Function,
    prox_G_dual::Function,
    forward_operator::Function,
    adjoint_linear_operator::Function;
    acceleration=0.05,
    dual_stepsize=1 / sqrt(8),
    Λ::Union{Function,Missing}=missing,
    primal_stepsize=1 / sqrt(8),
    relaxation=1.0,
    relax::Symbol=:primal,
    stopping_criterion::StoppingCriterion=StopAfterIteration(200),
    update_primal_base::Union{Function,Missing}=missing,
    update_dual_base::Union{Function,Missing}=missing,
    retraction_method::RM=ExponentialRetraction(),
    inverse_retraction_method::IRM=LogarithmicInverseRetraction(),
    vector_transport_method::VTM=ParallelTransport(),
    variant=ismissing(Λ) ? :exact : :linearized,
    return_options=false,
    kwargs...,
) where {
    P,
    Q,
    T,
    RM<:AbstractRetractionMethod,
    IRM<:AbstractInverseRetractionMethod,
    VTM<:AbstractVectorTransportMethod,
}
    p = PrimalDualProblem(
        M, N, cost, prox_F, prox_G_dual, forward_operator, adjoint_linear_operator, Λ
    )
    o = ChambollePockOptions(
        m,
        n,
        x,
        ξ,
        primal_stepsize,
        dual_stepsize;
        acceleration=acceleration,
        relaxation=relaxation,
        stopping_criterion=stopping_criterion,
        relax=relax,
        update_primal_base=update_primal_base,
        update_dual_base=update_dual_base,
        variant=variant,
        retraction_method=retraction_method,
        inverse_retraction_method=inverse_retraction_method,
        vector_transport_method=vector_transport_method,
    )
    o = decorate_options(o; kwargs...)
    resultO = solve(p, o)
    if return_options
        return resultO
    else
        return get_solver_result(resultO)
    end
end

function initialize_solver!(::PrimalDualProblem, ::ChambollePockOptions) end

function step_solver!(p::PrimalDualProblem, o::ChambollePockOptions, iter)
    primal_dual_step!(p, o, Val(o.relax))
    o.m = ismissing(o.update_primal_base) ? o.m : o.update_primal_base(p, o, iter)
    if !ismissing(o.update_dual_base)
        n_old = deepcopy(o.n)
        o.n = o.update_dual_base(p, o, iter)
        o.ξ = vector_transport_to(p.N, n_old, o.ξ, o.n, o.vector_transport_method)
        o.ξbar = vector_transport_to(p.N, n_old, o.ξbar, o.n, o.vector_transport_method)
    end
    return o
end
#
# Variant 1: primal relax
#
function primal_dual_step!(p::PrimalDualProblem, o::ChambollePockOptions, ::Val{:primal})
    dual_update!(p, o, o.xbar, Val(o.variant))
    ptξn = if ismissing(p.Λ)
        o.ξ
    else
        vector_transport_to(p.N, o.n, o.ξ, p.Λ(o.m), o.vector_transport_method)
    end
    xOld = o.x
    o.x = p.prox_F(
        p.M,
        o.m,
        o.primal_stepsize,
        retract(
            p.M,
            o.x,
            vector_transport_to(
                p.M,
                o.m,
                -o.primal_stepsize * (p.adjoint_linearized_operator(o.m, ptξn)),
                o.x,
                o.vector_transport_method,
            ),
            o.retraction_method,
        ),
    )
    update_prox_parameters!(o)
    o.xbar = retract(
        p.M,
        o.x,
        -o.relaxation * inverse_retract(p.M, o.x, xOld, o.inverse_retraction_method),
        o.retraction_method,
    )
    return o
end
#
# Variant 2: dual relax
#
function primal_dual_step!(p::PrimalDualProblem, o::ChambollePockOptions, ::Val{:dual})
    ptξbar = if ismissing(p.Λ)
        o.ξbar
    else
        vector_transport_to(p.N, o.n, o.ξbar, p.Λ(o.m), o.vector_transport_method)
    end
    o.x = p.prox_F(
        p.M,
        o.m,
        o.primal_stepsize,
        retract(
            p.M,
            o.x,
            vector_transport_to(
                p.M,
                o.m,
                -o.primal_stepsize * (p.adjoint_linearized_operator(o.m, ptξbar)),
                o.x,
                o.vector_transport_method,
            ),
            o.retraction_method,
        ),
    )
    ξ_old = deepcopy(o.ξ)
    dual_update!(p, o, o.x, Val(o.variant))
    update_prox_parameters!(o)
    o.ξbar = o.ξ + o.relaxation * (o.ξ - ξ_old)
    return o
end
#
# Dual step: linearized
# depending on whether its primal relaxed or dual relaxed we start from start=o.x or start=o.xbar here
#
function dual_update!(
    p::PrimalDualProblem, o::ChambollePockOptions, start::P, ::Val{:linearized}
) where {P}
    # (1) compute update direction
    ξ_update = p.forward_operator(
        o.m, inverse_retract(p.M, o.m, start, o.inverse_retraction_method)
    )
    # (2) if p.Λ is missing, we assume that n = Λ(m) and do  not PT, otherwise we do
    ξ_update = if ismissing(p.Λ)
        ξ_update
    else
        vector_transport_to(p.N, p.Λ(o.m), ξ_update, o.n, o.vector_transport_method)
    end
    # (3) to the dual update
    o.ξ = p.prox_G_dual(p.N, o.n, o.dual_stepsize, o.ξ + o.dual_stepsize * ξ_update)
    return o
end
#
# Dual step: exact
# depending on whether its primal relaxed or dual relaxed we start from start=o.x or start=o.xbar here
#
function dual_update!(
    p::PrimalDualProblem, o::ChambollePockOptions, start::P, ::Val{:exact}
) where {P}
    ξ_update = inverse_retract(
        p.N, o.n, p.forward_operator(start), o.inverse_retraction_method
    )
    t = p.prox_G_dual(p.N, o.n, o.dual_stepsize, o.ξ + o.dual_stepsize * ξ_update)
    o.ξ = t
    return o
end

@doc raw"""
    update_prox_parameters!(o)
update the prox parameters as described in Algorithm 2 of Chambolle, Pock, 2010, i.e.

1. ``\theta_{n} = \frac{1}{\sqrt{1+2\gamma\tau_n}}``
2. ``\tau_{n+1} = \theta_n\tau_n``
3. ``\sigma_{n+1} = \frac{\sigma_n}{\theta_n}``
"""
function update_prox_parameters!(o::O) where {O<:PrimalDualOptions}
    if o.acceleration > 0
        o.relaxation = 1 / sqrt(1 + 2 * o.acceleration * o.primal_stepsize)
        o.primal_stepsize = o.primal_stepsize * o.relaxation
        o.dual_stepsize = o.dual_stepsize / o.relaxation
    end
    return o
end


text:
export show_processing, show_src, show_writes

function print_log(notes::Array{String,1}, k::String)
  mm = 60
  println("")
  pl = string("| Time | ", titlecase(k), k == "processing" ? " | Description |\n|:-----|:---------|:------------|\n" : " |\n|:-----|:---------|\n")
  ee = true
  for i = 1:length(notes)
    nn = split(notes[i], " ¦ ", keepempty=true, limit=4)
    (length(nn) < 3) && continue
    L = lastindex(nn[3])
    if nn[2] == k
      (ee == true) && (ee = false)
      func_str = (L > mm) ? (nn[3][firstindex(nn[3]):prevind(nn[3], mm)] * "…") : nn[3]
      if k in ("processing", "write")
        pl *= string("| ", nn[1], "|`", func_str, "`|", nn[4], "|\n")
      else
        pl *= string("| ", nn[1], "|`", func_str, "`|\n")
      end
    end
  end

  if ee
    pl *= (k == "processing") ? "|      | (none)   |             |\n" : "|      | (none)   |\n"
  end
  show(stdout, MIME("text/plain"), Markdown.parse(pl))
  println("")
  return nothing
end

"""
    show_processing(S::GphysData)
    show_processing(S::GphysData, i::Int64)
    show_processing(C::GphysChannel)

Tabulate and print all processing steps in `:notes` to stdout in human-readable format.

See also: `show_src`, `show_writes`, `note!`, `clear_notes!`
"""
function show_processing(S::GphysData)
  for i in 1:S.n
    println("\nChannel ", i)
    print_log(S.notes[i], "processing")
  end
  return nothing
end
show_processing(S::GphysData, i::Int) = print_log(S.notes[i], "processing")
show_processing(C::GphysChannel) = print_log(C.notes, "processing")

"""
    show_src(S::GphysData)
    show_src(S::GphysData, i::Int64)
    show_src(C::GphysChannel)

Tabulate and print all data sources logged in `:notes` to stdout in human-readable format.

See also: `show_processing`, `show_writes`, `note!`, `clear_notes!`
"""
function show_src(S::GphysData)
  for i in 1:S.n
    println("\nChannel ", i)
    print_log(S.notes[i], "+source")
  end
  return nothing
end
show_src(S::GphysData, i::Int) = print_log(S.notes[i], "+source")
show_src(C::GphysChannel) = print_log(C.notes, "+source")

"""
    show_writes(S::GphysData)
    show_writes(S::GphysData, i::Int64)
    show_writes(C::GphysChannel)

Tabulate and print all data writes logged in `:notes` to stdout in human-readable format.

See also: `show_processing`, `show_src`, `note!`, `clear_notes!`
"""
function show_writes(S::GphysData)
  for i in 1:S.n
    println("\nChannel ", i)
    print_log(S.notes[i], "write")
  end
  return nothing
end
show_writes(S::GphysData, i::Int) = print_log(S.notes[i], "write")
show_writes(C::GphysChannel) = print_log(C.notes, "write")


text:
export cg, cg!

cg(A, b, Pl=1; kwargs...) =  cg!(zerox(A,b), A, b, Pl; kwargs...)

function cg!(x, A, b, Pl=1; tol::Real=size(A,2)*eps(), maxiter::Int=size(A,2))
    history = ConvergenceHistory()
    history[:tol] = tol
    reserve!(history,:resnorm, maxiter)

    K = KrylovSubspace(A, length(b), 1, Vector{Adivtype(A,b)}[])
    init!(K, x)
    cg_method!(history, x, K, b, Pl; tol=tol, maxiter=maxiter)
    x, history
end

function cg_method!(log::ConvergenceHistory, x, K::KrylovSubspace, b, Pl=1;
        tol::Real=size(K.A,2)*eps(), maxiter::Integer=size(K.A,2))

    tol = tol * norm(b)
    r = b - nextvec(K)
    p = z = Pl\r
    γ = dot(r, z)
    for iter=1:maxiter
        nextiter!(log, mvps=1)
        append!(K, p)
        q = nextvec(K)
        α = γ/dot(p, q)
        # α>=0 || throw(PosSemidefException("α=$α"))
        update!(x, α, p)
        r -= α*q
        resnorm = norm(r)
        push!(log,:resnorm,resnorm)
        resnorm < tol && break
        z = Pl\r
        oldγ = γ
        γ = dot(r, z)
        β = γ/oldγ
        p = z + β*p
    end
    shrink!(log)
    setconv(log, 0<=norm(r)<tol)
    x
end


text:
using Test
using LogicCircuits
using ProbabilisticCircuits
using DataFrames

@testset "heuristics tests" begin
    data = DataFrame(BitArray([1 0 1 0 1 0 1 0 1 0;
            1 1 1 1 1 1 1 1 1 1;
            0 0 0 0 0 0 0 0 0 0;
            0 1 1 0 1 0 0 1 0 1]))
    pc1, vtree1 = learn_chow_liu_tree_circuit(data)
    @test_throws Exception heuristic_loss(pc1, data;pick_edge="", pick_var="")
    @test_throws Exception heuristic_loss(pc1, data;pick_edge="eFlow", pick_var="")
    for (pick_edge, pick_var) in [("eFlow","vMI"), ("eFlow", "vRand"),
                                        ("eRand","vMI"), ("eRand", "vRand")]
        (or, and), var = heuristic_loss(pc1, data;pick_edge=pick_edge, pick_var=pick_var)
        @test or in linearize(pc1)
        @test issum(or)
        @test ismul(and) || isliteralgate(and)
        @test and in children(or)
        @test var in variables(and)
    end
end


text:
struct DCG_YugeStr <: AbstractProgramme end
typeProgramme = DCG_YugeStr()
nameProgramme = "DCG_YugeStr"
exerProg = Dict()

SquatT1 = Exercise(
    name = squat["low"],
    equipment = equipment["bb"],
    muscles = [
        muscles["quads"],
        muscles["glutes"],
        muscles["spinal erectors"],
        muscles["transverse abdominis"],
    ],
    trainingMax = 115,
)
SquatT2 = Exercise(
    name = squat["front"],
    equipment = equipment["bb"],
    muscles = [muscles["quads"], muscles["glutes"], muscles["core"]],
    trainingMax = 80,
)
push!(exerProg, "SquatT1" => (exercise = SquatT1, progression = deepcopy(DCGYugeStrT1)))
push!(exerProg, "SquatT2" => (exercise = SquatT2, progression = deepcopy(DCGYugeStrT2)))

BenchT1 = Exercise(
    name = bench["flat"],
    equipment = equipment["bb"],
    muscles = [muscles["pecs"], muscles["triceps"], muscles["front delts"]],
    trainingMax = 92.5,
)
BenchT2 = Exercise(
    name = bench["cg"],
    equipment = equipment["bb"],
    muscles = [muscles["pecs"], muscles["triceps"], muscles["front delts"]],
    trainingMax = 92.5,
)
push!(exerProg, "BenchT1" => (exercise = BenchT1, progression = deepcopy(DCGYugeStrT1)))
push!(exerProg, "BenchT2" => (exercise = BenchT2, progression = deepcopy(DCGYugeStrT2)))

DeadliftT1 = Exercise(
    name = deadlift["dl"],
    equipment = equipment["bb"],
    muscles = [
        muscles["hamstrings"],
        muscles["glutes"],
        muscles["spinal erectors"],
        muscles["quads"],
        muscles["traps"],
        muscles["lats"],
        muscles["finger flexors"],
    ],
    trainingMax = 162.5,
)
DeadliftT2 = Exercise(
    name = deadlift["sumo"],
    equipment = equipment["bb"],
    muscles = [
        muscles["quads"],
        muscles["glutes"],
        muscles["hamstrings"],
        muscles["spinal erectors"],
        muscles["traps"],
        muscles["lats"],
        muscles["finger flexors"],
    ],
    trainingMax = 140,
)
push!(
    exerProg,
    "DeadliftT1" => (exercise = DeadliftT1, progression = deepcopy(DCGYugeStrT1)),
)
push!(
    exerProg,
    "DeadliftT2" => (exercise = DeadliftT2, progression = deepcopy(DCGYugeStrT2)),
)

RowT1 = Exercise(
    name = row["pendlay"],
    modality = modality["cheat"],
    equipment = equipment["bb"],
    muscles = [
        muscles["lats"],
        muscles["biceps"],
        muscles["rear delts"],
        muscles["rhomboids"],
        muscles["finger flexors"],
    ],
    trainingMax = 132.5,
)
RowT2 = Exercise(
    name = row["kroc"],
    equipment = equipment["db"],
    muscles = [
        muscles["lats"],
        muscles["biceps"],
        muscles["rear delts"],
        muscles["rhomboids"],
        muscles["finger flexors"],
    ],
    trainingMax = 42.5,
)
push!(exerProg, "RowT1" => (exercise = RowT1, progression = deepcopy(DCGYugeStrT1)))
push!(exerProg, "RowT2" => (exercise = RowT2, progression = deepcopy(DCGYugeStrT2)))

OHPT1 = Exercise(;
    name = press["ohp"],
    equipment = equipment["bb"],
    muscles = [muscles["front delts"], muscles["triceps"], muscles["traps"]],
    trainingMax = 60,
)
OHPT2 = Exercise(;
    name = press["db ohp"],
    equipment = equipment["db"],
    muscles = [muscles["front delts"], muscles["triceps"], muscles["traps"]],
    trainingMax = 50,
)
push!(exerProg, "OHPT1" => (exercise = OHPT1, progression = deepcopy(DCGYugeStrT1)))
push!(exerProg, "OHPT2" => (exercise = OHPT2, progression = deepcopy(DCGYugeStrT2)))

PullupT1 = Exercise(;
    name = pull["pullup"],
    equipment = equipment["pullup bar"],
    muscles = [muscles["lats"], muscles["biceps"], muscles["brachialis"]],
    trainingMax = 12.5,
    roundBase = 1.25,
)
PullupT2 = Exercise(;
    name = pull["archer"],
    equipment = equipment["pullup bar"],
    muscles = [muscles["lats"], muscles["biceps"], muscles["brachialis"]],
    trainingMax = 0,
    roundBase = 1.25,
)
push!(exerProg, "PullupT1" => (exercise = PullupT1, progression = deepcopy(DCGYugeStrT3)))
push!(exerProg, "PullupT2" => (exercise = PullupT2, progression = deepcopy(DCGYugeStrT3)))

Dips = Exercise(;
    name = bench["dips"],
    equipment = [equipment["parallette"], equipment["chain belt"]],
    muscles = [muscles["pecs"], muscles["triceps"], muscles["front delts"]],
    trainingMax = 12.5,
    roundBase = 1.25,
)
InclineDBBench = Exercise(;
    name = bench["db incline"],
    equipment = equipment["db"],
    muscles = [muscles["pecs"], muscles["triceps"], muscles["front delts"]],
    trainingMax = 45,
)
push!(exerProg, "Dips" => (exercise = Dips, progression = deepcopy(DCGYugeStrT3)))
push!(
    exerProg,
    "InclineDBBench" => (exercise = InclineDBBench, progression = deepcopy(DCGYugeStrT3)),
)

HipThrustT1 = Exercise(;
    name = deadlift["hip thrust"],
    equipment = equipment["bb"],
    muscles = [
        muscles["hamstrings"],
        muscles["glutes"],
        muscles["spinal erectors"],
        muscles["quads"],
    ],
    trainingMax = 100,
)
HipThrustT2 = Exercise(;
    name = deadlift["single leg hip thrust"],
    equipment = equipment["bb"],
    muscles = [
        muscles["hamstrings"],
        muscles["glutes"],
        muscles["spinal erectors"],
        muscles["quads"],
    ],
    trainingMax = 10,
)
push!(
    exerProg,
    "HipThrustT1" => (exercise = HipThrustT1, progression = deepcopy(DCGYugeStrT1)),
)
push!(
    exerProg,
    "HipThrustT2" => (exercise = HipThrustT2, progression = deepcopy(DCGYugeStrT2)),
)

Lunges = Exercise(;
    name = squat["lunge"],
    modality = [modality["walking"], modality["backstep"], modality["dragon"]],
    equipment = equipment["db"],
    muscles = [muscles["quads"], muscles["glutes"]],
    trainingMax = 47.5,
)
BulgarianSplit = Exercise(;
    name = squat["bulgarian"],
    equipment = equipment["db"],
    muscles = [muscles["quads"], muscles["glutes"]],
    trainingMax = 37.5,
)
push!(exerProg, "Lunges" => (exercise = Lunges, progression = deepcopy(DCGYugeStrT3)))
push!(
    exerProg,
    "BulgarianSplit" => (exercise = BulgarianSplit, progression = deepcopy(DCGYugeStrT3)),
)

# Grip
WristExt = Exercise(;
    name = grip["wrist extension"],
    equipment = equipment["db"],
    muscles = muscles["wrist extensors"],
    trainingMax = 10,
    roundBase = 1,
)
WristCurl = Exercise(;
    name = grip["wrist curl"],
    equipment = equipment["db"],
    muscles = muscles["wrist extensors"],
    trainingMax = 10,
    roundBase = 1,
)
FingerCurl = Exercise(;
    name = grip["finger curl"],
    equipment = equipment["db"],
    muscles = muscles["finger flexors"],
    trainingMax = 20,
)
PlatePinch = Exercise(;
    name = grip["plate pinch"],
    equipment = equipment["plate"],
    muscles = muscles["hand"],
    trainingMax = 27.5,
)
push!(exerProg, "WristExt" => (exercise = wristExt, progression = deepcopy(DCGYugeStrT4)))
push!(exerProg, "WristCurl" => (exercise = wristCurl, progression = deepcopy(DCGYugeStrT4)))
push!(
    exerProg,
    "FingerCurl" => (exercise = FingerCurl, progression = deepcopy(DCGYugeStrT4)),
)
push!(
    exerProg,
    "PlatePinch" => (exercise = PlatePinch, progression = deepcopy(DCGYugeStrT4)),
)

# Core
DragonFlag = Exercise(;
    name = core["dragon flag"],
    equipment = equipment["none"],
    muscles = muscles["core"],
    trainingMax = 10,
)
ParaCrunch = Exercise(;
    name = core["v sit"],
    equipment = equipment["parallette"],
    muscles = muscles["core"],
    trainingMax = 10,
)
LegRaise = Exercise(;
    name = core["leg raise"],
    equipment = equipment["pullup bar"],
    muscles = muscles["core"],
    trainingMax = 12.5,
)
AbRollout = Exercise(;
    name = core["ab roll"],
    modality = modality["kneeling"],
    equipment = equipment["ab wheel"],
    muscles = muscles["core"],
    trainingMax = 10,
)
LandMine = Exercise(;
    name = core["landmine"],
    equipment = equipment["bb"],
    muscles = muscles["core"],
    trainingMax = 50,
)
ParaCrunch = Exercise(;
    name = core["parallette crunch"],
    equipment = [equipment["parallette"], equipment["db"]],
    muscles = muscles["core"],
    trainingMax = 15,
)
Vsit = Exercise(;
    name = core["v sit"],
    equipment = equipment["parallette"],
    muscles = muscles["core"],
    trainingMax = 10,
)
push!(exerProg, "DragonFlag" => (exercise = DragonFlag, progression = deepcopy(AMRAP_Cali)))
push!(exerProg, "ParaCrunch" => (exercise = ParaCrunch, progression = deepcopy(AMRAP_Cali)))
push!(exerProg, "LegRaise" => (exercise = LegRaise, progression = deepcopy(DCGYugeStrT3)))
push!(exerProg, "AbRollout" => (exercise = AbRollout, progression = deepcopy(AMRAP_Cali)))
push!(exerProg, "LandMine" => (exercise = LandMine, progression = deepcopy(DCGYugeStrT3)))
push!(
    exerProg,
    "ParaCrunch" => (exercise = ParaCrunch, progression = deepcopy(DCGYugeStrT3)),
)
push!(exerProg, "Vsit" => (exercise = Vsit, progression = deepcopy(AMRAP_Cali)))

# Calisthenics Isometrics
FrontLever = Exercise(;
    name = isometric["front"],
    modality = modality["one leg out"],
    equipment = equipment["pullup bar"],
    muscles = [muscles["lats"], muscles["core"]],
    trainingMax = 15,
)
BackLever = Exercise(;
    name = isometric["back"],
    modality = modality["straddle"],
    equipment = equipment["pullup bar"],
    muscles = [muscles["lats"], muscles["pecs"], muscles["front delts"], muscles["core"]],
    trainingMax = 15,
)
push!(exerProg, "FrontLever" => (exercise = FrontLever, progression = deepcopy(AMRAP_Cali)))
push!(exerProg, "BackLever" => (exercise = BackLever, progression = deepcopy(AMRAP_Cali)))

# Misc
FacePull = Exercise(;
    name = pull["face"],
    equipment = equipment["cable"],
    muscles = [muscles["rear delts"], muscles["rhomboids"], muscles["rot cuff"]],
    trainingMax = 30,
)
HammerCurl = Exercise(;
    name = pull["hammer"],
    equipment = [equipment["db"], equipment["fat grips"]],
    muscles = [muscles["brachialis"], muscles["biceps"]],
    trainingMax = 20,
)
ZottmanCurl = Exercise(;
    name = pull["zottman"],
    equipment = [equipment["db"], equipment["fat grips"]],
    muscles = [muscles["biceps"], muscles["brachialis"]],
    trainingMax = 12.5,
)
SideRaise = Exercise(;
    name = press["side raise"],
    equipment = equipment["db"],
    muscles = muscles["side delts"],
    trainingMax = 12.5,
)
CalfRaise = Exercise(;
    name = squat["calf raise"],
    modality = modality["standing"],
    equipment = equipment["db"],
    muscles = muscles["calves"],
    trainingMax = 20,
)
push!(exerProg, "FacePull" => (exercise = FacePull, progression = deepcopy(DCGYugeStrT4)))
push!(
    exerProg,
    "HammerCurl" => (exercise = HammerCurl, progression = deepcopy(DCGYugeStrT4)),
)
push!(
    exerProg,
    "ZottmanCurl" => (exercise = ZottmanCurl, progression = deepcopy(DCGYugeStrT4)),
)
push!(exerProg, "SideRaise" => (exercise = SideRaise, progression = deepcopy(DCGYugeStrT4)))
push!(exerProg, "CalfRaise" => (exercise = CalfRaise, progression = deepcopy(DCGYugeStrT4)))

import Lifting: makeDays
function makeDays(::DCG_YugeStr, dict::Dict)
    week1 = [[], [], [], [], [], [], []]
    week2 = [[], [], [], [], [], [], []]
    week3 = [[], [], [], [], [], [], []]
    week4 = [[], [], [], [], [], [], []]
    week5 = [[], [], [], [], [], [], []]
    week6 = [[], [], [], [], [], [], []]
    week7 = [[], [], [], [], [], [], []]
    week8 = [[], [], [], [], [], [], []]
    week9 = [[], [], [], [], [], [], []]

    #=
        Day 1
            bench 1
            wrist extension
            wrist flexion
            dragon flag
            row 2
            front lever
            lunges
            face pull
    =#
    push!(week1[1], dict["BenchT1"][1], dict["BenchT1"][2], 1)
    push!(week2[1], dict["BenchT1"][1], dict["BenchT1"][2], 2)
    push!(week3[1], dict["BenchT1"][1], dict["BenchT1"][2], 3)
    push!(week4[1], dict["BenchT1"][1], dict["BenchT1"][2], 4)
    push!(week5[1], dict["BenchT1"][1], dict["BenchT1"][2], 5)
    push!(week6[1], dict["BenchT1"][1], dict["BenchT1"][2], 6)
    push!(week7[1], dict["BenchT1"][1], dict["BenchT1"][2], 7)
    push!(week8[1], dict["BenchT1"][1], dict["BenchT1"][2], 8)
    push!(week9[1], dict["BenchT1"][1], dict["BenchT1"][2], 9)

    push!(week1[1], dict["WristExt"][1], dict["WristExt"][2], 2)
    push!(week2[1], dict["WristExt"][1], dict["WristExt"][2], 3)
    push!(week3[1], dict["WristExt"][1], dict["WristExt"][2], 4)
    push!(week4[1], dict["WristExt"][1], dict["WristExt"][2], 5)
    push!(week5[1], dict["WristExt"][1], dict["WristExt"][2], 6)
    push!(week6[1], dict["WristExt"][1], dict["WristExt"][2], 7)
    push!(week7[1], dict["WristExt"][1], dict["WristExt"][2], 8)
    push!(week8[1], dict["WristExt"][1], dict["WristExt"][2], 9)

    push!(week1[1], dict["WristCurl"][1], dict["WristCurl"][2], 2)
    push!(week2[1], dict["WristCurl"][1], dict["WristCurl"][2], 3)
    push!(week3[1], dict["WristCurl"][1], dict["WristCurl"][2], 4)
    push!(week4[1], dict["WristCurl"][1], dict["WristCurl"][2], 5)
    push!(week5[1], dict["WristCurl"][1], dict["WristCurl"][2], 6)
    push!(week6[1], dict["WristCurl"][1], dict["WristCurl"][2], 7)
    push!(week7[1], dict["WristCurl"][1], dict["WristCurl"][2], 8)
    push!(week8[1], dict["WristCurl"][1], dict["WristCurl"][2], 9)

    push!(week1[1], dict["RowT2"][1], dict["RowT2"][2], 1)
    push!(week2[1], dict["RowT2"][1], dict["RowT2"][2], 2)
    push!(week3[1], dict["RowT2"][1], dict["RowT2"][2], 3)
    push!(week4[1], dict["RowT2"][1], dict["RowT2"][2], 4)
    push!(week5[1], dict["RowT2"][1], dict["RowT2"][2], 5)
    push!(week6[1], dict["RowT2"][1], dict["RowT2"][2], 6)
    push!(week7[1], dict["RowT2"][1], dict["RowT2"][2], 7)
    push!(week8[1], dict["RowT2"][1], dict["RowT2"][2], 8)
    push!(week9[1], dict["RowT2"][1], dict["RowT2"][2], 9)

    push!(week1[1], dict["DragonFlag"][1], dict["DragonFlag"][2])
    push!(week2[1], dict["DragonFlag"][1], dict["DragonFlag"][2])
    push!(week3[1], dict["DragonFlag"][1], dict["DragonFlag"][2])
    push!(week4[1], dict["DragonFlag"][1], dict["DragonFlag"][2])
    push!(week5[1], dict["DragonFlag"][1], dict["DragonFlag"][2])
    push!(week6[1], dict["DragonFlag"][1], dict["DragonFlag"][2])
    push!(week7[1], dict["DragonFlag"][1], dict["DragonFlag"][2])
    push!(week8[1], dict["DragonFlag"][1], dict["DragonFlag"][2])

    push!(week1[1], dict["FrontLever"][1], dict["FrontLever"][2])
    push!(week2[1], dict["FrontLever"][1], dict["FrontLever"][2])
    push!(week3[1], dict["FrontLever"][1], dict["FrontLever"][2])
    push!(week4[1], dict["FrontLever"][1], dict["FrontLever"][2])
    push!(week5[1], dict["FrontLever"][1], dict["FrontLever"][2])
    push!(week6[1], dict["FrontLever"][1], dict["FrontLever"][2])
    push!(week7[1], dict["FrontLever"][1], dict["FrontLever"][2])
    push!(week8[1], dict["FrontLever"][1], dict["FrontLever"][2])

    push!(week1[1], dict["Lunges"][1], dict["Lunges"][2], 2)
    push!(week2[1], dict["Lunges"][1], dict["Lunges"][2], 3)
    push!(week3[1], dict["Lunges"][1], dict["Lunges"][2], 4)
    push!(week4[1], dict["Lunges"][1], dict["Lunges"][2], 5)
    push!(week5[1], dict["Lunges"][1], dict["Lunges"][2], 6)
    push!(week6[1], dict["Lunges"][1], dict["Lunges"][2], 7)
    push!(week7[1], dict["Lunges"][1], dict["Lunges"][2], 8)
    push!(week8[1], dict["Lunges"][1], dict["Lunges"][2], 9)

    push!(week1[1], dict["FacePull"][1], dict["FacePull"][2], 1)
    push!(week2[1], dict["FacePull"][1], dict["FacePull"][2], 2)
    push!(week3[1], dict["FacePull"][1], dict["FacePull"][2], 3)
    push!(week4[1], dict["FacePull"][1], dict["FacePull"][2], 4)
    push!(week5[1], dict["FacePull"][1], dict["FacePull"][2], 5)
    push!(week6[1], dict["FacePull"][1], dict["FacePull"][2], 6)
    push!(week7[1], dict["FacePull"][1], dict["FacePull"][2], 7)
    push!(week8[1], dict["FacePull"][1], dict["FacePull"][2], 8)
    push!(week9[1], dict["FacePull"][1], dict["FacePull"][2], 9)
    #=
        Day 2
            squat 2
            ohp 1
            pullup 2
            back lever
            calf raises
            finger curls
            leg raises
    =#
    push!(week1[2], dict["SquatT2"][1], dict["SquatT2"][2], 1)
    push!(week2[2], dict["SquatT2"][1], dict["SquatT2"][2], 2)
    push!(week3[2], dict["SquatT2"][1], dict["SquatT2"][2], 3)
    push!(week4[2], dict["SquatT2"][1], dict["SquatT2"][2], 4)
    push!(week5[2], dict["SquatT2"][1], dict["SquatT2"][2], 5)
    push!(week6[2], dict["SquatT2"][1], dict["SquatT2"][2], 6)
    push!(week7[2], dict["SquatT2"][1], dict["SquatT2"][2], 7)
    push!(week8[2], dict["SquatT2"][1], dict["SquatT2"][2], 8)
    push!(week9[2], dict["SquatT2"][1], dict["SquatT2"][2], 9)

    push!(week1[2], dict["OHPT1"][1], dict["OHPT1"][2], 1)
    push!(week2[2], dict["OHPT1"][1], dict["OHPT1"][2], 2)
    push!(week3[2], dict["OHPT1"][1], dict["OHPT1"][2], 3)
    push!(week4[2], dict["OHPT1"][1], dict["OHPT1"][2], 4)
    push!(week5[2], dict["OHPT1"][1], dict["OHPT1"][2], 5)
    push!(week6[2], dict["OHPT1"][1], dict["OHPT1"][2], 6)
    push!(week7[2], dict["OHPT1"][1], dict["OHPT1"][2], 7)
    push!(week8[2], dict["OHPT1"][1], dict["OHPT1"][2], 8)
    push!(week9[2], dict["OHPT1"][1], dict["OHPT1"][2], 9)

    push!(week1[2], dict["PullupT2"][1], dict["PullupT2"][2], 1)
    push!(week2[2], dict["PullupT2"][1], dict["PullupT2"][2], 2)
    push!(week3[2], dict["PullupT2"][1], dict["PullupT2"][2], 3)
    push!(week4[2], dict["PullupT2"][1], dict["PullupT2"][2], 4)
    push!(week5[2], dict["PullupT2"][1], dict["PullupT2"][2], 5)
    push!(week6[2], dict["PullupT2"][1], dict["PullupT2"][2], 6)
    push!(week7[2], dict["PullupT2"][1], dict["PullupT2"][2], 7)
    push!(week8[2], dict["PullupT2"][1], dict["PullupT2"][2], 8)
    push!(week9[2], dict["PullupT2"][1], dict["PullupT2"][2], 9)

    push!(week1[2], dict["BackLever"][1], dict["BackLever"][2])
    push!(week2[2], dict["BackLever"][1], dict["BackLever"][2])
    push!(week3[2], dict["BackLever"][1], dict["BackLever"][2])
    push!(week4[2], dict["BackLever"][1], dict["BackLever"][2])
    push!(week5[2], dict["BackLever"][1], dict["BackLever"][2])
    push!(week6[2], dict["BackLever"][1], dict["BackLever"][2])
    push!(week7[2], dict["BackLever"][1], dict["BackLever"][2])
    push!(week8[2], dict["BackLever"][1], dict["BackLever"][2])

    push!(week1[2], dict["CalfRaise"][1], dict["CalfRaise"][2], 2)
    push!(week2[2], dict["CalfRaise"][1], dict["CalfRaise"][2], 3)
    push!(week3[2], dict["CalfRaise"][1], dict["CalfRaise"][2], 4)
    push!(week4[2], dict["CalfRaise"][1], dict["CalfRaise"][2], 5)
    push!(week5[2], dict["CalfRaise"][1], dict["CalfRaise"][2], 6)
    push!(week6[2], dict["CalfRaise"][1], dict["CalfRaise"][2], 7)
    push!(week7[2], dict["CalfRaise"][1], dict["CalfRaise"][2], 8)
    push!(week8[2], dict["CalfRaise"][1], dict["CalfRaise"][2], 9)

    push!(week1[2], dict["FingerCurl"][1], dict["FingerCurl"][2], 2)
    push!(week2[2], dict["FingerCurl"][1], dict["FingerCurl"][2], 3)
    push!(week3[2], dict["FingerCurl"][1], dict["FingerCurl"][2], 4)
    push!(week4[2], dict["FingerCurl"][1], dict["FingerCurl"][2], 5)
    push!(week5[2], dict["FingerCurl"][1], dict["FingerCurl"][2], 6)
    push!(week6[2], dict["FingerCurl"][1], dict["FingerCurl"][2], 7)
    push!(week7[2], dict["FingerCurl"][1], dict["FingerCurl"][2], 8)
    push!(week8[2], dict["FingerCurl"][1], dict["FingerCurl"][2], 9)

    push!(week1[2], dict["LegRaise"][1], dict["LegRaise"][2], 2)
    push!(week2[2], dict["LegRaise"][1], dict["LegRaise"][2], 3)
    push!(week3[2], dict["LegRaise"][1], dict["LegRaise"][2], 4)
    push!(week4[2], dict["LegRaise"][1], dict["LegRaise"][2], 5)
    push!(week5[2], dict["LegRaise"][1], dict["LegRaise"][2], 6)
    push!(week6[2], dict["LegRaise"][1], dict["LegRaise"][2], 7)
    push!(week7[2], dict["LegRaise"][1], dict["LegRaise"][2], 8)
    push!(week8[2], dict["LegRaise"][1], dict["LegRaise"][2], 9)

    #=
        Day 3
            deadlift 1
            hip thrust 2
            side raise
            dips
            bulgarian split
            plate pinch
            hammer curls
            v-sit
    =#
    push!(week1[3], dict["DeadliftT1"][1], dict["DeadliftT1"][2], 1)
    push!(week2[3], dict["DeadliftT1"][1], dict["DeadliftT1"][2], 2)
    push!(week3[3], dict["DeadliftT1"][1], dict["DeadliftT1"][2], 3)
    push!(week4[3], dict["DeadliftT1"][1], dict["DeadliftT1"][2], 4)
    push!(week5[3], dict["DeadliftT1"][1], dict["DeadliftT1"][2], 5)
    push!(week6[3], dict["DeadliftT1"][1], dict["DeadliftT1"][2], 6)
    push!(week7[3], dict["DeadliftT1"][1], dict["DeadliftT1"][2], 7)
    push!(week8[3], dict["DeadliftT1"][1], dict["DeadliftT1"][2], 8)
    push!(week9[3], dict["DeadliftT1"][1], dict["DeadliftT1"][2], 9)

    push!(week1[3], dict["HipThrustT2"][1], dict["HipThrustT2"][2], 1)
    push!(week2[3], dict["HipThrustT2"][1], dict["HipThrustT2"][2], 2)
    push!(week3[3], dict["HipThrustT2"][1], dict["HipThrustT2"][2], 3)
    push!(week4[3], dict["HipThrustT2"][1], dict["HipThrustT2"][2], 4)
    push!(week5[3], dict["HipThrustT2"][1], dict["HipThrustT2"][2], 5)
    push!(week6[3], dict["HipThrustT2"][1], dict["HipThrustT2"][2], 6)
    push!(week7[3], dict["HipThrustT2"][1], dict["HipThrustT2"][2], 7)
    push!(week8[3], dict["HipThrustT2"][1], dict["HipThrustT2"][2], 8)
    push!(week9[3], dict["HipThrustT2"][1], dict["HipThrustT2"][2], 9)

    push!(week1[3], dict["SideRaise"][1], dict["SideRaise"][2], 2)
    push!(week2[3], dict["SideRaise"][1], dict["SideRaise"][2], 3)
    push!(week3[3], dict["SideRaise"][1], dict["SideRaise"][2], 4)
    push!(week4[3], dict["SideRaise"][1], dict["SideRaise"][2], 5)
    push!(week5[3], dict["SideRaise"][1], dict["SideRaise"][2], 6)
    push!(week6[3], dict["SideRaise"][1], dict["SideRaise"][2], 7)
    push!(week7[3], dict["SideRaise"][1], dict["SideRaise"][2], 8)
    push!(week8[3], dict["SideRaise"][1], dict["SideRaise"][2], 9)

    push!(week1[3], dict["Dips"][1], dict["Dips"][2], 1)
    push!(week2[3], dict["Dips"][1], dict["Dips"][2], 2)
    push!(week3[3], dict["Dips"][1], dict["Dips"][2], 3)
    push!(week4[3], dict["Dips"][1], dict["Dips"][2], 4)
    push!(week5[3], dict["Dips"][1], dict["Dips"][2], 5)
    push!(week6[3], dict["Dips"][1], dict["Dips"][2], 6)
    push!(week7[3], dict["Dips"][1], dict["Dips"][2], 7)
    push!(week8[3], dict["Dips"][1], dict["Dips"][2], 8)
    push!(week9[3], dict["Dips"][1], dict["Dips"][2], 9)

    push!(week1[3], dict["BulgarianSplit"][1], dict["BulgarianSplit"][2], 2)
    push!(week2[3], dict["BulgarianSplit"][1], dict["BulgarianSplit"][2], 3)
    push!(week3[3], dict["BulgarianSplit"][1], dict["BulgarianSplit"][2], 4)
    push!(week4[3], dict["BulgarianSplit"][1], dict["BulgarianSplit"][2], 5)
    push!(week5[3], dict["BulgarianSplit"][1], dict["BulgarianSplit"][2], 6)
    push!(week6[3], dict["BulgarianSplit"][1], dict["BulgarianSplit"][2], 7)
    push!(week7[3], dict["BulgarianSplit"][1], dict["BulgarianSplit"][2], 8)
    push!(week8[3], dict["BulgarianSplit"][1], dict["BulgarianSplit"][2], 9)

    push!(week1[3], dict["PlatePinch"][1], dict["PlatePinch"][2], 2)
    push!(week2[3], dict["PlatePinch"][1], dict["PlatePinch"][2], 3)
    push!(week3[3], dict["PlatePinch"][1], dict["PlatePinch"][2], 4)
    push!(week4[3], dict["PlatePinch"][1], dict["PlatePinch"][2], 5)
    push!(week5[3], dict["PlatePinch"][1], dict["PlatePinch"][2], 6)
    push!(week6[3], dict["PlatePinch"][1], dict["PlatePinch"][2], 7)
    push!(week7[3], dict["PlatePinch"][1], dict["PlatePinch"][2], 8)
    push!(week8[3], dict["PlatePinch"][1], dict["PlatePinch"][2], 9)

    push!(week1[3], dict["HammerCurl"][1], dict["HammerCurl"][2], 2)
    push!(week2[3], dict["HammerCurl"][1], dict["HammerCurl"][2], 3)
    push!(week3[3], dict["HammerCurl"][1], dict["HammerCurl"][2], 4)
    push!(week4[3], dict["HammerCurl"][1], dict["HammerCurl"][2], 5)
    push!(week5[3], dict["HammerCurl"][1], dict["HammerCurl"][2], 6)
    push!(week6[3], dict["HammerCurl"][1], dict["HammerCurl"][2], 7)
    push!(week7[3], dict["HammerCurl"][1], dict["HammerCurl"][2], 8)
    push!(week8[3], dict["HammerCurl"][1], dict["HammerCurl"][2], 9)

    push!(week1[3], dict["Vsit"][1], dict["Vsit"][2])
    push!(week2[3], dict["Vsit"][1], dict["Vsit"][2])
    push!(week3[3], dict["Vsit"][1], dict["Vsit"][2])
    push!(week4[3], dict["Vsit"][1], dict["Vsit"][2])
    push!(week5[3], dict["Vsit"][1], dict["Vsit"][2])
    push!(week6[3], dict["Vsit"][1], dict["Vsit"][2])
    push!(week7[3], dict["Vsit"][1], dict["Vsit"][2])
    push!(week8[3], dict["Vsit"][1], dict["Vsit"][2])

    #=
        Day 4
            Rest
    =#
    push!(week1[4], "Rest")
    push!(week2[4], "Rest")
    push!(week3[4], "Rest")
    push!(week4[4], "Rest")
    push!(week5[4], "Rest")
    push!(week6[4], "Rest")
    push!(week7[4], "Rest")
    push!(week8[4], "Rest")
    push!(week9[4], "Rest")

    #=
        Day 5
            bench 2
            land mine anti rotation
            row 1
            front lever
            lunges
            wrist extension
            wrist flexion
            face pull
    =#
    push!(week1[5], dict["BenchT2"][1], dict["BenchT2"][2], 1)
    push!(week2[5], dict["BenchT2"][1], dict["BenchT2"][2], 2)
    push!(week3[5], dict["BenchT2"][1], dict["BenchT2"][2], 3)
    push!(week4[5], dict["BenchT2"][1], dict["BenchT2"][2], 4)
    push!(week5[5], dict["BenchT2"][1], dict["BenchT2"][2], 5)
    push!(week6[5], dict["BenchT2"][1], dict["BenchT2"][2], 6)
    push!(week7[5], dict["BenchT2"][1], dict["BenchT2"][2], 7)
    push!(week8[5], dict["BenchT2"][1], dict["BenchT2"][2], 8)
    push!(week9[5], dict["BenchT2"][1], dict["BenchT2"][2], 9)

    push!(week1[5], dict["WristExt"][1], dict["WristExt"][2], 2)
    push!(week2[5], dict["WristExt"][1], dict["WristExt"][2], 3)
    push!(week3[5], dict["WristExt"][1], dict["WristExt"][2], 4)
    push!(week4[5], dict["WristExt"][1], dict["WristExt"][2], 5)
    push!(week5[5], dict["WristExt"][1], dict["WristExt"][2], 6)
    push!(week6[5], dict["WristExt"][1], dict["WristExt"][2], 7)
    push!(week7[5], dict["WristExt"][1], dict["WristExt"][2], 8)
    push!(week8[5], dict["WristExt"][1], dict["WristExt"][2], 9)

    push!(week1[5], dict["WristCurl"][1], dict["WristCurl"][2], 2)
    push!(week2[5], dict["WristCurl"][1], dict["WristCurl"][2], 3)
    push!(week3[5], dict["WristCurl"][1], dict["WristCurl"][2], 4)
    push!(week4[5], dict["WristCurl"][1], dict["WristCurl"][2], 5)
    push!(week5[5], dict["WristCurl"][1], dict["WristCurl"][2], 6)
    push!(week6[5], dict["WristCurl"][1], dict["WristCurl"][2], 7)
    push!(week7[5], dict["WristCurl"][1], dict["WristCurl"][2], 8)
    push!(week8[5], dict["WristCurl"][1], dict["WristCurl"][2], 9)

    push!(week1[5], dict["RowT1"][1], dict["RowT1"][2], 1)
    push!(week2[5], dict["RowT1"][1], dict["RowT1"][2], 2)
    push!(week3[5], dict["RowT1"][1], dict["RowT1"][2], 3)
    push!(week4[5], dict["RowT1"][1], dict["RowT1"][2], 4)
    push!(week5[5], dict["RowT1"][1], dict["RowT1"][2], 5)
    push!(week6[5], dict["RowT1"][1], dict["RowT1"][2], 6)
    push!(week7[5], dict["RowT1"][1], dict["RowT1"][2], 7)
    push!(week8[5], dict["RowT1"][1], dict["RowT1"][2], 8)
    push!(week9[5], dict["RowT1"][1], dict["RowT1"][2], 9)

    push!(week1[5], dict["LandMine"][1], dict["LandMine"][2], 2)
    push!(week2[5], dict["LandMine"][1], dict["LandMine"][2], 3)
    push!(week3[5], dict["LandMine"][1], dict["LandMine"][2], 4)
    push!(week4[5], dict["LandMine"][1], dict["LandMine"][2], 5)
    push!(week5[5], dict["LandMine"][1], dict["LandMine"][2], 6)
    push!(week6[5], dict["LandMine"][1], dict["LandMine"][2], 7)
    push!(week7[5], dict["LandMine"][1], dict["LandMine"][2], 8)
    push!(week8[5], dict["LandMine"][1], dict["LandMine"][2], 9)

    push!(week1[5], dict["FrontLever"][1], dict["FrontLever"][2])
    push!(week2[5], dict["FrontLever"][1], dict["FrontLever"][2])
    push!(week3[5], dict["FrontLever"][1], dict["FrontLever"][2])
    push!(week4[5], dict["FrontLever"][1], dict["FrontLever"][2])
    push!(week5[5], dict["FrontLever"][1], dict["FrontLever"][2])
    push!(week6[5], dict["FrontLever"][1], dict["FrontLever"][2])
    push!(week7[5], dict["FrontLever"][1], dict["FrontLever"][2])
    push!(week8[5], dict["FrontLever"][1], dict["FrontLever"][2])

    push!(week1[5], dict["Lunges"][1], dict["Lunges"][2], 2)
    push!(week2[5], dict["Lunges"][1], dict["Lunges"][2], 3)
    push!(week3[5], dict["Lunges"][1], dict["Lunges"][2], 4)
    push!(week4[5], dict["Lunges"][1], dict["Lunges"][2], 5)
    push!(week5[5], dict["Lunges"][1], dict["Lunges"][2], 6)
    push!(week6[5], dict["Lunges"][1], dict["Lunges"][2], 7)
    push!(week7[5], dict["Lunges"][1], dict["Lunges"][2], 8)
    push!(week8[5], dict["Lunges"][1], dict["Lunges"][2], 9)

    push!(week1[5], dict["FacePull"][1], dict["FacePull"][2], 1)
    push!(week2[5], dict["FacePull"][1], dict["FacePull"][2], 2)
    push!(week3[5], dict["FacePull"][1], dict["FacePull"][2], 3)
    push!(week4[5], dict["FacePull"][1], dict["FacePull"][2], 4)
    push!(week5[5], dict["FacePull"][1], dict["FacePull"][2], 5)
    push!(week6[5], dict["FacePull"][1], dict["FacePull"][2], 6)
    push!(week7[5], dict["FacePull"][1], dict["FacePull"][2], 7)
    push!(week8[5], dict["FacePull"][1], dict["FacePull"][2], 8)
    push!(week9[5], dict["FacePull"][1], dict["FacePull"][2], 9)

    #=
        Day 6
            squat 1
            ohp 2
            pullup
            back lever
            calf raises
            finger curls
            ab wheel
    =#
    push!(week1[6], dict["SquatT1"][1], dict["SquatT1"][2], 1)
    push!(week2[6], dict["SquatT1"][1], dict["SquatT1"][2], 2)
    push!(week3[6], dict["SquatT1"][1], dict["SquatT1"][2], 3)
    push!(week4[6], dict["SquatT1"][1], dict["SquatT1"][2], 4)
    push!(week5[6], dict["SquatT1"][1], dict["SquatT1"][2], 5)
    push!(week6[6], dict["SquatT1"][1], dict["SquatT1"][2], 6)
    push!(week7[6], dict["SquatT1"][1], dict["SquatT1"][2], 7)
    push!(week8[6], dict["SquatT1"][1], dict["SquatT1"][2], 8)
    push!(week9[6], dict["SquatT1"][1], dict["SquatT1"][2], 9)

    push!(week1[6], dict["OHPT2"][1], dict["OHPT2"][2], 1)
    push!(week2[6], dict["OHPT2"][1], dict["OHPT2"][2], 2)
    push!(week3[6], dict["OHPT2"][1], dict["OHPT2"][2], 3)
    push!(week4[6], dict["OHPT2"][1], dict["OHPT2"][2], 4)
    push!(week5[6], dict["OHPT2"][1], dict["OHPT2"][2], 5)
    push!(week6[6], dict["OHPT2"][1], dict["OHPT2"][2], 6)
    push!(week7[6], dict["OHPT2"][1], dict["OHPT2"][2], 7)
    push!(week8[6], dict["OHPT2"][1], dict["OHPT2"][2], 8)
    push!(week9[6], dict["OHPT2"][1], dict["OHPT2"][2], 9)

    push!(week1[6], dict["PullupT1"][1], dict["PullupT1"][2], 1)
    push!(week2[6], dict["PullupT1"][1], dict["PullupT1"][2], 2)
    push!(week3[6], dict["PullupT1"][1], dict["PullupT1"][2], 3)
    push!(week4[6], dict["PullupT1"][1], dict["PullupT1"][2], 4)
    push!(week5[6], dict["PullupT1"][1], dict["PullupT1"][2], 5)
    push!(week6[6], dict["PullupT1"][1], dict["PullupT1"][2], 6)
    push!(week7[6], dict["PullupT1"][1], dict["PullupT1"][2], 7)
    push!(week8[6], dict["PullupT1"][1], dict["PullupT1"][2], 8)
    push!(week9[6], dict["PullupT1"][1], dict["PullupT1"][2], 9)

    push!(week1[6], dict["BackLever"][1], dict["BackLever"][2])
    push!(week2[6], dict["BackLever"][1], dict["BackLever"][2])
    push!(week3[6], dict["BackLever"][1], dict["BackLever"][2])
    push!(week4[6], dict["BackLever"][1], dict["BackLever"][2])
    push!(week5[6], dict["BackLever"][1], dict["BackLever"][2])
    push!(week6[6], dict["BackLever"][1], dict["BackLever"][2])
    push!(week7[6], dict["BackLever"][1], dict["BackLever"][2])
    push!(week8[6], dict["BackLever"][1], dict["BackLever"][2])

    push!(week1[6], dict["CalfRaise"][1], dict["CalfRaise"][2], 2)
    push!(week2[6], dict["CalfRaise"][1], dict["CalfRaise"][2], 3)
    push!(week3[6], dict["CalfRaise"][1], dict["CalfRaise"][2], 4)
    push!(week4[6], dict["CalfRaise"][1], dict["CalfRaise"][2], 5)
    push!(week5[6], dict["CalfRaise"][1], dict["CalfRaise"][2], 6)
    push!(week6[6], dict["CalfRaise"][1], dict["CalfRaise"][2], 7)
    push!(week7[6], dict["CalfRaise"][1], dict["CalfRaise"][2], 8)
    push!(week8[6], dict["CalfRaise"][1], dict["CalfRaise"][2], 9)

    push!(week1[6], dict["FingerCurl"][1], dict["FingerCurl"][2], 2)
    push!(week2[6], dict["FingerCurl"][1], dict["FingerCurl"][2], 3)
    push!(week3[6], dict["FingerCurl"][1], dict["FingerCurl"][2], 4)
    push!(week4[6], dict["FingerCurl"][1], dict["FingerCurl"][2], 5)
    push!(week5[6], dict["FingerCurl"][1], dict["FingerCurl"][2], 6)
    push!(week6[6], dict["FingerCurl"][1], dict["FingerCurl"][2], 7)
    push!(week7[6], dict["FingerCurl"][1], dict["FingerCurl"][2], 8)
    push!(week8[6], dict["FingerCurl"][1], dict["FingerCurl"][2], 9)

    push!(week1[6], dict["AbRollout"][1], dict["AbRollout"][2])
    push!(week2[6], dict["AbRollout"][1], dict["AbRollout"][2])
    push!(week3[6], dict["AbRollout"][1], dict["AbRollout"][2])
    push!(week4[6], dict["AbRollout"][1], dict["AbRollout"][2])
    push!(week5[6], dict["AbRollout"][1], dict["AbRollout"][2])
    push!(week6[6], dict["AbRollout"][1], dict["AbRollout"][2])
    push!(week7[6], dict["AbRollout"][1], dict["AbRollout"][2])
    push!(week8[6], dict["AbRollout"][1], dict["AbRollout"][2])

    #=
        Day 7
            deadlift 2
            hip thrust 1
            side raise
            incline db bench
            bulgarian split
            plate pinch
            zottman curls
            paralette crunch
    =#
    push!(week1[7], dict["DeadliftT2"][1], dict["DeadliftT2"][2], 1)
    push!(week2[7], dict["DeadliftT2"][1], dict["DeadliftT2"][2], 2)
    push!(week3[7], dict["DeadliftT2"][1], dict["DeadliftT2"][2], 3)
    push!(week4[7], dict["DeadliftT2"][1], dict["DeadliftT2"][2], 4)
    push!(week5[7], dict["DeadliftT2"][1], dict["DeadliftT2"][2], 5)
    push!(week6[7], dict["DeadliftT2"][1], dict["DeadliftT2"][2], 6)
    push!(week7[7], dict["DeadliftT2"][1], dict["DeadliftT2"][2], 7)
    push!(week8[7], dict["DeadliftT2"][1], dict["DeadliftT2"][2], 8)
    push!(week9[7], dict["DeadliftT2"][1], dict["DeadliftT2"][2], 9)

    push!(week1[7], dict["HipThrustT1"][1], dict["HipThrustT1"][2], 1)
    push!(week2[7], dict["HipThrustT1"][1], dict["HipThrustT1"][2], 2)
    push!(week3[7], dict["HipThrustT1"][1], dict["HipThrustT1"][2], 3)
    push!(week4[7], dict["HipThrustT1"][1], dict["HipThrustT1"][2], 4)
    push!(week5[7], dict["HipThrustT1"][1], dict["HipThrustT1"][2], 5)
    push!(week6[7], dict["HipThrustT1"][1], dict["HipThrustT1"][2], 6)
    push!(week7[7], dict["HipThrustT1"][1], dict["HipThrustT1"][2], 7)
    push!(week8[7], dict["HipThrustT1"][1], dict["HipThrustT1"][2], 8)
    push!(week9[7], dict["HipThrustT1"][1], dict["HipThrustT1"][2], 9)

    push!(week1[7], dict["SideRaise"][1], dict["SideRaise"][2], 2)
    push!(week2[7], dict["SideRaise"][1], dict["SideRaise"][2], 3)
    push!(week3[7], dict["SideRaise"][1], dict["SideRaise"][2], 4)
    push!(week4[7], dict["SideRaise"][1], dict["SideRaise"][2], 5)
    push!(week5[7], dict["SideRaise"][1], dict["SideRaise"][2], 6)
    push!(week6[7], dict["SideRaise"][1], dict["SideRaise"][2], 7)
    push!(week7[7], dict["SideRaise"][1], dict["SideRaise"][2], 8)
    push!(week8[7], dict["SideRaise"][1], dict["SideRaise"][2], 9)

    push!(week1[7], dict["InclineDBBench"][1], dict["InclineDBBench"][2], 1)
    push!(week2[7], dict["InclineDBBench"][1], dict["InclineDBBench"][2], 2)
    push!(week3[7], dict["InclineDBBench"][1], dict["InclineDBBench"][2], 3)
    push!(week4[7], dict["InclineDBBench"][1], dict["InclineDBBench"][2], 4)
    push!(week5[7], dict["InclineDBBench"][1], dict["InclineDBBench"][2], 5)
    push!(week6[7], dict["InclineDBBench"][1], dict["InclineDBBench"][2], 6)
    push!(week7[7], dict["InclineDBBench"][1], dict["InclineDBBench"][2], 7)
    push!(week8[7], dict["InclineDBBench"][1], dict["InclineDBBench"][2], 8)
    push!(week9[7], dict["InclineDBBench"][1], dict["InclineDBBench"][2], 9)

    push!(week1[7], dict["BulgarianSplit"][1], dict["BulgarianSplit"][2], 2)
    push!(week2[7], dict["BulgarianSplit"][1], dict["BulgarianSplit"][2], 3)
    push!(week3[7], dict["BulgarianSplit"][1], dict["BulgarianSplit"][2], 4)
    push!(week4[7], dict["BulgarianSplit"][1], dict["BulgarianSplit"][2], 5)
    push!(week5[7], dict["BulgarianSplit"][1], dict["BulgarianSplit"][2], 6)
    push!(week6[7], dict["BulgarianSplit"][1], dict["BulgarianSplit"][2], 7)
    push!(week7[7], dict["BulgarianSplit"][1], dict["BulgarianSplit"][2], 8)
    push!(week8[7], dict["BulgarianSplit"][1], dict["BulgarianSplit"][2], 9)

    push!(week1[7], dict["PlatePinch"][1], dict["PlatePinch"][2], 2)
    push!(week2[7], dict["PlatePinch"][1], dict["PlatePinch"][2], 3)
    push!(week3[7], dict["PlatePinch"][1], dict["PlatePinch"][2], 4)
    push!(week4[7], dict["PlatePinch"][1], dict["PlatePinch"][2], 5)
    push!(week5[7], dict["PlatePinch"][1], dict["PlatePinch"][2], 6)
    push!(week6[7], dict["PlatePinch"][1], dict["PlatePinch"][2], 7)
    push!(week7[7], dict["PlatePinch"][1], dict["PlatePinch"][2], 8)
    push!(week8[7], dict["PlatePinch"][1], dict["PlatePinch"][2], 9)

    push!(week1[7], dict["ZottmanCurl"][1], dict["ZottmanCurl"][2], 2)
    push!(week2[7], dict["ZottmanCurl"][1], dict["ZottmanCurl"][2], 3)
    push!(week3[7], dict["ZottmanCurl"][1], dict["ZottmanCurl"][2], 4)
    push!(week4[7], dict["ZottmanCurl"][1], dict["ZottmanCurl"][2], 5)
    push!(week5[7], dict["ZottmanCurl"][1], dict["ZottmanCurl"][2], 6)
    push!(week6[7], dict["ZottmanCurl"][1], dict["ZottmanCurl"][2], 7)
    push!(week7[7], dict["ZottmanCurl"][1], dict["ZottmanCurl"][2], 8)
    push!(week8[7], dict["ZottmanCurl"][1], dict["ZottmanCurl"][2], 9)

    push!(week1[7], dict["ParaCrunch"][1], dict["ParaCrunch"][2], 2)
    push!(week2[7], dict["ParaCrunch"][1], dict["ParaCrunch"][2], 3)
    push!(week3[7], dict["ParaCrunch"][1], dict["ParaCrunch"][2], 4)
    push!(week4[7], dict["ParaCrunch"][1], dict["ParaCrunch"][2], 5)
    push!(week5[7], dict["ParaCrunch"][1], dict["ParaCrunch"][2], 6)
    push!(week6[7], dict["ParaCrunch"][1], dict["ParaCrunch"][2], 7)
    push!(week7[7], dict["ParaCrunch"][1], dict["ParaCrunch"][2], 8)
    push!(week8[7], dict["ParaCrunch"][1], dict["ParaCrunch"][2], 9)

    return vcat(week1, week2, week3, week4, week5, week6, week7, week8, week9)
end

days = makeDays(typeProgramme, exerProg)
push!(
    Lifting_Programmes,
    nameProgramme => Programme(typeProgramme, nameProgramme, exerProg, days),
)

max_stars_repo_name: dcelisgarza/Lifting


text:
# https://juliaci.github.io/PkgBenchmark.jl/stable/run_benchmarks/
# https://github.com/JuliaCI/PkgBenchmark.jl/blob/master/benchmark/benchmarks.jl
using BenchmarkTools
using Distributions
using HMMBase
using PyCall

# Helpers

# Recursively create benchmark groups
function mkgrp!(base, path)
    path = map(string, path)
    curr = base
    for k in path
        if !(k in keys(curr))
            curr[k] = BenchmarkGroup()
        end
        curr = curr[k]
    end
    curr
end

function mkbench!(f, base, path)
    grp = mkgrp!(base, path[1:end-1])
    grp[path[end]] = f()
end

# Dataset

function randhmm(K)
    A = randtransmat(K)
    B = [Normal(rand() * 100, rand() * 10) for _ in 1:K]
    HMM(A, B)
end

HMMs = Dict([K => randhmm(K) for K in 2:2:10])
Ys = Dict([K => rand(HMMs[K], 5000) for K in keys(HMMs)])

# Suite
const SUITE = BenchmarkGroup()

# HMMBase

for f in (forward, backward, viterbi), logl in [true, false], (K, hmm) in HMMs
    L = likelihoods(hmm, Ys[K], logl = logl)
    mkbench!(SUITE, ["hmmbase", f, logl, K]) do
        @benchmarkable ($f)($hmm.a, $hmm.A, $L, logl = $logl)
    end
end

# pyhsmm

pyhsmm = pyimport("pyhsmm")
pyhsmmi = pyimport("pyhsmm.internals.hmm_messages_interface")

# fs = [
#     :_messages_forwards_normalized,
#     :_messages_backwards_normalized
# ]

fs = [
    pyhsmm.internals.hmm_states.HMMStatesPython._messages_forwards_normalized,
    pyhsmm.internals.hmm_states.HMMStatesPython._messages_backwards_normalized
]

for f in fs, (K, hmm) in HMMs
    L = likelihoods(hmm, Ys[K], logl = true)
    mkbench!(SUITE, ["pyhsmm", f, K]) do
        @benchmarkable ($f)($hmm.A, $hmm.a, $L)
    end
end

for (K, hmm) in HMMs
    L = likelihoods(hmm, Ys[K], logl = true)
    buf = zeros(Int32, size(L, 1))

    a = hmm.a
    A = PyReverseDims(permutedims(hmm.A))
    L = PyReverseDims(permutedims(L))

    mkbench!(SUITE, ["pyhsmm", "viterbi", K]) do
        @benchmarkable pyhsmmi.viterbi($A, $L, $a, $buf)
    end
end

# hmmlearn

hmml = pyimport("hmmlearn.hmm")

function tohmmlearn(hmm::AbstractHMM)
    K = size(hmm, 1)
    model = hmml.GaussianHMM(n_components = K, covariance_type = :diag)
    model.startprob_ = hmm.a
    model.transmat_ = hmm.A # TODO: PermuteDims ?
    model.means_ = [d.μ for d in hmm.B]
    model.covars_ = reshape([d.σ^2 for d in hmm.B], (1, 1, :))
    model
end

fs = [:_do_forward_pass, :_do_backward_pass, :_do_viterbi_pass]

for f in fs, (K, hmm) in HMMs
    model = tohmmlearn(hmm)
    L = likelihoods(hmm, Ys[K], logl = true)
    mkbench!(SUITE, ["hmmlearn", f, K]) do
        @benchmarkable $(model).$(f)($L)
    end
end

max_stars_repo_name: Shu-Lea-Lai/module-MultiHMMBase


text:
module Sandbox
using Preferences, Scratch, Artifacts, Tar, TOML, Libdl

import Base: run
export SandboxExecutor, DockerExecutor, UserNamespacesExecutor, SandboxConfig,
       preferred_executor, executor_available, probe_executor, run

# Include some utilities for things like file manipulation, uname() parsing, etc...
include("utils.jl")

"""
    SandboxExecutor

This represents the base type for all execution backends within this package.
Valid concrete subtypes must implement at least the following methods:

* `executor_available(::DataType{T})`: Checks whether executor type `T` is available
  on this system.  For example, `UserNamespacesExecutor`s are only available on
  Linux, and even then only on certain kernels.  Availablility checks may run a
  program to determine whether that executor is actually available.

* `build_executor_command(exe::T, config::SandboxConfig, cmd::Cmd)`: Builds the
  `Cmd` object that, when run, executes the user's desired command within the given
  sandbox.  The `config` object contains all necessary metadata such as shard
  mappings, environment variables, `stdin`/`stdout`/`stderr` redirection, etc...
"""
abstract type SandboxExecutor; end

include("SandboxConfig.jl")

# Load the Docker executor
include("Docker.jl")

# Load the UserNamespace executor
include("UserNamespaces.jl")

all_executors = Type{<:SandboxExecutor}[
    # We always prefer the UserNamespaces executor, if we can use it,
    # and the unprivileged one most of all.  Only after that do we try `docker`.
    UnprivilegedUserNamespacesExecutor,
    PrivilegedUserNamespacesExecutor,
    DockerExecutor,
]

function select_executor(verbose::Bool)
    # If we have a preference set, use that unconditionally.
    executor = @load_preference("executor")
    if executor !== nothing
        executor = lowercase(executor)
        if executor ∈ ("unprivilegedusernamespacesexecutor", "unprivileged", "userns")
            return UnprivilegedUserNamespacesExecutor
        elseif executor ∈ ("privilegedusernamespacesexecutor", "privileged")
            return PrivilegedUserNamespacesExecutor
        elseif executor ∈ ("dockerexecutor", "docker")
            return DockerExecutor
        end
    end

    # Otherwise, just try them in priority order
    for executor in all_executors
        if executor_available(executor; verbose)
            return executor
        end
    end
    error("Could not find any available executors!")
end

_preferred_executor = nothing
function preferred_executor(;verbose::Bool = false)
    # If we've already asked this question, return the old answer
    global _preferred_executor
    if _preferred_executor === nothing
        _preferred_executor = select_executor(verbose)
    end
    return _preferred_executor()
end

# Helper function for warning about privileged execution trying to invoke `sudo`
prompted_userns_privileged = false
function warn_priviledged(::PrivilegedUserNamespacesExecutor)
    global prompted_userns_privileged
    if !prompted_userns_privileged
        @info("Running privileged container via `sudo`, may ask for your password:")
        prompted_userns_privileged = true
    end
end
warn_priviledged(::SandboxExecutor) = nothing

function run(exe::SandboxExecutor, config::SandboxConfig, user_cmd::Cmd)
    cmd = pipeline(build_executor_command(exe, config, user_cmd); config.stdin, config.stdout, config.stderr)
    if config.verbose
        @info("Running sandboxed command", cmd.exec)
    end
    warn_priviledged(exe)
    return success(run(cmd))
end

function probe_executor(executor::SandboxExecutor; verbose::Bool = false, test_read_only_map=false, test_read_write_map=false)
    mktempdir() do tmpdir
        read_only_maps = Dict{String,String}(
            "/" => artifact"AlpineRootfs",
        )
        read_write_maps = Dict{String,String}()

        # The simplest test is to see if stdout capturing works
        inner_cmd = "echo 'hello julia'"

        # Test whether we can use read-only mappings
        if test_read_only_map
            map_dir = joinpath(tmpdir, "read_only_map")
            mkdir(map_dir)
            open(joinpath(map_dir, "foo"), write=true) do io
                write(io, "read-only mapping successful")
            end

            # Mount that directory in /read_only
            read_only_maps["/read_only"] = map_dir

            # Read from the foo file
            inner_cmd = "$(inner_cmd) && cat /read_only/foo"
        end

        # Test whether we can use read-write mappings
        if test_read_write_map
            workspace_dir = joinpath(tmpdir, "read_write_map")
            mkdir(workspace_dir)

            # Mount that directory in /read_write
            read_write_maps["/read_write"] = workspace_dir

            # Write to the foo file
            inner_cmd = "$(inner_cmd) && echo read-write mapping successful >> /read_write/foo"
        end

        cmd_stdout = IOBuffer()
        cmd_stderr = IOBuffer()
        config = SandboxConfig(
            read_only_maps,
            read_write_maps,
            Dict("PATH" => "/bin:/usr/bin");
            stdout=cmd_stdout,
            stderr=cmd_stderr,
        )
        user_cmd = `/bin/sh -c "$(inner_cmd)"`

        if verbose
            tests = String[]
            if test_read_only_map
                push!(tests, "read-only")
            end
            if test_read_write_map
                push!(tests, "read-write")
            end
            @info("Testing $(executor) ($(join(tests, ", ")))")
        end

        # Command should execute successfully
        if !run(executor, config, user_cmd)
            if verbose
                cmd_stdout = String(take!(cmd_stdout))
                cmd_stderr = String(take!(cmd_stderr))
                @warn("Unable to run `sandbox` itself", cmd_stdout, cmd_stderr)
            end
            return false
        end

        # No stderr output (unless we're running in verbose mode)
        if !verbose
            stderr_output = String(take!(cmd_stderr))
            if !isempty(stderr_output)
                @warn(" -> Non-empty stderr output", stderr_output)
                return false
            end
        end

        # stdout shuold contain "hello julia" as its own line
        cmd_stdout = String(take!(cmd_stdout))
        stdout_lines = split(cmd_stdout, "\n")
        if !("hello julia" in stdout_lines)
            if verbose
                @warn(" -> Basic stdout sentinel missing!", stdout_lines)
            end
            return false
        end

        if test_read_only_map
            if !("read-only mapping successful" in stdout_lines)
                if verbose
                    @warn(" -> Read-only mapping sentinel missing!")
                end
                return false
            end
        end

        if test_read_write_map
            foo_file = joinpath(joinpath(tmpdir, "read_write_map", "foo"))
            if !isfile(foo_file)
                if verbose
                    @warn(" -> Read-write mapping sentinel file missing!")
                end
                return false
            end

            foo_file_contents = String(read(foo_file))
            if foo_file_contents != "read-write mapping successful\n"
                if verbose
                    @warn(" -> Read-write mapping data corrupted", foo_file_contents)
                end
                return false
            end
        end
        return true
    end
end

# Convenience function for other users who want to do some testing
alpine_rootfs() = artifact"AlpineRootfs"

end # module