text:
function makeAncillaryIs(A::fPEPS, L::Environments, R::Environments, col::Int)
    Ny, Nx         = size(A)
    is_cu          = is_gpu(A)
    dummy          = is_cu    ? cuITensor(1.0)  : ITensor(1.0)
    col_site_inds  = [firstind(x, "Site") for x in A[:, col]]
    ops            = map(x -> spinI(x; is_gpu=is_cu), col_site_inds)
    ancIs          = Vector{ITensor}(undef, Ny)
    ancI           = copy(dummy)
    for row in reverse(1:Ny)
        ancI      *= L.I[row] * A[row, col]
        ancI      *= ops[row] * dag(prime(A[row, col]))
        ancI      *= R.I[row]
        ancIs[row] = copy(ancI)
    end
    reverse!(ancIs)
    return ancIs
end

function makeAncillaryIsHorizontal(A::fPEPS, Aj::Vector{ITensor}, L::Environments, R::Environments, col::Int, ncol::Int)
    Ny, Nx         = size(A)
    is_cu          = is_gpu(A)
    dummy          = is_cu    ? cuITensor(1.0)  : ITensor(1.0)
    col_site_inds  = [firstind(x, "Site") for x in A[:, col]]
    ncol_site_inds = [firstind(x, "Site") for x in A[:, ncol]]
    ops            = map(x -> spinI(x; is_gpu=is_cu), col_site_inds) .* map(x -> spinI(x; is_gpu=is_cu), ncol_site_inds)
    ancIs          = Vector{ITensor}(undef, Ny)
    ancI           = copy(dummy)
    for row in reverse(1:Ny)
        ancI      *= L.I[row]
        ancI      *= Aj[row]
        ancI      *= ops[row] * dag(prime(Aj[row]))
        ancI      *= R.I[row]
        ancIs[row] = copy(ancI)
    end
    reverse!(ancIs)
    return ancIs
end

function updateAncillaryIs(A::fPEPS, Ibelow::Vector{ITensor}, L::Environments, R::Environments, row::Int, col::Int )
    Ny, Nx  = size(A)
    is_cu   = is_gpu(A)
    dummy   = is_cu   ? cuITensor(1.0)  : ITensor(1.0)
    op      = spinI(firstind(A[row, col], "Site"); is_gpu=is_cu)
    AA      = row > 1 ? copy(Ibelow[row - 1]) : dummy
    AA     *= L.I[row] * A[row, col]
    AA     *= op * dag(prime(A[row, col]))
    AA     *= R.I[row]
    Ibelow[row] = AA
    return Ibelow
end

function updateAncillaryIsHorizontal(A::fPEPS, Aj::Vector{ITensor}, Ibelow::Vector{ITensor}, L::Environments, R::Environments, row::Int, col::Int, ncol::Int )
    Ny, Nx  = size(A)
    is_cu   = is_gpu(A)
    dummy   = is_cu   ? cuITensor(1.0)  : ITensor(1.0)
    op      = spinI(firstind(A[row, col], "Site"); is_gpu=is_cu)
    op     *= spinI(firstind(A[row, ncol], "Site"); is_gpu=is_cu)
    AA      = row > 1 ? copy(Ibelow[row - 1]) : dummy
    AA     *= L.I[row] * Aj[row]
    AA     *= op * dag(prime(Aj[row]))
    AA     *= R.I[row]
    Ibelow[row] = AA
    return Ibelow
end

function makeAncillaryFs(A::fPEPS, L::Environments, R::Environments, H, col::Int)
    Ny, Nx   = size(A)
    is_cu    = is_gpu(A)
    dummy    = is_cu    ? cuITensor(1.0)  : ITensor(1.0)
    Fabove   = [Vector{ITensor}(undef, Ny) for opcode in 1:length(H)]
    col_site_inds = [firstind(x, "Site") for x in A[:, col]]
    for opcode in 1:length(H)
        op_row      = H[opcode].sites[1][1]
        ops         = map(x -> spinI(x; is_gpu=is_cu), col_site_inds)
        ops[op_row] = replaceind!(copy(H[opcode].ops[1]), H[opcode].site_ind, col_site_inds[op_row])
        ops[op_row] = replaceind!(copy(ops[op_row]), H[opcode].site_ind', col_site_inds[op_row]')
        ancFs       = Vector{ITensor}(undef, Ny)
        ancF        = copy(dummy)
        for row in reverse(1:Ny)
            ancF   *= L.I[row] * A[row, col]
            ancF   *= ops[row] * dag(prime(A[row, col]))
            ancF   *= R.I[row]
            ancFs[row] = copy(ancF)
        end
        Fabove[opcode] .= reverse(ancFs)
    end
    return Fabove
end

function makeAncillaryFsHorizontal(A::fPEPS, Aj::Vector{ITensor}, L::Environments, R::Environments, H, col::Int, ncol::Int)
    Ny, Nx   = size(A)
    is_cu    = is_gpu(A)
    dummy    = is_cu    ? cuITensor(1.0)  : ITensor(1.0)
    Fabove   = [Vector{ITensor}(undef, Ny) for opcode in 1:length(H)]
    col_site_inds = [firstind(x, "Site") for x in A[:, col]]
    ncol_site_inds = [firstind(x, "Site") for x in A[:, ncol]]
    for opcode in 1:length(H)
        op_row      = H[opcode].sites[1][1]
        op_col      = H[opcode].sites[1][2]
        ops         = map(x -> spinI(x; is_gpu=is_cu), col_site_inds) .* map(x -> spinI(x; is_gpu=is_cu), ncol_site_inds)
        if op_col == col
            op_a = replaceind!(copy(H[opcode].ops[1]), H[opcode].site_ind, col_site_inds[op_row])
            op_a = replaceind!(op_a, H[opcode].site_ind', col_site_inds[op_row]')
            ops[op_row] = op_a * spinI(ncol_site_inds[op_row]; is_gpu=is_cu)
        else
            op_b = replaceind!(copy(H[opcode].ops[1]), H[opcode].site_ind, ncol_site_inds[op_row])
            op_b = replaceind!(op_b, H[opcode].site_ind', ncol_site_inds[op_row]')
            ops[op_row] = op_b * spinI(col_site_inds[op_row]; is_gpu=is_cu)
        end
        ancFs       = Vector{ITensor}(undef, Ny)
        ancF        = copy(dummy)
        for row in reverse(1:Ny)
            ancF   *= L.I[row] * Aj[row]
            ancF   *= ops[row] * dag(prime(Aj[row]))
            ancF   *= R.I[row]
            ancFs[row] = copy(ancF)
        end
        Fabove[opcode] .= reverse(ancFs)
    end
    return Fabove
end

function updateAncillaryFs(A::fPEPS, Fbelow, Ibelow::Vector{ITensor}, L::Environments, R::Environments, H, row::Int, col::Int)
    Ny, Nx   = size(A)
    is_cu    = is_gpu(A)
    dummy    = is_cu    ? cuITensor(1.0)  : ITensor(1.0)
    col_site_inds = [firstind(x, "Site") for x in A[:, col]]
    for opcode in 1:length(H)
        op_row      = H[opcode].sites[1][1]
        ops         = ITensor[spinI(spin_ind; is_gpu=is_cu) for spin_ind in col_site_inds]
        if op_row == row
            ops[op_row] = replaceind!(copy(H[opcode].ops[1]), H[opcode].site_ind, col_site_inds[op_row])
            ops[op_row] = replaceind!(ops[op_row], H[opcode].site_ind', col_site_inds[op_row]')
        end
        ancF  = length(Fbelow[opcode]) > 0 ? copy(Fbelow[opcode][end]) : copy(dummy)
        ancF *= L.I[row]
        ancF *= A[row, col]
        ancF *= ops[row]
        ancF *= prime(dag(A[row, col]))
        ancF *= R.I[row]
        push!(Fbelow[opcode], ancF)
    end
    return Fbelow
end

function updateAncillaryFsHorizontal(A::fPEPS, Aj::Vector{ITensor}, Fbelow, Ibelow::Vector{ITensor}, L::Environments, R::Environments, H, row::Int, col::Int, ncol::Int)
    Ny, Nx   = size(A)
    is_cu    = is_gpu(A)
    dummy    = is_cu    ? cuITensor(1.0)  : ITensor(1.0)
    col_site_inds = [firstind(x, "Site") for x in A[:, col]]
    ncol_site_inds = [firstind(x, "Site") for x in A[:, ncol]]
    for opcode in 1:length(H)
        op_row      = H[opcode].sites[1][1]
        op_col      = H[opcode].sites[1][2]
        ops         = map(x -> spinI(x; is_gpu=is_cu), col_site_inds) .* map(x -> spinI(x; is_gpu=is_cu), ncol_site_inds)
        if op_row == row
            if op_col == col
                op_a = replaceind!(copy(H[opcode].ops[1]), H[opcode].site_ind, col_site_inds[op_row])
                op_a = replaceind!(op_a, H[opcode].site_ind', col_site_inds[op_row]')
                ops[op_row] = op_a*spinI(ncol_site_inds[op_row]; is_gpu=is_cu)
            else
                op_b = replaceind!(copy(H[opcode].ops[1]), H[opcode].site_ind, ncol_site_inds[op_row])
                op_b = replaceind!(op_b, H[opcode].site_ind', ncol_site_inds[op_row]')
                ops[op_row] = op_b*spinI(col_site_inds[op_row]; is_gpu=is_cu)
            end
        end
        ancF  = length(Fbelow[opcode]) > 0 ? copy(Fbelow[opcode][end]) : copy(dummy)
        ancF *= L.I[row]
        ancF *= Aj[row]
        ancF *= ops[row]
        ancF *= prime(dag(Aj[row]))
        ancF *= R.I[row]
        push!(Fbelow[opcode], ancF)
    end
    return Fbelow
end

function makeAncillaryIntsHorizontal(A::fPEPS, Aj::Vector{ITensor}, L::Environments, R::Environments, H, col::Int, ncol::Int)
    Ny, Nx   = size(A)
    is_cu    = is_gpu(A)
    dummy    = is_cu    ? cuITensor(1.0)  : ITensor(1.0)
    ITabove   = [Vector{ITensor}(undef, Ny) for opcode in 1:length(H)]
    col_site_inds = [firstind(x, "Site") for x in A[:, col]]
    ncol_site_inds = [firstind(x, "Site") for x in A[:, ncol]]
    for opcode in 1:length(H)
        op_row      = H[opcode].sites[1][1]
        op_col_a    = H[opcode].sites[1][2]
        op_col_b    = H[opcode].sites[2][2]
        ops         = map(x -> spinI(x; is_gpu=is_cu), col_site_inds) .* map(x -> spinI(x; is_gpu=is_cu), ncol_site_inds)
        if op_col_a == col
            op_a = replaceind!(copy(H[opcode].ops[1]), H[opcode].site_ind, col_site_inds[op_row])
            op_a = replaceind!(op_a, H[opcode].site_ind', col_site_inds[op_row]')
            op_b = replaceind!(copy(H[opcode].ops[2]), H[opcode].site_ind, ncol_site_inds[op_row])
            op_b = replaceind!(op_b, H[opcode].site_ind', ncol_site_inds[op_row]')
            ops[op_row] = op_a*op_b
        else
            op_a = replaceind!(copy(H[opcode].ops[1]), H[opcode].site_ind, ncol_site_inds[op_row])
            op_a = replaceind!(op_a, H[opcode].site_ind', ncol_site_inds[op_row]')
            op_b = replaceind!(copy(H[opcode].ops[2]), H[opcode].site_ind, col_site_inds[op_row])
            op_b = replaceind!(op_b, H[opcode].site_ind', col_site_inds[op_row]')
            ops[op_row] = op_a*op_b
        end
        ancITs       = Vector{ITensor}(undef, Ny)
        ancIT        = copy(dummy)
        for row in reverse(1:Ny)
            ancIT   *= L.I[row] * Aj[row]
            ancIT   *= ops[row] * dag(prime(Aj[row]))
            ancIT   *= R.I[row]
            ancITs[row] = copy(ancIT)
        end
        ITabove[opcode] .= reverse(ancITs)
    end
    return ITabove
end

function updateAncillaryIntsHorizontal(A::fPEPS, Aj::Vector{ITensor}, ITbelow, Ibelow::Vector{ITensor}, L::Environments, R::Environments, H, row::Int, col::Int, ncol::Int)
    Ny, Nx   = size(A)
    is_cu    = is_gpu(A)
    dummy    = is_cu    ? cuITensor(1.0)  : ITensor(1.0)
    col_site_inds = [firstind(x, "Site") for x in A[:, col]]
    ncol_site_inds = [firstind(x, "Site") for x in A[:, ncol]]
    for opcode in 1:length(H)
        op_row      = H[opcode].sites[1][1]
        op_col_a    = H[opcode].sites[1][2]
        op_col_b    = H[opcode].sites[2][2]
        ops         = map(x -> spinI(x; is_gpu=is_cu), col_site_inds) .* map(x -> spinI(x; is_gpu=is_cu), ncol_site_inds)
        if op_row == row
            if op_col_a == col
                op_a = replaceind!(copy(H[opcode].ops[1]), H[opcode].site_ind, col_site_inds[op_row])
                op_a = replaceind!(op_a, H[opcode].site_ind', col_site_inds[op_row]')
                op_b = replaceind!(copy(H[opcode].ops[2]), H[opcode].site_ind, ncol_site_inds[op_row])
                op_b = replaceind!(op_b, H[opcode].site_ind', ncol_site_inds[op_row]')
                ops[op_row] = op_a*op_b
            else
                op_a = replaceind!(copy(H[opcode].ops[1]), H[opcode].site_ind, ncol_site_inds[op_row])
                op_a = replaceind!(op_a, H[opcode].site_ind', ncol_site_inds[op_row]')
                op_b = replaceind!(copy(H[opcode].ops[2]), H[opcode].site_ind, col_site_inds[op_row])
                op_b = replaceind!(op_b, H[opcode].site_ind', col_site_inds[op_row]')
                ops[op_row] = op_a*op_b
            end
        end
        ancIT  = length(ITbelow[opcode]) > 0 ? copy(ITbelow[opcode][end]) : copy(dummy)
        ancIT *= L.I[row]
        ancIT *= Aj[row]
        ancIT *= ops[row]
        ancIT *= prime(dag(Aj[row]))
        ancIT *= R.I[row]
        push!(ITbelow[opcode], ancIT)
    end
    return ITbelow
end

function makeAncillaryVs(A::fPEPS, L::Environments, R::Environments, H, col::Int)
    Ny, Nx   = size(A)
    is_cu    = is_gpu(A)
    dummy    = is_cu    ? cuITensor(1.0)  : ITensor(1.0)
    col_site_inds = [firstind(x, "Site") for x in A[:, col]]
    Vabove   = [Vector{ITensor}(undef, Ny) for opcode in 1:length(H)]
    for opcode in 1:length(H)
        op_row_a      = H[opcode].sites[1][1]
        op_row_b      = H[opcode].sites[2][1]
        ops           = map(x->spinI(x; is_gpu=is_cu), col_site_inds)
        ops[op_row_a] = replaceind!(copy(H[opcode].ops[1]), H[opcode].site_ind, col_site_inds[op_row_a])
        ops[op_row_a] = replaceind!(ops[op_row_a], H[opcode].site_ind', col_site_inds[op_row_a]')
        ops[op_row_b] = replaceind!(copy(H[opcode].ops[2]), H[opcode].site_ind, col_site_inds[op_row_b])
        ops[op_row_b] = replaceind!(ops[op_row_b], H[opcode].site_ind', col_site_inds[op_row_b]')
        ancVs         = Vector{ITensor}(undef, Ny)
        ancV          = copy(dummy)
        for row in reverse(1:Ny)
            ancV   *= L.I[row] * A[row, col]
            ancV   *= ops[row] * dag(prime(A[row, col]))
            ancV   *= R.I[row]
            ancVs[row] = copy(ancV)
        end
        Vabove[opcode] = reverse(ancVs)
    end
    return Vabove
end

function makeAncillaryVsHorizontal(A::fPEPS, Aj::Vector{ITensor}, L::Environments, R::Environments, H, col::Int, ncol::Int)
    Ny, Nx   = size(A)
    is_cu    = is_gpu(A)
    dummy    = is_cu    ? cuITensor(1.0)  : ITensor(1.0)
    col_site_inds = [firstind(x, "Site") for x in A[:, col]]
    ncol_site_inds = [firstind(x, "Site") for x in A[:, ncol]]
    Vabove   = [Vector{ITensor}(undef, Ny) for opcode in 1:length(H)]
    for opcode in 1:length(H)
        op_row_a      = H[opcode].sites[1][1]
        op_row_b      = H[opcode].sites[2][1]
        op_col        = H[opcode].sites[1][2]
        ops           = map(x -> spinI(x; is_gpu=is_cu), col_site_inds) .* map(x -> spinI(x; is_gpu=is_cu), ncol_site_inds)
        if op_col == col
            op_a = replaceind!(copy(H[opcode].ops[1]), H[opcode].site_ind, col_site_inds[op_row_a])
            op_a = replaceind!(op_a, H[opcode].site_ind', col_site_inds[op_row_a]')
            op_b = replaceind!(copy(H[opcode].ops[2]), H[opcode].site_ind, col_site_inds[op_row_b])
            op_b = replaceind!(op_b, H[opcode].site_ind', col_site_inds[op_row_b]')
            ops[op_row_a] = op_a*spinI(ncol_site_inds[op_row_a]; is_gpu=is_cu)
            ops[op_row_b] = op_b*spinI(ncol_site_inds[op_row_b]; is_gpu=is_cu)
        else
            op_a = replaceind!(copy(H[opcode].ops[1]), H[opcode].site_ind, ncol_site_inds[op_row_a])
            op_a = replaceind!(op_a, H[opcode].site_ind', ncol_site_inds[op_row_a]')
            op_b = replaceind!(copy(H[opcode].ops[2]), H[opcode].site_ind, ncol_site_inds[op_row_b])
            op_b = replaceind!(op_b, H[opcode].site_ind', ncol_site_inds[op_row_b]')
            ops[op_row_a] = op_a*spinI(col_site_inds[op_row_a]; is_gpu=is_cu)
            ops[op_row_b] = op_b*spinI(col_site_inds[op_row_b]; is_gpu=is_cu)
        end
        ancVs         = Vector{ITensor}(undef, Ny)
        ancV          = copy(dummy)
        for row in reverse(1:Ny)
            ancV   *= L.I[row] * Aj[row]
            ancV   *= ops[row] * dag(prime(Aj[row]))
            ancV   *= R.I[row]
            ancVs[row] = copy(ancV)
        end
        Vabove[opcode] = reverse(ancVs)
    end
    return Vabove
end

function updateAncillaryVs(A::fPEPS, Vbelow, Ibelow::Vector{ITensor}, L::Environments, R::Environments, H, row::Int, col::Int)
    Ny, Nx   = size(A)
    is_cu    = is_gpu(A)
    dummy    = is_cu    ? cuITensor(1.0)  : ITensor(1.0)
    col_site_inds = [firstind(x, "Site") for x in A[:, col]]
    for opcode in 1:length(H)
        op_row_a      = H[opcode].sites[1][1]
        op_row_b      = H[opcode].sites[2][1]
        ops           = map(x->spinI(x; is_gpu=is_cu), col_site_inds)
        ops[op_row_a] = replaceind!(copy(H[opcode].ops[1]), H[opcode].site_ind, col_site_inds[op_row_a])
        ops[op_row_a] = replaceind!(ops[op_row_a], H[opcode].site_ind', col_site_inds[op_row_a]')
        ops[op_row_b] = replaceind!(copy(H[opcode].ops[2]), H[opcode].site_ind, col_site_inds[op_row_b])
        ops[op_row_b] = replaceind!(ops[op_row_b], H[opcode].site_ind', col_site_inds[op_row_b]')
        if op_row_b < row
            ancV  = copy(Vbelow[opcode][end])
            ancV *= L.I[row] * A[row, col]
            ancV *= ops[row] * dag(prime(A[row, col]))
            ancV *= R.I[row]
            push!(Vbelow[opcode], ancV)
        elseif op_row_b == row
            ancV  = op_row_a > 1 ? copy(Ibelow[op_row_a - 1]) : dummy
            for row_ in (op_row_a, op_row_b)
                ancV   *= L.I[row_] * A[row_, col]
                ancV   *= ops[row_] * dag(prime(A[row_, col]))
                ancV   *= R.I[row_]
            end
            push!(Vbelow[opcode], ancV)
        end
    end
    return Vbelow
end

function updateAncillaryVsHorizontal(A::fPEPS, Aj::Vector{ITensor}, Vbelow, Ibelow::Vector{ITensor}, L::Environments, R::Environments, H, row::Int, col::Int, ncol::Int)
    Ny, Nx   = size(A)
    is_cu    = is_gpu(A)
    dummy    = is_cu    ? cuITensor(1.0)  : ITensor(1.0)
    col_site_inds = [firstind(x, "Site") for x in A[:, col]]
    ncol_site_inds = [firstind(x, "Site") for x in A[:, ncol]]
    for opcode in 1:length(H)
        op_row_a      = H[opcode].sites[1][1]
        op_row_b      = H[opcode].sites[2][1]
        op_col        = H[opcode].sites[2][2]
        ops           = map(x->spinI(x; is_gpu=is_cu), col_site_inds) .* map(x->spinI(x; is_gpu=is_cu), ncol_site_inds)
        if op_col == col
            op_a = replaceind!(copy(H[opcode].ops[1]), H[opcode].site_ind, col_site_inds[op_row_a])
            op_a = replaceind!(op_a, H[opcode].site_ind', col_site_inds[op_row_a]')
            op_b = replaceind!(copy(H[opcode].ops[2]), H[opcode].site_ind, col_site_inds[op_row_b])
            op_b = replaceind!(op_b, H[opcode].site_ind', col_site_inds[op_row_b]')
            ops[op_row_a] = op_a * spinI(ncol_site_inds[op_row_a]; is_gpu=is_cu)
            ops[op_row_b] = op_b * spinI(ncol_site_inds[op_row_b]; is_gpu=is_cu)
        else
            op_a = replaceind!(copy(H[opcode].ops[1]), H[opcode].site_ind, ncol_site_inds[op_row_a])
            op_a = replaceind!(op_a, H[opcode].site_ind', ncol_site_inds[op_row_a]')
            op_b = replaceind!(copy(H[opcode].ops[2]), H[opcode].site_ind, ncol_site_inds[op_row_b])
            op_b = replaceind!(op_b, H[opcode].site_ind', ncol_site_inds[op_row_b]')
            ops[op_row_a] = op_a * spinI(col_site_inds[op_row_a]; is_gpu=is_cu)
            ops[op_row_b] = op_b * spinI(col_site_inds[op_row_b]; is_gpu=is_cu)
        end
        if op_row_b < row
            ancV  = copy(Vbelow[opcode][end])
            ancV *= L.I[row] * Aj[row]
            ancV *= ops[row] * dag(prime(Aj[row]))
            ancV *= R.I[row]
            push!(Vbelow[opcode], ancV)
        elseif op_row_b == row
            ancV  = op_row_a > 1 ? copy(Ibelow[op_row_a - 1]) : dummy
            for row_ in (op_row_a, op_row_b)
                ancV   *= L.I[row_] * Aj[row_]
                ancV   *= ops[row_] * dag(prime(Aj[row_]))
                ancV   *= R.I[row_]
            end
            push!(Vbelow[opcode], ancV)
        end
    end
    return Vbelow
end

function makeAncillarySide(A::fPEPS, EnvIP::Environments, EnvIdent::Environments, H, col::Int, side::Symbol)
    Ny, Nx   = size(A)
    is_cu    = is_gpu(A)
    col_site_inds = [firstind(x, "Site") for x in A[:, col]]
    Sabove   = [Vector{ITensor}(undef, Ny) for opcode in 1:length(H)]
    next_col = side == :left ? col + 1 : col - 1
    dummy    = is_cu    ? cuITensor(1.0)  : ITensor(1.0)
    for opcode in 1:length(H)
        op_row      = side == :left ? H[opcode].sites[2][1] : H[opcode].sites[1][1]
        ops         = map(x->spinI(x; is_gpu=is_cu), col_site_inds)
        this_op     = side == :left ? H[opcode].ops[2] : H[opcode].ops[1]
        ops[op_row] = replaceind!(copy(this_op), H[opcode].site_ind, col_site_inds[op_row])
        ops[op_row] = replaceind!(ops[op_row], H[opcode].site_ind', col_site_inds[op_row]')
        AAs         = Vector{ITensor}(undef, Ny)
        AA          = copy(dummy)
        for row in reverse(1:Ny)
            AA *= EnvIdent.I[row]
            AA *= A[row, col]
            AA *= ops[row]
            AA *= dag(prime(A[row, col]))
            AA *= EnvIP.InProgress[row, opcode]
            AAs[row] = copy(AA)
        end
        Sabove[opcode] = reverse(AAs)
    end
    return Sabove
end

function makeAncillarySideHorizontal(A::fPEPS, Aj::Vector{ITensor}, EnvIP::Environments, EnvIdent::Environments, H, col::Int, ncol::Int, side::Symbol)
    Ny, Nx   = size(A)
    is_cu    = is_gpu(A)
    col_site_inds = [firstind(x, "Site") for x in A[:, col]]
    ncol_site_inds = [firstind(x, "Site") for x in A[:, ncol]]
    Sabove   = [Vector{ITensor}(undef, Ny) for opcode in 1:length(H)]
    next_col = side == :left ? col + 1 : col - 1
    dummy    = is_cu    ? cuITensor(1.0)  : ITensor(1.0)
    for opcode in 1:length(H)
        op_row      = side == :left ? H[opcode].sites[2][1] : H[opcode].sites[1][1]
        op_col      = side == :left ? H[opcode].sites[2][2] : H[opcode].sites[1][2]
        ops         = map(x->spinI(x; is_gpu=is_cu), col_site_inds) .* map(x->spinI(x; is_gpu=is_cu), ncol_site_inds)
        this_op     = side == :left ? H[opcode].ops[2] : H[opcode].ops[1]
        if col == op_col
            op_a = replaceind!(copy(this_op), H[opcode].site_ind, col_site_inds[op_row])
            op_a = replaceind!(op_a, H[opcode].site_ind', col_site_inds[op_row]')
            ops[op_row] = op_a * spinI(ncol_site_inds[op_row]; is_gpu=is_cu)
        else
            op_b = replaceind!(copy(this_op), H[opcode].site_ind, ncol_site_inds[op_row])
            op_b = replaceind!(op_b, H[opcode].site_ind', ncol_site_inds[op_row]')
            ops[op_row] = op_b * spinI(col_site_inds[op_row]; is_gpu=is_cu)
        end
        AAs         = Vector{ITensor}(undef, Ny)
        AA          = copy(dummy)
        for row in reverse(1:Ny)
            AA *= EnvIdent.I[row]
            AA *= Aj[row]
            AA *= ops[row]
            AA *= dag(prime(Aj[row]))
            AA *= EnvIP.InProgress[row, opcode]
            AAs[row] = copy(AA)
        end
        Sabove[opcode] = reverse(AAs)
    end
    return Sabove
end

function updateAncillarySide(A::fPEPS, Sbelow, Ibelow::Vector{ITensor}, EnvIP::Environments, EnvIdent::Environments, H, row::Int, col::Int, side::Symbol)
    Ny, Nx   = size(A)
    is_cu    = is_gpu(A)
    col_site_inds = [firstind(x, "Site") for x in A[:, col]]
    next_col = side == :left ? col + 1 : col - 1
    prev_col = side == :left ? col + 1 : col - 1
    dummy    = is_cu    ? cuITensor(1.0)  : ITensor(1.0)
    for opcode in 1:length(H)
        op_row      = side == :left ? H[opcode].sites[2][1] : H[opcode].sites[1][1]
        ops         = ITensor[spinI(spin_ind; is_gpu=is_cu) for spin_ind in col_site_inds]
        this_op     = side == :left ? H[opcode].ops[2] : H[opcode].ops[1]
        ops[op_row] = replaceind!(copy(this_op), H[opcode].site_ind, col_site_inds[op_row])
        ops[op_row] = replaceind!(ops[op_row], H[opcode].site_ind', col_site_inds[op_row]')
        AA          = row >= 2 ? copy(Sbelow[opcode][row-1]) : copy(dummy)
        AA         *= EnvIdent.I[row]
        AA         *= A[row, col]
        AA         *= ops[row]
        AA         *= prime(dag(A[row, col]))
        AA         *= EnvIP.InProgress[row, opcode]
        push!(Sbelow[opcode], AA)
    end
    return Sbelow
end

function updateAncillarySideHorizontal(A::fPEPS, Aj::Vector{ITensor}, Sbelow, Ibelow::Vector{ITensor}, EnvIP::Environments, EnvIdent::Environments, H, row::Int, col::Int, ncol::Int, side::Symbol)
    Ny, Nx   = size(A)
    is_cu    = is_gpu(A)
    col_site_inds = [firstind(x, "Site") for x in A[:, col]]
    ncol_site_inds = [firstind(x, "Site") for x in A[:, ncol]]
    next_col = side == :left ? col + 1 : col - 1
    prev_col = side == :left ? col + 1 : col - 1
    dummy    = is_cu    ? cuITensor(1.0)  : ITensor(1.0)
    for opcode in 1:length(H)
        op_row      = side == :left ? H[opcode].sites[2][1] : H[opcode].sites[1][1]
        op_col      = side == :left ? H[opcode].sites[2][2] : H[opcode].sites[1][2]
        ops         = ITensor[spinI(spin_ind; is_gpu=is_cu) for spin_ind in col_site_inds] .* ITensor[spinI(spin_ind; is_gpu=is_cu) for spin_ind in ncol_site_inds]
        this_op     = side == :left ? H[opcode].ops[2] : H[opcode].ops[1]
        if op_col == col
            op_a = replaceind!(copy(this_op), H[opcode].site_ind, col_site_inds[op_row])
            op_a = replaceind!(op_a, H[opcode].site_ind', col_site_inds[op_row]')
            ops[op_row] = op_a * spinI(ncol_site_inds[op_row]; is_gpu=is_cu)
        else
            op_b = replaceind!(copy(this_op), H[opcode].site_ind, ncol_site_inds[op_row])
            op_b = replaceind!(op_b, H[opcode].site_ind', ncol_site_inds[op_row]')
            ops[op_row] = op_b * spinI(col_site_inds[op_row]; is_gpu=is_cu)
        end
        AA          = row >= 2 ? copy(Sbelow[opcode][row-1]) : copy(dummy)
        AA         *= EnvIdent.I[row]
        AA         *= Aj[row]
        AA         *= ops[row]
        AA         *= prime(dag(Aj[row]))
        AA         *= EnvIP.InProgress[row, opcode]
        push!(Sbelow[opcode], AA)
    end
    return Sbelow
end

function makeAncillarySideDiag(A::fPEPS, EnvIP::Environments, EnvIdent::Environments, H, col::Int, side::Symbol)
    Ny, Nx   = size(A)
    is_cu    = is_gpu(A)
    col_site_inds = [firstind(x, "Site") for x in A[:, col]]
    Sabove   = [Vector{ITensor}(undef, Ny) for opcode in 1:length(H)]
    next_col = side == :left ? col + 1 : col - 1
    dummy    = is_cu    ? cuITensor(1.0)  : ITensor(1.0)
    for opcode in 1:length(H)
        op_row      = side == :left ? H[opcode].sites[2][1] : H[opcode].sites[1][1]
        ops         = map(x->spinI(x; is_gpu=is_cu), col_site_inds)
        this_op     = side == :left ? H[opcode].ops[2] : H[opcode].ops[1]
        ops[op_row] = replaceind!(copy(this_op), H[opcode].site_ind, col_site_inds[op_row])
        ops[op_row] = replaceind!(ops[op_row], H[opcode].site_ind', col_site_inds[op_row]')
        AAs         = Vector{ITensor}(undef, Ny)
        AA          = copy(dummy)
        for row in reverse(1:Ny)
            AA *= EnvIdent.I[row]
            AA *= A[row, col]
            AA *= ops[row]
            AA *= dag(prime(A[row, col]))
            AA *= EnvIP.DiagInProgress[row, opcode]
            AAs[row] = copy(AA)
        end
        Sabove[opcode] = reverse(AAs)
    end
    return Sabove
end

function updateAncillarySideDiag(A::fPEPS, Dbelow, Ibelow::Vector{ITensor}, EnvIP::Environments, EnvIdent::Environments, H, row::Int, col::Int, side::Symbol)
    Ny, Nx   = size(A)
    is_cu    = is_gpu(A)
    next_col = side == :left ? col + 1 : col - 1
    prev_col = side == :left ? col + 1 : col - 1
    dummy    = is_cu    ? cuITensor(1.0)  : ITensor(1.0)
    col_site_inds = [firstind(x, "Site") for x in A[:, col]]
    for opcode in 1:length(H)
        op_row_a      = H[opcode].sites[1][1]
        op_row_b      = H[opcode].sites[2][1]
        op_col_a      = H[opcode].sites[1][2]
        op_col_b      = H[opcode].sites[2][2]
        ops           = map(x->spinI(x; is_gpu=is_cu), col_site_inds)
        if op_col_a == col
            ops[op_row_a] = replaceind!(copy(H[opcode].ops[1]), H[opcode].site_ind, col_site_inds[op_row_a])
            ops[op_row_a] = replaceind!(ops[op_row_a], H[opcode].site_ind', col_site_inds[op_row_a]')
        else
            ops[op_row_b] = replaceind!(copy(H[opcode].ops[2]), H[opcode].site_ind, col_site_inds[op_row_b])
            ops[op_row_b] = replaceind!(ops[op_row_b], H[opcode].site_ind', col_site_inds[op_row_b]')
        end
        ancD  = row >= 2 ? copy(Dbelow[opcode][row-1]) : copy(dummy)
        ancD *= EnvIdent.I[row]
        ancD *= A[row, col]
        ancD *= ops[row] * dag(prime(A[row, col]))
        ancD *= EnvIP.DiagInProgress[row, opcode]
        push!(Dbelow[opcode], ancD)
    end
    return Dbelow
end


text:
using BidirectionalMaps
using Test

@testset "main" begin
    U = ["A", "D", "B"]
    V = ["F", "G", "B"]

    b = ImmutableBimap{String,String}(U, V)
    @test b.left["A"] == "F"
    @test_throws MethodError b.left["A"] = "Z"
    @test b.left["D"] == "G"
    @test b.right["G"] == "D"
end

ext: jl


text:
# TODO
# remove Lx, Ly or save and use them consistently

struct TriangularLattice <: AbstractLattice
    Lx::Int
    Ly::Int
    sites::Int
    # row = up, down, left/right; col = siteidx
    neighs::Matrix{Int}
    # row (1) = up, down, left/right; cols (2,3) = cartesian siteidx
    neighs_cartesian::Array{Int, 3}
    # not really 2nd neighbors, but R_i + 2a_j
    ext_neighs::Matrix{Int}
    ext_neighs_cartesian::Array{Int, 3}

    lattice::Matrix{Int}
    isAsite::Vector{Bool}

    # for generic checkerboard decomposition
    n_bonds::Int
    bonds::Matrix{Int} # src, trg, type
end

function TriangularLattice(L::Int; Lx=L, Ly=L)
    sites = Lx * Ly
    lattice = convert(Array, reshape(1:sites, (Lx, Ly)))
    neighs, neighs_cartesian = build_neighbortable(TriangularLattice, lattice, Lx, Ly)
    ext_neighs, ext_neighs_cartesian = build_ext_neighbortable(
        TriangularLattice, lattice, Lx, Ly
    )

    # for generic checkerboard decomposition
    n_bonds = 6*sites
    bonds = zeros(n_bonds, 3)
    bondid = 1

    for src in lattice
        for trg in neighs[1:3, src]
            bonds[bondid, 1] = src
            bonds[bondid, 2] = trg
            bonds[bondid, 3] = 0
            bondid += 1
        end
        for trg in ext_neighs[1:3, src]
            bonds[bondid, 1] = src
            bonds[bondid, 2] = trg
            bonds[bondid, 3] = 0
            bondid += 1
        end
    end

    isAsite = [iseven(i) for i in 1:Lx, j in 1:Ly][:]

    return TriangularLattice(
        Lx, Ly, sites,
        neighs, neighs_cartesian,
        ext_neighs, ext_neighs_cartesian,
        lattice, isAsite, n_bonds, bonds
    )
end

function build_neighbortable(::Type{TriangularLattice}, lattice, Lx, Ly)
    up          = circshift(lattice, (-1,  0))
    upright     = circshift(lattice, (-1, -1))
    right       = circshift(lattice, ( 0, -1))
    down        = circshift(lattice, ( 1,  0))
    downleft    = circshift(lattice, ( 1,  1))
    left        = circshift(lattice, ( 0,  1))

    neighs = vcat(up[:]', upright[:]', right[:]', down[:]', downleft[:]', left[:]')

    neighs_cartesian = Array{Int, 3}(undef, 6, Lx, Ly)
    neighs_cartesian[1,:,:] = up
    neighs_cartesian[2,:,:] = upright
    neighs_cartesian[3,:,:] = right
    neighs_cartesian[4,:,:] = down
    neighs_cartesian[5,:,:] = downleft
    neighs_cartesian[6,:,:] = left

    return neighs, neighs_cartesian
end

function build_ext_neighbortable(::Type{TriangularLattice}, lattice, Lx, Ly)
    up          = circshift(lattice, (-2,  0))
    upright     = circshift(lattice, (-2, -2))
    right       = circshift(lattice, ( 0, -2))
    down        = circshift(lattice, ( 2,  0))
    downleft    = circshift(lattice, ( 2,  2))
    left        = circshift(lattice, ( 0,  2))

    ext_neighs = vcat(up[:]', upright[:]', right[:]', down[:]', downleft[:]', left[:]')

    ext_neighs_cartesian = Array{Int, 3}(undef, 6, Lx, Ly)
    ext_neighs_cartesian[1,:,:] = up
    ext_neighs_cartesian[2,:,:] = upright
    ext_neighs_cartesian[3,:,:] = right
    ext_neighs_cartesian[4,:,:] = down
    ext_neighs_cartesian[5,:,:] = downleft
    ext_neighs_cartesian[6,:,:] = left

    return ext_neighs, ext_neighs_cartesian
end

# Implement AbstractLattice interface: mandatory
@inline Base.length(l::TriangularLattice) = l.sites
@inline Base.size(l::TriangularLattice) = (l.Lx, l.Ly)

# Implement AbstractLattice interface: optional
@inline neighbors_lookup_table(l::TriangularLattice) = copy(l.neighs)

# HasNeighborsTable and HasBondsTable traits
has_neighbors_table(::TriangularLattice) = HasNeighborsTable()
has_bonds_table(::TriangularLattice) = HasBondsTable()

function positions(l::TriangularLattice)
    idxs = l.lattice |> CartesianIndices .|> Tuple .|> collect
    [[0.5, 0.8660254037844386] * idx[1] + [1, 0] * idx[2] for idx in idxs]
end
lattice_vectors(l::TriangularLattice) = [[0.5, 0.8660254037844386] * l.Lx, [l.Ly, 0]]


text:
import MathProgBase
const MPB = MathProgBase.SolverInterface
import SemidefiniteModels
const SDM = SemidefiniteModels

export CSDPMathProgModel, CSDPSolver

function checkoptions(d::Dict{Symbol, Any})
    for key in keys(d)
        if !(key in ALLOWED_OPTIONS)
            error("Option $key is not not a valid CSDP option. The valid options are $ALLOWED_OPTIONS.")
        end
    end
    return d
end

struct CSDPSolver <: MPB.AbstractMathProgSolver
    options::Dict{Symbol,Any}
end
CSDPSolver(; kwargs...) = CSDPSolver(checkoptions(Dict{Symbol,Any}(kwargs)))

mutable struct CSDPMathProgModel <: SDM.AbstractSDModel
    blockdims::Vector{CSDP_INT}
    b::Vector{Cdouble}
    entries::Vector{Tuple{CSDP_INT,CSDP_INT,CSDP_INT,CSDP_INT,Cdouble}}
    C::blockmatrix
    problem::Union{Nothing, LoadingProblem}
    X::blockmatrix
    y::Union{Nothing, Vector{Cdouble}}
    Z::blockmatrix
    status::CSDP_INT
    pobj::Cdouble
    dobj::Cdouble
    options::Dict{Symbol,Any}
    function CSDPMathProgModel(; kwargs...)
        new(CSDP_INT[], Cdouble[], Tuple{CSDP_INT,CSDP_INT,CSDP_INT,CSDP_INT,Cdouble}[],
            blockmatrix(), nothing, blockmatrix(), nothing, blockmatrix(),
            -1, 0.0, 0.0, checkoptions(Dict{Symbol, Any}(kwargs)))
    end
end
SDM.SDModel(s::CSDPSolver) = CSDPMathProgModel(; s.options...)
MPB.ConicModel(s::CSDPSolver) = SDM.SDtoConicBridge(SDM.SDModel(s))
MPB.LinearQuadraticModel(s::CSDPSolver) = MPB.ConicToLPQPBridge(MPB.ConicModel(s))

MPB.supportedcones(s::CSDPSolver) = [:Free,:Zero,:NonNeg,:NonPos,:SOC,:RSOC,:SDP]
function MPB.setvartype!(m::CSDPMathProgModel, vtype, blk, i, j)
    if vtype != :Cont
        error("Unsupported variable type $vtype by CSDP")
    end
end

function MPB.loadproblem!(m::CSDPMathProgModel, filename::AbstractString)
    if m.problem !== nothing
        if m.y !== nothing
            free_loaded_prob(m.problem, m.X, m.y, m.Z)
        end
        free_loading_prob(m.problem)
    end
    m.problem = nothing
    m.y = nothing
    if endswith(filename,".dat-s")
        m.problem = load_prob_from_file(filename, Ref(m.C))
    else
        error("unrecognized input format extension in $filename")
    end
end
#writeproblem(m, filename::String)
function MPB.loadproblem!(m::CSDPMathProgModel, blkdims::Vector{Int}, constr::Int)
    if m.problem !== nothing
        if m.y !== nothing
            free_loaded_prob(m.problem, m.X, m.y, m.Z)
        end
        free_loading_prob(m.problem)
    end
    m.problem = nothing
    m.y = nothing
    m.blockdims = blkdims
    m.b = zeros(Cdouble, constr)
    empty!(m.entries)
end

function SDM.setconstrB!(m::CSDPMathProgModel, val, constr::Integer)
    m.b[constr] = val
end
function SDM.setconstrentry!(m::CSDPMathProgModel, coef, constr::Integer, blk::Integer, i::Integer, j::Integer)
    push!(m.entries, (constr, blk, i, j, coef))
end
function SDM.setobjentry!(m::CSDPMathProgModel, coef, blk::Integer, i::Integer, j::Integer)
    push!(m.entries, (0, blk, i, j, coef))
end

function MPB.optimize!(m::CSDPMathProgModel)
    if m.problem === nothing
        # `m.problem` is not `nothing` if it was loaded from a file.
        m.C.nblocks = length(m.blockdims)
        num_entries = zeros(CSDP_INT, length(m.b), length(m.blockdims))
        for entry in m.entries
            if entry[1] > 0
                num_entries[entry[1], entry[2]] += 1
            end
        end
        m.problem = allocate_loading_prob(Ref(m.C), m.blockdims, length(m.b), num_entries, 3)
        for (i, x) in enumerate(m.b)
            setconstant(m.problem, i, x)
        end
        for entry in m.entries
            duplicate = addentry(m.problem, entry..., true)
            @assert !duplicate
        end
    end
    m.y = loaded_initsoln(m.problem, length(m.b), Ref(m.X), Ref(m.Z))
    m.status, m.pobj, m.dobj = loaded_sdp(m.problem, Ref(m.X), m.y, Ref(m.Z), m.options)
end

function MPB.status(m::CSDPMathProgModel)
    status = m.status
    if status == 0
        return :Optimal
    elseif status == 1
        return :Infeasible
    elseif status == 2
        return :Unbounded
    elseif status == 3
        return :Suboptimal
    elseif status == 4
        return :UserLimit
    elseif 5 <= status <= 9
        return :Error
    elseif status == -1
        return :Uninitialized
    else
        error("Internal library error: status=$status")
    end
end

function MPB.getobjval(m::CSDPMathProgModel)
    m.pobj
end
function MPB.getsolution(m::CSDPMathProgModel)
    m.X
end
function MPB.getdual(m::CSDPMathProgModel)
    m.y
end
function MPB.getvardual(m::CSDPMathProgModel)
    m.Z
end


text:
"""
SensitiveFloat objects for expectations involving r_s and c_s.

Args:
vs: A vector of variational parameters

Attributes:
Each matrix has one row for each color and a column for
star / galaxy.  Row 3 is the gamma distribute baseline brightness,
and all other rows are lognormal offsets.
- E_l_a: A B x Ia matrix of expectations and derivatives of
  color terms.  The rows are bands, and the columns
  are star / galaxy.
- E_ll_a: A B x Ia matrix of expectations and derivatives of
  squared color terms.  The rows are bands, and the columns
  are star / galaxy.
"""
immutable SourceBrightness{NumType <: Number}
    # [E[l|a=0], E[l]|a=1]]
    E_l_a::Matrix{SensitiveFloat{NumType}}

    # [E[l^2|a=0], E[l^2]|a=1]]
    E_ll_a::Matrix{SensitiveFloat{NumType}}
end

function SourceBrightness{NumType <: Number}(vs::Vector{NumType};
                                             calculate_gradient=true,
                                             calculate_hessian=true)
    r1 = vs[ids.r1]
    r2 = vs[ids.r2]
    c1 = vs[ids.c1]
    c2 = vs[ids.c2]

    # E_l_a has a row for each of the five colors and columns
    # for star / galaxy.
    E_l_a = Array(SensitiveFloat{NumType}, B, Ia)
    E_ll_a = Array(SensitiveFloat{NumType}, B, Ia)

    for i = 1:Ia
        ids_band_3 = Int[bids.r1, bids.r2]
        ids_color_1 = Int[bids.c1[1], bids.c2[1]]
        ids_color_2 = Int[bids.c1[2], bids.c2[2]]
        ids_color_3 = Int[bids.c1[3], bids.c2[3]]
        ids_color_4 = Int[bids.c1[4], bids.c2[4]]

        for b = 1:B
            E_l_a[b, i] = SensitiveFloat{NumType}(length(BrightnessParams), 1,
                                       calculate_gradient, calculate_hessian)
        end

        E_l_a[3, i].v[] = exp(r1[i] + 0.5 * r2[i])
        E_l_a[4, i].v[] = exp(c1[3, i] + .5 * c2[3, i])
        E_l_a[5, i].v[] = exp(c1[4, i] + .5 * c2[4, i])
        E_l_a[2, i].v[] = exp(-c1[2, i] + .5 * c2[2, i])
        E_l_a[1, i].v[] = exp(-c1[1, i] + .5 * c2[1, i])

        if calculate_gradient
            # band 3 is the reference band, relative to which the colors are
            # specified.
            # It is denoted r_s and has a lognormal expectation.
            E_l_a[3, i].d[bids.r1] = E_l_a[3, i].v[]
            E_l_a[3, i].d[bids.r2] = E_l_a[3, i].v[] * .5

            if calculate_hessian
                set_hess!(E_l_a[3, i], bids.r1, bids.r1, E_l_a[3, i].v[])
                set_hess!(E_l_a[3, i], bids.r1, bids.r2, E_l_a[3, i].v[] * 0.5)
                set_hess!(E_l_a[3, i], bids.r2, bids.r2, E_l_a[3, i].v[] * 0.25)
            end

            # The remaining indices involve c_s and have lognormal
            # expectations times E_c_3.

            # band 4 = band 3 * color 3.
            E_l_a[4, i].d[bids.c1[3]] = E_l_a[4, i].v[]
            E_l_a[4, i].d[bids.c2[3]] = E_l_a[4, i].v[] * .5
            if calculate_hessian
                set_hess!(E_l_a[4, i], bids.c1[3], bids.c1[3], E_l_a[4, i].v[])
                set_hess!(E_l_a[4, i], bids.c1[3], bids.c2[3], E_l_a[4, i].v[] * 0.5)
                set_hess!(E_l_a[4, i], bids.c2[3], bids.c2[3], E_l_a[4, i].v[] * 0.25)
            end
            multiply_sfs!(E_l_a[4, i], E_l_a[3, i])

            # Band 5 = band 4 * color 4.
            E_l_a[5, i].d[bids.c1[4]] = E_l_a[5, i].v[]
            E_l_a[5, i].d[bids.c2[4]] = E_l_a[5, i].v[] * .5
            if calculate_hessian
                set_hess!(E_l_a[5, i], bids.c1[4], bids.c1[4], E_l_a[5, i].v[])
                set_hess!(E_l_a[5, i], bids.c1[4], bids.c2[4], E_l_a[5, i].v[] * 0.5)
                set_hess!(E_l_a[5, i], bids.c2[4], bids.c2[4], E_l_a[5, i].v[] * 0.25)
            end
            multiply_sfs!(E_l_a[5, i], E_l_a[4, i])

            # Band 2 = band 3 * color 2.
            E_l_a[2, i].d[bids.c1[2]] = E_l_a[2, i].v[] * -1.
            E_l_a[2, i].d[bids.c2[2]] = E_l_a[2, i].v[] * .5
            if calculate_hessian
                set_hess!(E_l_a[2, i], bids.c1[2], bids.c1[2], E_l_a[2, i].v[])
                set_hess!(E_l_a[2, i], bids.c1[2], bids.c2[2], E_l_a[2, i].v[] * -0.5)
                set_hess!(E_l_a[2, i], bids.c2[2], bids.c2[2], E_l_a[2, i].v[] * 0.25)
            end
            multiply_sfs!(E_l_a[2, i], E_l_a[3, i])

            # Band 1 = band 2 * color 1.
            E_l_a[1, i].d[bids.c1[1]] = E_l_a[1, i].v[] * -1.
            E_l_a[1, i].d[bids.c2[1]] = E_l_a[1, i].v[] * .5
            if calculate_hessian
                set_hess!(E_l_a[1, i], bids.c1[1], bids.c1[1], E_l_a[1, i].v[])
                set_hess!(E_l_a[1, i], bids.c1[1], bids.c2[1], E_l_a[1, i].v[] * -0.5)
                set_hess!(E_l_a[1, i], bids.c2[1], bids.c2[1], E_l_a[1, i].v[] * 0.25)
            end
            multiply_sfs!(E_l_a[1, i], E_l_a[2, i])
        else
            # Simply update the values if not calculating derivatives.
            E_l_a[4, i].v[] *= E_l_a[3, i].v[]
            E_l_a[5, i].v[] *= E_l_a[4, i].v[]
            E_l_a[2, i].v[] *= E_l_a[3, i].v[]
            E_l_a[1, i].v[] *= E_l_a[2, i].v[]
        end # Derivs

        ################################
        # Squared terms.

        for b = 1:B
            E_ll_a[b, i] = SensitiveFloat{NumType}(length(BrightnessParams), 1,
                                           calculate_gradient, calculate_hessian)
        end

        E_ll_a[3, i].v[] = exp(2 * r1[i] + 2 * r2[i])
        E_ll_a[4, i].v[] = exp(2 * c1[3, i] + 2 * c2[3, i])
        E_ll_a[5, i].v[] = exp(2 * c1[4, i] + 2 * c2[4, i])
        E_ll_a[2, i].v[] = exp(-2 * c1[2, i] + 2 * c2[2, i])
        E_ll_a[1, i].v[] = exp(-2 * c1[1, i] + 2 * c2[1, i])

        if calculate_gradient
            # Band 3, the reference band.
            E_ll_a[3, i].d[bids.r1] = 2 * E_ll_a[3, i].v[]
            E_ll_a[3, i].d[bids.r2] = 2 * E_ll_a[3, i].v[]
            if calculate_hessian
                for hess_ids in [(bids.r1, bids.r1),
                                 (bids.r1, bids.r2),
                                 (bids.r2, bids.r2)]
                    set_hess!(E_ll_a[3, i], hess_ids..., 4.0 * E_ll_a[3, i].v[])
                end
            end

            # Band 4 = band 3 * color 3.
            E_ll_a[4, i].d[bids.c1[3]] = E_ll_a[4, i].v[] * 2.
            E_ll_a[4, i].d[bids.c2[3]] = E_ll_a[4, i].v[] * 2.
            if calculate_hessian
                for hess_ids in [(bids.c1[3], bids.c1[3]),
                                 (bids.c1[3], bids.c2[3]),
                                 (bids.c2[3], bids.c2[3])]
                    set_hess!(E_ll_a[4, i], hess_ids..., E_ll_a[4, i].v[] * 4.0)
                end
            end
            multiply_sfs!(E_ll_a[4, i], E_ll_a[3, i])

            # Band 5 = band 4 * color 4.
            tmp4 = exp(2 * c1[4, i] + 2 * c2[4, i])
            E_ll_a[5, i].d[bids.c1[4]] = E_ll_a[5, i].v[] * 2.
            E_ll_a[5, i].d[bids.c2[4]] = E_ll_a[5, i].v[] * 2.
            if calculate_hessian
                for hess_ids in [(bids.c1[4], bids.c1[4]),
                                 (bids.c1[4], bids.c2[4]),
                                 (bids.c2[4], bids.c2[4])]
                    set_hess!(E_ll_a[5, i], hess_ids..., E_ll_a[5, i].v[] * 4.0)
                end
            end
            multiply_sfs!(E_ll_a[5, i], E_ll_a[4, i])

            # Band 2 = band 3 * color 2
            tmp2 = exp(-2 * c1[2, i] + 2 * c2[2, i])
            E_ll_a[2, i].d[bids.c1[2]] = E_ll_a[2, i].v[] * -2.
            E_ll_a[2, i].d[bids.c2[2]] = E_ll_a[2, i].v[] * 2.
            if calculate_hessian
                for hess_ids in [(bids.c1[2], bids.c1[2]),
                                 (bids.c2[2], bids.c2[2])]
                    set_hess!(E_ll_a[2, i], hess_ids..., E_ll_a[2, i].v[] * 4.0)
                end
                set_hess!(E_ll_a[2, i], bids.c1[2], bids.c2[2],
                          E_ll_a[2, i].v[] * -4.0)
            end
            multiply_sfs!(E_ll_a[2, i], E_ll_a[3, i])

            # Band 1 = band 2 * color 1
            E_ll_a[1, i].d[bids.c1[1]] = E_ll_a[1, i].v[] * -2.
            E_ll_a[1, i].d[bids.c2[1]] = E_ll_a[1, i].v[] * 2.
            if calculate_hessian
                for hess_ids in [(bids.c1[1], bids.c1[1]),
                                 (bids.c2[1], bids.c2[1])]
                    set_hess!(E_ll_a[1, i], hess_ids..., E_ll_a[1, i].v[] * 4.0)
                end
                set_hess!(E_ll_a[1, i], bids.c1[1], bids.c2[1],
                          E_ll_a[1, i].v[] * -4.0)
            end
            multiply_sfs!(E_ll_a[1, i], E_ll_a[2, i])
        else
            # Simply update the values if not calculating derivatives.
            E_ll_a[4, i].v[] *= E_ll_a[3, i].v[]
            E_ll_a[5, i].v[] *= E_ll_a[4, i].v[]
            E_ll_a[2, i].v[] *= E_ll_a[3, i].v[]
            E_ll_a[1, i].v[] *= E_ll_a[2, i].v[]
        end # calculate_gradient
    end

    SourceBrightness(E_l_a, E_ll_a)
end

"""
A convenience function for getting only the brightness parameters
from model parameters.

Args:
  ea: Model parameters

Returns:
  An array of E_l_a and E_ll_a for each source.
"""
function get_brightness{NumType <: Number}(ea::ElboArgs{NumType})
    brightness = [SourceBrightness(ea.vp[s]) for s in ea.S]
    brightness_vals = [ Float64[b.E_l_a[i, j].v[] for
        i=1:size(b.E_l_a, 1), j=1:size(b.E_l_a, 2)] for b in brightness]
    brightness_squares = [ Float64[b.E_l_a[i, j].v[] for
        i=1:size(b.E_ll_a, 1), j=1:size(b.E_ll_a, 2)] for b in brightness]

    brightness_vals, brightness_squares
end

"""
Load the source brightnesses for these model params.  Each SourceBrightness
object has information for all bands and object types.

Returns:
  - An array of SourceBrightness objects for each object in 1:ea.S.  Only
    sources in ea.active_sources will have derivative information.
"""
function load_source_brightnesses{NumType <: Number}(
                    ea::ElboArgs{NumType};
                    calculate_gradient::Bool=true,
                    calculate_hessian::Bool=true)
    sbs = Array(SourceBrightness{NumType}, ea.S)

    for s in 1:ea.S
        this_deriv = (s in ea.active_sources) && calculate_gradient
        this_hess = (s in ea.active_sources) && calculate_hessian
        sbs[s] = SourceBrightness(ea.vp[s];
                                  calculate_gradient=this_deriv,
                                  calculate_hessian=this_hess)
    end

    sbs
end


text:
# ![invalid cover](nonexistence.jpg)
# <!-- ![not cover](../not_cover.png) -->
# ![valid cover](../logo.png)
# <!-- ![not cover](../not_cover.png) -->

# First valid image link is the cover


text:
function gitignore()
    return "# Temporary
*.jl.cov
*.DS_Store
*.tar.gz
*.jl.*.cov
*.jl.mem

# Data
*.vtk
*.dat
*.csv
*.vtu
*.pvtu
*.swp
*.png
*.nc

# System-specific files and directories generated by the BinaryProvider and BinDeps packages
# They contain absolute paths specific to the host computer, and so should not be committed
deps/deps.jl
deps/build.log
deps/downloads/
deps/usr/
deps/src/

# Build artifacts for creating documentation generated by the Documenter package
!docs/src/assets/*.png
docs/build/
docs/site/

# File generated by Pkg, the package manager, based on a corresponding Project.toml
# It records a fixed state of all packages used by the project. As such, it should not be
# committed for packages, but should be committed for applications that require a static
# environment.
Manifest.toml
"

end

max_stars_repo_name: charleskawczynski/PkgTemplates

max_issues_repo_name: charleskawczynski/PkgTemplates

max_forks_repo_name: charleskawczynski/PkgTemplates


text:
# This file investigates the effect of number of points on the 2D hidden integration error.

using FractalTools
using Makie

# Construct interpolation data
f(x, y) = [
    x^2 + y^2 + 1,
    x^y - y^2
]
vtx = [
    BigFloat.([0.0, 0.0]),
    BigFloat.([1.0, 0.0]),
    BigFloat.([0.5, 1.0])
    ]
freevar = 0.001
npts    = 50 : 5 : 200
ntpts   = 2 * npts[end]

# Compute errors
fval = 35 / 48
mse = map(npts) do npt
    @info npt
    # Construct interpolation data
    pts = getdata(f, vtx, npt)

    # Construct interpolant
    ival = integrate(pts, HInterp2D(fill(freevar, 2, 2)))

    # Compute error
    abs(fval - ival)
end

# Plot mse
fig = Figure()
ax = fig[1, 1] = Axis(fig, xlabel="Number of Points", ylabel="Error", title="2D Hidden Integration Error")
stem!(ax, npts, mse, color=:black)
save(joinpath(@__DIR__, "hinteg2d_mse.png"), fig)
display(fig)


text:
using DrWatson
@quickactivate "StatReth"

# %%
using StatsBase
using Distributions
using Plots
using StatsPlots
using Turing  # for MCMCChains

# %% 3.2
n = 1000
p_grid = range(0, 1, length = n)
prob_p = ones(n)
prob_data = @. pdf(Binomial(9, p_grid), 6)
posterior = prob_data .* prob_p
posterior ./= sum(posterior)

plot(p_grid, posterior)

# %% 3.3 - 3.5
weights = pweights(posterior)
samples = sample(p_grid, weights, 10_000)

scatter(samples)
density(samples)
# plot!(p_grid, posterior * n)

# %% 3.6, 3.7
sum(posterior[p_grid .< 0.5])

sum(samples .< 0.5) / 10_000    # either
mean(samples .< 0.5)            # or

# %% 3.8
mean(0.5 .< samples .< 0.75)    # either
mean(samples) do s              # or (this will be faster than above)
    0.5 < s < 0.75
end
mean(samples) do s              # or (can do multi-line; return variable on the last line)
    condition1 = 0.5 < s
    condition2 = s < 0.75
    both_conditions = condition1 & condition2
    both_conditions
end

# %% 3.9, 3.10
quantile(samples, 0.8)
quantile(samples, (0.1, 0.9))

# %% 3.11
p_grid = range(0, 1, length = n)
prior = ones(n)
likelihood = @. pdf(Binomial(3, p_grid), 3)
posterior = likelihood .* prior
posterior ./= sum(posterior)
samples = sample(p_grid, pweights(posterior), 10_000)

# %% 3.12, 3.13
quantile(samples, (0.25, 0.75))

chn = Chains(samples)
hpd(chn; alpha = 0.5)

# %% 3.14, 3.16, TODO: not sure about 3.15
p_grid[argmax(posterior)]

mean(samples)
median(samples)

# %% 3.17
sum(posterior .* abs.(0.5 .- p_grid))

# %% 3.18, 3.19
loss = map(p_grid) do d
    sum(posterior .* abs.(d .- p_grid))
end
plot(loss)

p_grid[argmin(loss)]


text:
# using TuringModels
using CSV, DataFrames, Query, Turing, Optim, StatsBase, StatsPlots

# define precis from StatisticalRethinking module
using NamedArrays
function precis(df::DataFrame; digits=3, depth=Inf, alpha=0.11)
  m = zeros(length(names(df)), 5)
  for (indx, col) in enumerate(names(df))
    m[indx, 1] = mean(df[:, col])
    m[indx, 2] = std(df[:, col])
    q = quantile(df[:, col], [alpha/2, 0.5, 1-alpha/2])
    m[indx, 3] = q[1]
    m[indx, 4] = q[2]
    m[indx, 5] = q[3]
  end
  NamedArray(
    round.(m, digits=digits),
    (names(df), ["mean", "sd", "5.5%", "50%", "94.5%"]),
    ("Rows", "Cols"))
end

# ### snippet 5.3

# TODO: use population standard deviation to be consistent with Rethinking?
scale(x) = (x .- mean(x)) ./ std(x)

df = CSV.read(joinpath(@__DIR__, "data", "milk.csv"))
categorical!(df, :clade)

df = df |> @mutate(clade_id = levelcode(_.clade)) |> DataFrame

df.K = scale(df.kcal_per_g)

K = df.K
clade = df.clade_id

@model model5_9(K, clade) = begin
    n_clades = length(unique(clade))
    σ ~ Exponential(1)
    α ~ filldist(Normal(0, 0.5), n_clades)
    μ = α[clade]
    @. K ~ Normal(μ, σ)
end

model = model5_9(K, clade)

m5_9 = sample(model, NUTS(), 2000)
loglik = logprob"K = K | chain = m5_9, clade = clade, model = model5_9"

describe(m5_9)
plot(m5_9)

params = ["α["] .* map(string, 1:4) .* "]"
dict = Dict(zip(params, levels(df.clade)))

post = DataFrame(m5_9)
rename!(post, dict)
post_long = DataFrames.stack(post)

using VegaLite

post_long |>
  @filter(_.variable != "σ") |>

  @vlplot(y={"variable:o", title = ""}) +
  @vlplot(
    mark={
        :point,
        filled=true
    },
    x={
        "mean(value)",
        scale={zero=false},
        title="expected kcal (std)"},
    color={value=:black}) +
  @vlplot(
    mark={
        :errorbar,
        extent=:stdev
     },
     x={:value, title="expected kcal (std)"})

using Random
Random.seed!(63);

df.house = sample(repeat(1:4, inner = 8), nrow(df))

@model model5_10(K, clade, house) = begin
    n_clades = length(unique(clade))
    n_houses = length(unique(house))
    σ ~ Exponential(1)
    α ~ filldist(Normal(0, 0.5), n_clades)
    h ~ filldist(Normal(0, 0.5), n_houses)
    μ = α[clade] + h[house]
    @. K ~ Normal(μ, σ)
end

K, clade, house = df.K, df.clade_id, df.house
model = model5_10(K, clade, house)

m5_10 = sample(model, NUTS(), 2000)

post = DataFrame(m5_10)

function recover_levels(var::String, factor::CategoricalArray)
    flevels = levels(factor)
    nlevels = length(flevels)
    params = [var * "["] .* map(string, 1:nlevels) .* "]"
    dict = Dict(zip(params, flevels))
    return dict
end

clades = recover_levels("α", df.clade)
rename!(post, clades)

houses = categorical(["Gryffindor", "Ravenclaw", "Hufflepuff", "Slytherin"])
houses = recover_levels("h", houses)
rename!(post, houses)

post_long = DataFrames.stack(post)

post_long |>
  @filter(_.variable != "σ") |>

  @vlplot(y={"variable:o", title = ""}) +
  @vlplot(
    mark={
        :point,
        filled=true
    },
    x={
        "mean(value)",
        scale={zero=false},
        title="expected kcal (std)"},
    color={value=:black}) +
  @vlplot(
    mark={
        :errorbar,
        extent=:stdev
     },
     x={:value, title="expected kcal (std)"})


text:
"""
Chu, C.S., Tsang, I.W. and Kwok, J.T., 2004, July.
Scaling up support vector data description by using core-sets.
In 2004 IEEE International Joint Conference on Neural Networks (IEEE Cat. No. 04CH37541) (Vol. 1, pp. 425-430). IEEE.

We dub the method core-set SVDD (CSSVDD)

Parameter guideline from paper:
num_init_sample = 20
k = 10
delta = 0.01*eps
eps none (experiments eps ∈ [0.15, 0.2, 0.3, 0.4, 0.5])
init_strategy: any strategy for gamma, C=1

Implementation notes:
1) We do not use fast incremental SVDD updates but fit a model from scratch each
   iteration.
2) We fit the initial SVDD not only on the approximal center and farthest point
   from the initial sample but on all points of the initial sample to avoid
   solver errors.
"""
mutable struct CSSVDD <: Sampler
    num_init_sample::Int
    k::Int
    delta::Real
    eps::Real
    outlier_percentage::Real
    init_strategy::SVDD.InitializationStrategy
    solver::OptimizerFactory
    function CSSVDD(num_init_sample::Int, k::Int, delta::Real, eps::Real, outlier_percentage::Real,
                    init_strategy::SVDD.InitializationStrategy, solver::OptimizerFactory)
        num_init_sample > 0 || throw(ArgumentError("Non-positive size of initial sample $(num_init_sample)."))
        k > 0 || throw(ArgumentError("Non-positive value k to initialize R $(k)."))
        delta > 0 || throw(ArgumentError("Non-positive value for delta $(delta)."))
        eps > 0 || throw(ArgumentError("Non-positive value for eps $(eps)."))
        outlier_percentage >= 0 || outlier_percentage <= 1 || throw(ArgumentError("Outlier percentage not in [0, 1] $(outlier_percentage)."))
        return new(num_init_sample, k, delta, eps, outlier_percentage, init_strategy, solver)
    end
end

function CSSVDD(eps::Real, outlier_percentage::Real, init_strategy::SVDD.InitializationStrategy, solver::OptimizerFactory,
                num_init_sample::Int=20, k::Int=10)
    return CSSVDD(num_init_sample, k, 0.01*eps, eps, outlier_percentage, init_strategy, solver)
end

function sample(sampler::CSSVDD, data::DataSet, labels::Vector{Symbol})::BitArray{1}
    n = size(data, 2)
    sampler.num_init_sample = min(n, sampler.num_init_sample)
    # Setup model
    model = SVDD.VanillaSVDD(data)
    SVDD.initialize!(model, sampler.init_strategy)
    SVDD.set_adjust_K!(model, true)
    K = copy(model.K_adjusted)
    # Step 1: Initialize R and c
    # Choose initial random sample and learn SVDD
    initial_sample = randperm(size(data, 2))[1:sampler.num_init_sample]
    SVDD.set_data!(model, data[:, initial_sample])
    try
        SVDD.fit!(model, sampler.solver)
    catch e
        throw(SamplingException("Failed fitting model on data set of size $(size(model.data)) due to $e."))
    end
    # Find closest point to center
    c = initial_sample[findmin(SVDD.predict(model, data[:, initial_sample]))[2]]
    # Set working set S (see implementation notes)
    S = copy(initial_sample)
    # Find furthest point from the random point (point 1) in the sample to init R
    # Distance is computed in the kernel space
    R_prev = minimum(K[:, initial_sample[1]]) / sampler.k
    while true
        # Step 2: Find set of points outside of (1+eps)-ball
        if length(S) == 1
            P = K[:, c] .> (1 + sampler.eps) * R_prev
        else
            model.R = (1 + sampler.eps) * model.R
            P = SVDD.predict(model, data) .> 0
            R_prev = model.R
        end
        # Step 3: check termination criteria
        if count(P) <= sampler.outlier_percentage * n
            break
        end
        # Step 4: Expand core-set with closest point to c in P that is not in S
        P_indices = [i for i in findall(P) if i ∉ S]
        if isempty(P_indices)
            break
        end
        p_close_out = P_indices[findmin(K[P_indices, c])[2]]
        push!(S, p_close_out)
        # Step 5: refit SVDD
        SVDD.set_data!(model, data[:, S])
        try
            SVDD.fit!(model, sampler.solver)
        catch e
            throw(SamplingException("Failed fitting model on data set of size $(size(model.data)) due to $e."))
        end
        # Step 6: enforce radius constraint
        if model.R < (1 + sampler.delta * sampler.eps) * R_prev
            model.R += sampler.delta * sampler.eps
        end
    end
    sample_mask = falses(n)
    sample_mask[S] .= true
    return sample_mask
end


text:
module DateExtensions

using Dates
# A Timestamp type couples an Instant type with a Calendar type
# to provide convenient human-conversion rules carried out
# by multiple dispatch.
immutable Timestamp{T<:Instant,C<:Calendar} <: Dates.TimeType
    instant::T
    Timestamp(x::T) = new(x)
end

# Holidays
newyears(y) = (y,1,1)
isnewyears(dt::Dates.TimeType) = yearmonthday(dt) == newyears(year(dt))
independenceday(y) = (y,7,4)
isindependenceday(dt::Dates.TimeType) = yearmonthday(dt) == independenceday(year(dt))
veteransday(y) = (y,11,11)
isveteransday(dt::Dates.TimeType) = yearmonthday(dt) == veteransday(year(dt))
christmas(y) = (y,12,25)
ischristmas(dt::Dates.TimeType) = yearmonthday(dt) == christmas(year(dt))

function ismartinlutherking(dt::Dates.TimeType)
    # 3rd Monday of January
    dayofweek(dt) == Mon && month(dt) == Jan && dayofweekofmonth(dt) == 3
end
function ispresidentsday(dt::Dates.TimeType)
    # 3rd Monday of February
   dayofweek(dt) == Mon && month(dt) == Feb && dayofweekofmonth(dt) == 3
end
function ismemorialday(dt::Dates.TimeType)
    # Last Monday of May
   dayofweek(dt) == Mon && month(dt) == May && dayofweekofmonth(dt) == daysofweekinmonth(dt)
end
function islaborday(dt::Dates.TimeType)
    # First Monday of September
    dayofweek(dt) == Mon && month(dt) == Sep && dayofweekofmonth(dt) == 1
end
function iscolumbusday(dt::Dates.TimeType)
    # Second Monday of October
    dayofweek(dt) == Mon && month(dt) == Oct && dayofweekofmonth(dt) == 2
end
function isthanksgiving(dt::Dates.TimeType)
    # Fourth Thursday in November
    dayofweek(dt) == Thu && month(dt) == Nov && dayofweekofmonth(dt) == 4
end

function easter(y)
    # Butcher's Algorithm: http://www.smart.net/~mmontes/butcher.html
    a=y%19
    b=div(y,100)
    c=y%100
    d=div(b,4)
    e=b%4
    f=div(b+8,25)
    g=div(b-f+1,3)
    h=(19*a+b-d-g+15)%30
    i=div(c,4)
    k=c%4
    l=(32+2*e+2*i-h-k)%7
    m=div(a+11*h+22*l,451)
    month=div(h+l-7*m+114,31)
    p=(h+l-7*m+114)%31
    return (y,month,p+1)
end
iseaster(dt::Dates.TimeType) = yearmonthday(dt) == easter(year(dt))

const HOLIDAYS = x->isnewyears(x) || isindependenceday(x) ||
                    isveteransday(x) || ischristmas(x) ||
                    ismartinlutherking(x) || ispresidentsday(x) ||
                    ismemorialday(x) || islaborday(x) ||
                    iscolumbusday(x) || isthanksgiving(x)
const OBSERVEDHOLIDAYS = x->begin
    # If the holiday is on a weekday
    if HOLIDAYS(x) && dayofweek(x) < Saturday
        return true
    # Holiday is observed Monday if falls on Sunday
    elseif dayofweek(x) == 1 && HOLIDAYS(x-Day(1))
        return true
    # Holiday is observed Friday if falls on Saturday
    elseif dayofweek(x) == 5 && HOLIDAYS(x+Day(1))
        return true
    else
        return false
    end
end

end # module


text:
# ---
# title: 1679. Max Number of K-Sum Pairs
# id: problem1679
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Hash Table
# link: <https://leetcode.com/problems/max-number-of-k-sum-pairs/description/>
# hidden: true
# ---
#
# You are given an integer array `nums` and an integer `k`.
#
# In one operation, you can pick two numbers from the array whose sum equals `k`
# and remove them from the array.
#
# Return _the maximum number of operations you can perform on the array_.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,2,3,4], k = 5
#     Output: 2
#     Explanation: Starting with nums = [1,2,3,4]:
#     - Remove numbers 1 and 4, then nums = [2,3]
#     - Remove numbers 2 and 3, then nums = []
#     There are no more pairs that sum up to 5, hence a total of 2 operations.
#
# **Example 2:**
#
#
#
#     Input: nums = [3,1,3,4,3], k = 6
#     Output: 1
#     Explanation: Starting with nums = [3,1,3,4,3]:
#     - Remove the first two 3's, then nums = [1,4,3]
#     There are no more pairs that sum up to 6, hence a total of 1 operation.
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 105`
#   * `1 <= nums[i] <= 109`
#   * `1 <= k <= 109`
#
#
## @lc code=start
using LeetCode

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


text:
using Reachability, MathematicalSystems, MathematicalPredicates, HybridSystems,
      SparseArrays

function print_dynamics(A, b, location_name)
    println("dynamics of location $location_name:")
    for i in 1:size(A, 1)-1  # ignore the last dimension (time)
        print("x_$i' = ")
        for j in 1:size(A, 2)
            if !iszero(A[i,j])
                print("$(A[i,j]) x_$j + ")
            end
        end
        println("$(b[i])\n")
    end
end

"""
    powertrain(θ::Int)

Return the hybrid system that models a mechanical system with backlash
(powertrain) from an automotive problem.

### Input

- `θ`        -- (optional, default: `1`) number of rotating masses
- `X0_scale` -- (optional, default: `1.0`) scaling factor in ``(0, 1]`` for
                modifying the initial states; for value `1.0` we do not modify
                the initial states and let `X0` be just a line segment;
                otherwise we let `X0' = δ ⋅ X0 + ⊕ {(1-δ) · center(X0)}` where
                `δ` is the factor

### Output

Hybrid system representing the powertrain model.

### Notes

The model is based on [1], which is an extension with a parametric number of
rotating masses of the automotive powertrain model presented in [2].
The parameters of the system's additional masses, which can be interpreted as
rotating elements in a gearbox and further powertrain elements, are taken from
[3].

[1] M. Althoff and B. H. Krogh, Avoiding Geometric Intersection Operations in
Reachability Analysis of Hybrid Systems [HSCC '12 Proceedings of the 15th ACM
international conference on Hybrid Systems: Computation and Control Pages
45-54](https://dl.acm.org/citation.cfm?id=2185643).

[2] A. Lagerberg. A benchmark on hybrid control of an automotive powertrain with
backlash. Technical Report R005/2007, Signals and Systems, Chalmers University
of Technology, 2007.

[3] E.-A. M. A. Rabeih. Torsional Vibration Analysis of Automotive Drivelines.
PhD thesis, University of Leeds, 1997.
"""
function powertrain(θ::Int=1; X0_scale::Float64=1.0)
    @assert θ > 0 "θ must be positive, but was $θ"
    @assert (X0_scale > 0.0 && X0_scale <= 1.0) "scale $X0_scale ∉ (0, 1]"

    # activate for printing the dynamics
    display_dynamics = false

    # dimension of state space (last dimension is time)
    n = 2 * θ + 7 + 1

    # constants
    # indices 'm' (ₘ) and 'l' (ₗ) refer to the motor and the load
    # indices 'i' (ᵢ) refer to the numbering of additional rotating masses
    t_init = 0.2  # time to stay in the initial location
    α = 0.03  # backlash size (half of the gap width)
    τ_eng = 0.1  # engine time constant
    γ = 12.0 # gearbox ratio (dimensionless)
    u = 5.0 # requested engine torque
    # moments of inertia [kg m²]
    Jₗ = 140.
    Jₘ = 0.3
    Jᵢ = 100.0  # TODO the paper says 0.01, the SpaceEx model uses 100
    # viscous friction constants [Nm s/rad]
    bₗ = 5.6
    bₘ = 0.0
    bᵢ = 1.0
    # shaft stiffness [Nm/rad]
    kᵢ = 1e5
    kₛ = 1e4
    # PID parameters
    k_P = 0.5  # [Nms/rad]
    k_I = 0.5  # [Nm/rad]
    k_D = 0.5  # [Nms²/rad]

    # additional bloating to get a non-flat guard intersection
    guard_bloating = sqrt(eps(Float64))

    function get_dynamics(kₛ, α, u)
        # physical units of the state variables
        # [x₁] = rad
        # [x₂] = Nm
        # [x₃] = rad
        # [x₄] = rad/s
        # [x₅] = rad
        # [x₆] = rad/s
        # [x₇] = rad/s
        # [x₈] = rad
        # [x₉] = rad/s
        # ...
        # [x_(2θ+6)] = rad
        # [x_(2θ+7)] = rad/s

        # linear dynamics
        A = spzeros(n, n)

        A[1, 7] = 1.0 / γ
        A[1, 9] = -1.0

        A[2, 1] = (-k_I * γ + k_D * kₛ / (γ * Jₘ)) / τ_eng
        A[2, 2] = (-k_D / Jₘ - 1.0) / τ_eng
        A[2, 3] = k_I * γ / τ_eng
        A[2, 4] = k_P * γ / τ_eng
        A[2, 7] = (-k_P + k_D * bₘ /Jₘ) / τ_eng
        A[2, 8] = -k_I * γ / τ_eng

        A[3, 4] = 1.0

        A[5, 6] = 1.0

        A[6, 5] = -kᵢ / Jₗ
        A[6, 6] = -bₗ / Jₗ
        A[6, 2*θ+6] = kᵢ / Jₗ

        A[7, 1] = -kₛ / (Jₘ * γ)
        A[7, 2] = 1.0 / Jₘ
        A[7, 7] = -bₘ / Jₘ

        i = 8
        while i < n-1
            A[i, i+1] = 1.0

            if i == 8
                # x9 has special dynamics
                A[i+1, 1] = kₛ / Jᵢ
                A[i+1, i] = -kᵢ / Jᵢ
            else
                A[i+1, i-2] = kᵢ / Jᵢ
                A[i+1, i] = -2. * kᵢ / Jᵢ
            end
            A[i+1, i+1] = -bᵢ / Jᵢ

            # wrap-around to x5 in the last step
            j = (i == n-2) ? 5 : i+2
            A[i+1, j] = kᵢ / Jᵢ

            i += 2
        end

        # affine vector
        b = spzeros(n)
        b[2] = k_D * (γ * u - kₛ * α / (Jₘ * γ)) / τ_eng
        b[4] = u
        b[7] = kₛ * α / (Jₘ * γ)
        b[9] = -kₛ * α / Jᵢ
        b[n] = 1.0  # time

        return A, b
    end

    # hybrid automaton
    automaton = LightAutomaton(4)

    # negAngle
    A, b = get_dynamics(kₛ, -α, u)
    X = HalfSpace(sparsevec([1], [1.], n), -α)  # x1 <= -α
    m_negAngle = CACS(A, b, X)
    if display_dynamics
        print_dynamics(A, b, "negAngle")
    end

    # deadzone
    A, b = get_dynamics(0., -α, u)
    X = HPolyhedron([HalfSpace(sparsevec([1], [-1.], n), α),  # x1 >= -α
                     HalfSpace(sparsevec([1], [1.], n), α)])  # x1 <= α
    m_deadzone = CACS(A, b, X)
    if display_dynamics
        print_dynamics(A, b, "deadzone")
    end

    # posAngle
    A, b = get_dynamics(kₛ, α, u)
    X = HalfSpace(sparsevec([1], [-1.], n), -α + guard_bloating)  # x1 >= α - ε
    m_posAngle = CACS(A, b, X)
    if display_dynamics
        print_dynamics(A, b, "posAngle")
    end

    # negAngleInit
    A, b = get_dynamics(kₛ, -α, -u)
    X = HalfSpace(sparsevec([n], [1.], n), t_init)  # t <= t_init
    m_negAngleInit = CACS(A, b, X)
    if display_dynamics
        print_dynamics(A, b, "negAngleInit")
    end

    # modes
    modes = [m_negAngle, m_deadzone, m_posAngle, m_negAngleInit]

    # transition negAngleInit -> negAngle
    add_transition!(automaton, 4, 1, 1)
    guard = HalfSpace(sparsevec([n], [-1.], n), -t_init + guard_bloating)  # t >= t_init - ε
    r_41 = ConstrainedIdentityMap(n, guard)

    # transition negAngle -> deadzone
    add_transition!(automaton, 1, 2, 2)
    guard = HalfSpace(sparsevec([1], [-1.], n), α)  # x1 >= -α
    r_12 = ConstrainedIdentityMap(n, guard)

    # transition deadzone -> posAngle
    add_transition!(automaton, 2, 3, 3)
    guard = HalfSpace(sparsevec([1], [-1.], n), -α)  # x1 >= α
    r_23 = ConstrainedIdentityMap(n, guard)

    # TODO the SpaceEx model does not contain the following transitions
#     # transition deadzone -> negAngle
#     add_transition!(automaton, 2, 1, 4)
#     guard = HalfSpace(sparsevec([1], [1.], n), -α)  # x1 <= -α
#     r_21 = ConstrainedIdentityMap(n, guard)
#     # transition posAngle -> deadzone
#     add_transition!(automaton, 3, 2, 5)
#     guard = HalfSpace(sparsevec([1], [1.], n), α)  # x1 <= α
#     r_32 = ConstrainedIdentityMap(n, guard)

    # transition annotations
    resetmaps = [r_41, r_12, r_23]

    # switching
    switchings = [HybridSystems.AutonomousSwitching()]

    ℋ = HybridSystem(automaton, modes, resetmaps, switchings)

    # initial condition in mode 1
    c = Vector{Float64}(undef, n)
    g = Vector{Float64}(undef, n)
    c[1:7] = [-0.0432, -11., 0., 30., 0., 30., 360.]
    g[1:7] = [0.0056, 4.67, 0., 10., 0., 10., 120.]
    i = 8
    while i < n
        c[i] = -0.0013
        g[i] = 0.0006
        i += 1
        c[i] = 30.
        g[i] = 10.
        i += 1
    end
    c[n] = 0.0
    g[n] = 0.0
    if X0_scale < 1.0
        g = X0_scale * g
    end
    X0 = Zonotope(c, hcat(g))
    initial_condition = [(4, X0)]

    system = InitialValueProblem(ℋ, initial_condition)

    # safety property
    property_2 = is_disjoint_from(HalfSpace(sparsevec([1], [1.], n), -α))  # x1 <= -α
    property_3 = is_disjoint_from(HalfSpace(sparsevec([1], [1.], n), α))   # x1 <= α
    property = Dict(2 => property_2, 3 => property_3)

    # default options
    options = Options(:T=>2.0, :property=>property)

    return (system, options)
end

function run_powertrain(system, options)
    opC = BFFPSV18(:δ => 0.0005, :assume_sparse => true)
    opD = LazyDiscretePost(:lazy_R⋂I => true, :lazy_R⋂G => true)
    options[:mode] = "check"
    options[:plot_vars] = [1, 3]

    solve(system, options, opC, opD)
end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
## IMPORTS ##
using Pkg
Pkg.activate(joinpath(@__DIR__,"../.."))
using Dates
include(joinpath(@__DIR__,"model.jl"))


text:
using HDF5
using Test

@testset "external" begin

# roughly following https://www.hdfgroup.org/ftp/HDF5/current/src/unpacked/examples/h5_extlink.c
fn1 = tempname()
fn2 = tempname()

source_file = h5open(fn1, "w")
agroup = g_create(source_file, "agroup")
target_file = h5open(fn2, "w")
target_group = g_create(target_file, "target_group")
target_group["abc"] = "abc"
target_group["1"] = 1
target_group["1.1"] = 1.1
close(target_file)

# create external link such that source_file["ext_link"] points to target_file["target_group"]
# test both an HDF5.File and an HDF5.Group for first argument
HDF5.create_external(source_file, "ext_link", target_file.filename, "target_group")
HDF5.create_external(agroup, "ext_link", target_file.filename, "target_group")
# write some things via the external link
new_group = g_create(source_file["ext_link"], "new_group")
new_group["abc"] = "abc"
new_group["1"] = 1
new_group["1.1"] = 1.1

# read things from target_group via exernal link created with HDF5File argument
group = source_file["ext_link"]
@test read(group["abc"]) == "abc"
@test read(group["1"]) == 1
@test read(group["1.1"]) == 1.1
# read things from target_group via the external link created with HDF5.Group argument
groupalt = source_file["agroup/ext_link"]
@test read(groupalt["abc"]) == "abc"
@test read(groupalt["1"]) == 1
@test read(groupalt["1.1"]) == 1.1
close(source_file)

##### tests that should be included but don't work
# when ggggggggg restarts julia and keeps track of target_file.filename,
# these tests succeed
# reopening the target_file crashes due to "file close degree doesn't match"
# target_file = h5open(target_file.filename, "r")
# group2 = target_file["target_group"]["new_group"]
# @test read(group2["abc"])=="abc"
# @test read(group2["1"])==1
# @test read(group2["1.1"])==1.1

rm(fn1)
# rm(fn2)

end # testset external


text:
@require CMAEvolutionStrategy = "8d3b24bd-414e-49e0-94fb-163cc3a3e411" begin
    mutable struct CMA
        x0
        sigma0
        lower
        upper
        kwargs
    end
    cma(prior; kwargs...) = cma(lower(prior), upper(prior); kwargs...)
    """
        cma(lower, upper;
            sigma0 = 0.25,
            kwargs...)
        cma(prior; kwargs...) = pycma(lower(prior), upper(prior); kwargs...)

    CMA optimizer for [`PointEstimator`](@ref), available when
    `using CMAEvolutionStrategy`.
    `kwargs` are passed to `CMAEvolutionStrategy.minimize`.
    """
    function cma(lower, upper;
                 sigma0 = 0.25,
                 kwargs...)
        d = length(lower)
        if d == 1
            error("CMA-ES is not recommended to be used in 1 dimension.")
        end
        CMA(rand(d) .* (upper .- lower) + lower, sigma0,
            lower, upper, kwargs)
    end
    export cma
    function _optimize(o::CMA, f; maxfevals, verbose = 1)
        res = CMAEvolutionStrategy.minimize(f, o.x0, o.sigma0;
                                            lower = o.lower, upper = o.upper,
                                            o.kwargs...)
        (x = xbest(res), f = fbest(res), res = res)
    end
end
@require BlackBoxOptim = "a134a8b2-14d6-55f6-9291-3336d3ab0209" begin
    mutable struct BBO
        options
    end
    bbo(prior; kwargs...) = bbo(lower(prior), upper(prior); kwargs...)
    """
        bbo(lower, upper; kwargs...)
        bbo(prior; kwargs...) = bbo(lower(prior), upper(prior); kwargs...)

    BlackBoxOptim optimizer for [`PointEstimator`](@ref), available
    when `using BlackBoxOptim`. `kwargs` are passed to `BlackBoxOptim.bboptimize`.
    """
    function bbo(lower, upper; kwargs...)
        srange = [(l, u) for (l, u) in zip(lower, upper)]
        BBO((SearchRange = srange, pairs(kwargs)...))
    end
    export bbo
    function _optimize(o::BBO, f; maxfevals, verbose)
        verbose = typeof(verbose) <: Bool ?
                    (verbose ? :verbose : :silent) : verbose
        res = BlackBoxOptim.bboptimize(f; MaxFuncEvals = maxfevals,
                                       TraceMode = verbose,
                                       o.options...)
        (x = BlackBoxOptim.best_candidate(res),
         f = BlackBoxOptim.best_fitness(res),
         res = res)
    end
end
@require BayesianOptimization = "4c6ed407-134f-591c-93fa-e0f7c164a0ec" begin
    const BO = BayesianOptimization
    mutable struct BOTmp
        model
        modeloptimizer
        acquisition
        lb
        ub
        options
        opt
    end
    bo(prior; kwargs...) = bo(lower(prior), upper(prior); kwargs...)
    """
        bo(lower, upper;
                acquisition = BO.UpperConfidenceBound(),
                capacity = 10^3,
                mean = BO.GaussianProcesses.MeanConst(0.),
                kernel = BO.GaussianProcesses.SEArd(zeros(length(lower)), 5.),
                logNoise = 0,
                modeloptimizer_options = NamedTuple(),
                kwargs...)
        bo(prior; kwargs...) = bo(lower(prior), upper(prior); kwargs...)

    BayesianOptimization optimizer for [`PointEstimator`](@ref), available
    when `using BayesianOptimization`. `kwargs` are passed to `BayesianOptimization.BOpt`.
    """
    function bo(lower, upper;
                acquisition = BO.UpperConfidenceBound(),
                capacity = 10^3,
                mean = BO.GaussianProcesses.MeanConst(0.),
                kernel = BO.GaussianProcesses.SEArd(zeros(length(lower)), 5.),
                logNoise = 0,
                modeloptimizer_options = NamedTuple(),
                kwargs...)
        model = BO.GaussianProcesses.ElasticGPE(length(lower),
                           mean = mean,
                           kernel = kernel,
                           logNoise = 0.,
                           capacity = capacity)
        modeloptimizer_options = merge((every = 125,
                                        noisebounds = [-4, 3],
                                        kernbounds = [[log.((upper .- lower)./10^2); 0],
                                                      [log.(2*(upper .- lower)); 10]],
                                        maxeval = 40),
                                       modeloptimizer_options)
        modeloptimizer = BO.MAPGPOptimizer(; modeloptimizer_options...)
        BOTmp(model, modeloptimizer, acquisition, lower, upper,
               merge((repetitions = 1, maxiterations = 500,
                      acquisitionoptions = (restarts = 20, maxtime = 2),
                      sense = BO.Min), kwargs), nothing)
    end
    export bo
    function _optimize(o::BOTmp, f; maxfevals, verbose)
        verbose = verbose ? BO.Progress : BO.Silent
        opt = BO.BOpt(f, o.model, o.acquisition, o.modeloptimizer, o.lb, o.ub;
                   o.options..., verbosity = verbose,
                   maxiterations = maxfevals)
        o.opt = opt
        res = BO.boptimize!(opt)
        (x = res.model_optimizer, f = res.model_optimum, res = opt)
    end
end
@require SimultaneousPerturbationStochasticApproximation = "1663b253-47ef-444c-a4fb-e919f25dc38d" begin
    struct SPSATmp
        s
    end
    spsa(prior; kwargs...) = spsa(lower(prior), upper(prior); kwargs...)
    """
        spsa(lower, upper; kwargs...)
        spsa(prior; kwargs...) = spsa(lower(prior), upper(prior); kwargs...)

    SimultaneousPerturbationStochasticApproximation optimizer for [`PointEstimator`](@ref), available when `using SimultaneousPerturbationStochasticApproximation`. `kwargs` are passed to `SimultaneousPerturbationStochasticApproximation.SPSA`.
    """
    function spsa(lower, upper; kwargs...)
        SPSATmp(SimultaneousPerturbationStochasticApproximation.SPSA(; lower = lower, upper = upper, kwargs...))
    end
    export spsa
    function _optimize(o::SPSATmp, f; maxfevals, verbose)
        res = SimultaneousPerturbationStochasticApproximation.minimize!(o.s, f, maxfevals = maxfevals, verbose = verbose)
        (x = res[1].x, f = res[1].f, res = res)
    end
end
@require CMAEvolutionStrategy = "8d3b24bd-414e-49e0-94fb-163cc3a3e411" begin
    cma(prior; kwargs...) = cma(lower(prior), upper(prior); kwargs...)
    function cma(lower, upper;
                 K = Ref(0), minK = 2, maxK = 10^3, alphaK = 1.5,
                 maxfevals = typemax(Int),
                 kwargs...)
        n = length(lower)
        x0 = (upper .- lower) .* rand(n) .+ lower
        cb = s -> begin
            if s > 0
                if K[] < maxK
                    K[] = min(maxK, round(Int, K[] * alphaK))
                    return true
                else
                    return true
                end
            else
                if K[] > minK
                    K[] = max(minK, round(Int, K[] / alphaK^.25))
                end
                return false
            end
        end
        noise = CMAEvolutionStrategy.NoiseHandling(n, callback = cb)
        CMAEvolutionStrategy.Optimizer(x0, .25; noise_handling = noise,
                                       lower = lower, upper = upper,
                                       maxfevals = maxfevals, kwargs...)
    end
    export cma
    function _optimize(o::CMAEvolutionStrategy.Optimizer, f; maxfevals, verbose)
        o.stop.maxfevals = maxfevals
        o.logger.verbosity = verbose
        res = CMAEvolutionStrategy.run!(o, f)
        (x = CMAEvolutionStrategy.population_mean(res),
         f = CMAEvolutionStrategy.fbest(res), res = res)
    end
end


text:
mutable struct SubOcean
    master_in_flds    :: InputFields
    master_ocn        :: Ocean
    worker_ocn        :: Ocean
    block_id          :: Integer

    pull_fr_rng      :: UnitRange
    push_fr_rng      :: UnitRange
    push_to_rng      :: UnitRange

    pull_fr_rng_bnd  :: AbstractArray{Union{UnitRange, Nothing}}
    pull_to_rng_bnd  :: AbstractArray{Union{UnitRange, Nothing}}
    push_fr_rng_bnd  :: AbstractArray{Union{UnitRange, Nothing}}
    push_to_rng_bnd  :: AbstractArray{Union{UnitRange, Nothing}}

end

function calParallizationRange(;
    N = Integer,     # Total grids
    P = Integer,     # Number of procs
    L = Integer,     # Overlapping grids
)

    if ! (N >= max(1, L) * P)
        throw(ErrorException("Condition must be satisfied: N >= max(1, L) * P"))
    end

    n̄ = floor(Integer, N / P)
    R = N - n̄ * P

    pull_fr_rngs = Array{Union{UnitRange, Nothing}}(undef, P)
    push_fr_rngs = Array{Union{UnitRange, Nothing}}(undef, P)
    push_to_rngs = Array{Union{UnitRange, Nothing}}(undef, P)

    # 1: lower latitude side (south), 2: higher latitude side (north)
    pull_fr_rngs_bnd = Array{Union{UnitRange, Nothing}}(undef, P, 2)
    pull_to_rngs_bnd = Array{Union{UnitRange, Nothing}}(undef, P, 2)
    push_fr_rngs_bnd = Array{Union{UnitRange, Nothing}}(undef, P, 2)
    push_to_rngs_bnd = Array{Union{UnitRange, Nothing}}(undef, P, 2)

    cnt = 1
    for p = 1:P
        m = (p <= R) ? n̄ + 1 : n̄  # assigned grids

        pull_fr_rngs[p] = cnt-L:cnt+m-1+L
        push_fr_rngs[p] = L+1:L+m
        push_to_rngs[p] = cnt:cnt+m-1

        # Boundary
        pull_fr_rngs_bnd[p, 1] = cnt-L:cnt-1
        pull_fr_rngs_bnd[p, 2] = cnt+m:cnt+m+L-1

        pull_to_rngs_bnd[p, 1] = 1:L
        pull_to_rngs_bnd[p, 2] = L+m+1:L+m+L

        push_fr_rngs_bnd[p, 1] = L+1:L+L
        push_fr_rngs_bnd[p, 2] = L+m-L+1:L+m

        push_to_rngs_bnd[p, 1] = cnt:cnt+L-1
        push_to_rngs_bnd[p, 2] = cnt+m-L:cnt+m-1

        cnt += m
    end

    # South pole and north pole do not have boundaries
    pull_fr_rngs_bnd[1, 1] = nothing
    pull_to_rngs_bnd[1, 1] = nothing
    push_fr_rngs_bnd[1, 1] = nothing
    push_to_rngs_bnd[1, 1] = nothing

    pull_fr_rngs_bnd[end, 2] = nothing
    pull_to_rngs_bnd[end, 2] = nothing
    push_fr_rngs_bnd[end, 2] = nothing
    push_to_rngs_bnd[end, 2] = nothing

    # Adjust the first and last range (south pole and north pole)
    pull_fr_rngs[1] = (pull_fr_rngs[1][1]+L):pull_fr_rngs[1][end]
    pull_fr_rngs[end] = pull_fr_rngs[end][1]:(pull_fr_rngs[end][end]-L)
    push_fr_rngs[1] = 1:length(push_fr_rngs[1])

    # Change range because to southmost boundary is trimmed
    if P > 1
        pull_to_rngs_bnd[1, 2] = pull_to_rngs_bnd[1, 2] .- L
        push_fr_rngs_bnd[1, 2] = push_fr_rngs_bnd[1, 2] .- L
    end

    return pull_fr_rngs, push_fr_rngs, push_to_rngs, pull_fr_rngs_bnd, pull_to_rngs_bnd, push_fr_rngs_bnd, push_to_rngs_bnd

end

function makeSubOcean(;
    master_ocn   :: Ocean,
    block_id     :: Integer,
    pull_fr_rng  :: UnitRange,
    push_fr_rng  :: UnitRange,
    push_to_rng  :: UnitRange,
    pull_fr_rng_bnd :: AbstractArray{Union{UnitRange, Nothing}},
    pull_to_rng_bnd :: AbstractArray{Union{UnitRange, Nothing}},
    push_fr_rng_bnd :: AbstractArray{Union{UnitRange, Nothing}},
    push_to_rng_bnd :: AbstractArray{Union{UnitRange, Nothing}},

#    nblocks      :: Integer,
)
#=
    overlap_grids = 2

    println(format("{:03d} Entering makeSubOcean.", block_id))

    touch_southpole = block_id == 1
    touch_northpole = block_id == nblocks

    sub_Ny_wo_ghost = ceil(Integer, master_ocn.Ny / nblocks)

    if block_id != nblocks
        sub_Ny = sub_Ny_wo_ghost
    else
        sub_Ny = master_ocn.Ny - (nblocks-1) * sub_Ny_wo_ghost
    end

    if sub_Ny <= 0
        throw(ErrorException("sub_Ny <= 0. Please check your resolution and nblocks"))
    end

    pull_fr_beg_y = push_to_beg_y = sub_Ny_wo_ghost * (block_id - 1) + 1
    pull_fr_end_y = push_to_end_y = pull_fr_beg_y + sub_Ny - 1

    push_fr_beg_y = 1
    push_fr_end_y = sub_Ny

    if ! touch_southpole
        # expand south boundary
        pull_fr_beg_y -= 1
        sub_Ny += 1

        # fix push from range.
        # We want to skip the expanded latitude
        push_fr_beg_y += 1
        push_fr_end_y += 1
    end

    if ! touch_northpole
        # expand north boundary
        pull_fr_end_y += 1
        sub_Ny += 1

        # No need to fix push from range.
        # It is not affected.
    end

=#

    sub_Ny = length(pull_fr_rng)

    println("pull_fr_rng: ", pull_fr_rng)
    println("push_to_rng: ", push_to_rng)
    println("push_fr_rng: ", push_fr_rng)
    println("pull_fr_rng_bnd: ", pull_fr_rng_bnd)
    println("pull_to_rng_bnd: ", pull_to_rng_bnd)
    println("push_fr_rng_bnd: ", push_fr_rng_bnd)
    println("push_to_rng_bnd: ", push_to_rng_bnd)

    #=
    println("### rng3: ")
    println("pull_fr_rng3: ", pull_fr_rng3)
    println("push_to_rng3: ", push_to_rng3)
    println("push_fr_rng3: ", push_fr_rng3)
    =#

#=
    if length(pull_fr_rng2[2]) != sub_Ny
        throw(ErrorException("Pull from dimension does not match sub_Ny"))
    end

    if length(push_fr_rng2[2]) != length(push_to_rng2[2])
        throw(ErrorException("Push from and push to dimensions do not match."))
    end
=#

    return SubOcean(
        SubInputFields(master_ocn.in_flds, :, pull_fr_rng),
        master_ocn,
        Ocean(
            id             = block_id,
            gridinfo_file  = master_ocn.gi_file,
            sub_yrng       = pull_fr_rng,
            Nx             = master_ocn.Nx,
            Ny             = sub_Ny,
            zs_bone        = master_ocn.zs_bone,
            Ts             = master_ocn.Ts[:, :, pull_fr_rng],
            Ss             = master_ocn.Ss[:, :, pull_fr_rng],
            K_v            = master_ocn.K_v,
            Dh_T           = master_ocn.Dh_T,
            Dv_T           = master_ocn.Dv_T,
            Dh_S           = master_ocn.Dh_S,
            Dv_S           = master_ocn.Dv_S,
            T_ML           = master_ocn.T_ML[:, pull_fr_rng],
            S_ML           = master_ocn.S_ML[:, pull_fr_rng],
            h_ML           = master_ocn.h_ML[:, pull_fr_rng],
            h_ML_min       = master_ocn.h_ML_min[:, pull_fr_rng],
            h_ML_max       = master_ocn.h_ML_max[:, pull_fr_rng],
            we_max         = master_ocn.we_max,
            R              = master_ocn.R,
            ζ              = master_ocn.ζ,
            Ts_clim_relax_time = master_ocn.Ts_clim_relax_time,
            Ss_clim_relax_time = master_ocn.Ss_clim_relax_time,
            Ts_clim        = ( master_ocn.Ts_clim != nothing ) ? master_ocn.Ts_clim[:, :, pull_fr_rng] : nothing,
            Ss_clim        = ( master_ocn.Ss_clim != nothing ) ? master_ocn.Ss_clim[:, :, pull_fr_rng] : nothing,
            topo           = master_ocn.topo[:, pull_fr_rng],
            fs             = master_ocn.fs[:, pull_fr_rng],
            ϵs             = master_ocn.ϵs[:, pull_fr_rng],
            in_flds        = InputFields(:local, master_ocn.Nx, sub_Ny),
            arrange        = :zxy,
        ),
        block_id,
        pull_fr_rng,
        push_fr_rng,
        push_to_rng,
        pull_fr_rng_bnd,
        pull_to_rng_bnd,
        push_fr_rng_bnd,
        push_to_rng_bnd,
    )

end

function syncToMaster!(subocn::SubOcean;
        vars2 :: Any,
        vars3 :: Any,
)

    (subocn.worker_ocn.id == 0) && throw(ErrorException("`id` should not be 0 (master)."))

    master_ocn = subocn.master_ocn
    worker_ocn = subocn.worker_ocn

    w_rng = subocn.push_fr_rng
    m_rng = subocn.push_to_rng

    for var in vars2
        getfield(master_ocn, var)[:, m_rng] = view(getfield(worker_ocn, var), :, w_rng)
    end

    for var in vars3
        getfield(master_ocn, var)[:, :, m_rng] = view(getfield(worker_ocn, var), :, :, w_rng)
    end

end

function syncFromMaster!(
        subocn::SubOcean;
        vars2 :: Any,
        vars3 :: Any,
)

    (subocn.worker_ocn.id == 0) && throw(ErrorException("`id` should not be 0 (master)."))

    master_ocn = subocn.master_ocn
    worker_ocn = subocn.worker_ocn

    m_rng = subocn.pull_fr_rng
    w_rng = Colon()

    for var in vars2
        getfield(worker_ocn, var)[:, w_rng] = view(getfield(master_ocn, var), :, m_rng)
    end

    for var in vars3
        getfield(worker_ocn, var)[:, :, w_rng] = view(getfield(master_ocn, var), :, :, m_rng)
    end

end

function syncForcingFromMaster!(
    subocn::SubOcean
)

    (subocn.worker_ocn.id == 0) && throw(ErrorException("`id` should not be 0 (master)."))

    master_ocn = subocn.master_ocn
    worker_ocn = subocn.worker_ocn

    copyfrom!(worker_ocn.in_flds, subocn.master_in_flds)
end

function syncBoundaryToMaster!(subocn::SubOcean;
        vars2 :: Tuple = (),
        vars3 :: Tuple = (),
)

    (subocn.worker_ocn.id == 0) && throw(ErrorException("`id` should not be 0 (master)."))

    master_ocn = subocn.master_ocn
    worker_ocn = subocn.worker_ocn

    w_rng = subocn.push_fr_rng_bnd
    m_rng = subocn.push_to_rng_bnd

    for r=1:2
        if m_rng[r] != nothing

            for var in vars2
                getfield(master_ocn, var)[:, m_rng[r]]    = view(getfield(worker_ocn, var), :, w_rng[r])
            end

            for var in vars3
                getfield(master_ocn, var)[:, :, m_rng[r]] = view(getfield(worker_ocn, var), :, :, w_rng[r])
            end

        end
    end
end

function syncBoundaryFromMaster!(subocn::SubOcean;
        vars2 :: Tuple = (),
        vars3 :: Tuple = (),
)

    (subocn.worker_ocn.id == 0) && throw(ErrorException("`id` should not be 0 (master)."))

    master_ocn = subocn.master_ocn
    worker_ocn = subocn.worker_ocn

    m_rng = subocn.pull_fr_rng_bnd
    w_rng = subocn.pull_to_rng_bnd

    for r=1:2
        if m_rng[r] != nothing

            for var in vars2
                getfield(worker_ocn, var)[:, w_rng[r]]    = view(getfield(master_ocn, var), :, m_rng[r])
            end

            for var in vars3
                getfield(worker_ocn, var)[:, :, w_rng[r]] = view(getfield(master_ocn, var), :, :, m_rng[r])
            end

        end
    end

end

function init(ocn::Ocean)

    global  wkrs  =  workers()
    nwkrs = length(wkrs)

    println("Number of all workers: ", length(wkrs))

    (ocn.id == 0) || throw(ErrorException("`id` is not 0 (master). Id received: " * string(ocn.id)))

    # Assign `nothing` to prevent from sending large array of `View` objects.
    tmp_cols = ocn.cols
    tmp_lays = ocn.lays

    ocn.cols = nothing
    ocn.lays = nothing

    pull_fr_rngs, push_fr_rngs, push_to_rngs, pull_fr_rngs_bnd, pull_to_rngs_bnd, push_fr_rngs_bnd, push_to_rngs_bnd = calParallizationRange(N=ocn.Ny, P=nwkrs, L=2)

    @sync for (i, p) in enumerate(wkrs)
            # We have P processors, N workers, N blocks
            # Block ids are numbered from 1 to N
            @spawnat p let
                global subocn = makeSubOcean(
                    master_ocn  = ocn,
                    block_id    = i,
                    pull_fr_rng      = pull_fr_rngs[i],
                    push_fr_rng      = push_fr_rngs[i],
                    push_to_rng      = push_to_rngs[i],
                    pull_fr_rng_bnd  = pull_fr_rngs_bnd[i, :],
                    pull_to_rng_bnd  = pull_to_rngs_bnd[i, :],
                    push_fr_rng_bnd  = push_fr_rngs_bnd[i, :],
                    push_to_rng_bnd  = push_to_rngs_bnd[i, :],
                )
            end

    end

    # Restore `View` objects
    ocn.cols = ocn.cols
    ocn.lays = ocn.lays

end

function run!(
    ocn :: Ocean;
    Δt  :: Float64,
    substeps :: Integer,
    cfgs...
)

    (ocn.id == 0) || throw(ErrorException("`id` is not 0 (master). Id received: " * string(ocn.id)))

    dt = Δt / substeps

    cost_prep = @elapsed @sync for (i, p) in enumerate(wkrs)
        @spawnat p let
            syncForcingFromMaster!(subocn)
            cleanQflx2atm!(subocn.worker_ocn)
            stepOcean_prepare!(subocn.worker_ocn; cfgs...)
        end
    end

    sync_bnd_vars2 = (:T_ML, :S_ML, :h_ML, :FLDO)
    sync_bnd_vars3 = (:Ts,   :Ss)

    sync_to_master_vars2 = (:T_ML, :S_ML, :b_ML, :h_ML, :FLDO, :h_MO, :fric_u, :qflx2atm, :qflx2atm_pos, :qflx2atm_neg, :τx, :τy, :TSAS_clim, :SSAS_clim, :TFLUX_DIV_implied, :SFLUX_DIV_implied, :TEMP, :dTEMPdt, :SALT, :dSALTdt, :dTdt_ent, :dSdt_ent, :TFLUX_bot, :SFLUX_bot, :SFLUX_top, :seaice_nudge_energy)
    sync_to_master_vars3 = (:Ts, :Ss, :bs, :u, :v, :w_bnd, :TFLUX_CONV, :SFLUX_CONV, :TFLUX_DEN_z, :SFLUX_DEN_z, :div, :Ts_mixed, :Ss_mixed, :qflx_T_correction, :qflx_S_correction)

    #accumulative_vars2 = (:dTdt_ent, :dSdt_ent)
    #accumulative_vars3 = (:TFLUX_CONV, :T_vflux_ML, :SFLUX_CONV, :S_vflux_ML)

    cost_main = @elapsed for substep = 1:substeps
        #println("substep: ", substep)

        @sync for (i, p) in enumerate(wkrs)
            @spawnat p let
                syncBoundaryFromMaster!(subocn; vars3 = sync_bnd_vars3, vars2 = sync_bnd_vars2)
                calFLDOPartition!(subocn.worker_ocn)
                stepOcean_Flow!(subocn.worker_ocn; Δt = dt, cfgs...)
                stepOcean_MLDynamics!(subocn.worker_ocn; Δt = dt, cfgs...)
            end
        end

        @sync for (i, p) in enumerate(wkrs)
            @spawnat p let
                syncBoundaryToMaster!(subocn; vars3 = sync_bnd_vars3, vars2 = sync_bnd_vars2)
                accumulate!(subocn.worker_ocn)
            end
        end

    end

    cost_final = @elapsed @sync for (i, p) in enumerate(wkrs)
        @spawnat p let

            stepOcean_slowprocesses!(subocn.worker_ocn; Δt = Δt, cfgs...)

            if cfgs[:do_qflx_finding]
                calFlxCorrection!(subocn.worker_ocn; Δt = Δt, τ=10*86400.0, cfgs...)
            end

            if cfgs[:do_seaice_nudging]
                nudgeSeaice!(subocn.worker_ocn; Δt = Δt, τ=5*86400.0, cfgs...)
            end

            calLatentHeatReleaseOfFreezing!(subocn.worker_ocn; Δt=Δt, do_convadjust=cfgs[:do_convadjust])

            avg_accumulate!(subocn.worker_ocn; count=substeps)

            calDirect∂TEMP∂t!(subocn.worker_ocn; Δt=Δt)
            calDirect∂SALT∂t!(subocn.worker_ocn; Δt=Δt)
            calImplied∂TEMP∂t!(subocn.worker_ocn; cfgs...)
            calImplied∂SALT∂t!(subocn.worker_ocn; cfgs...)

            # Special: calculate Ts_mixed and Ss_mixed
            calTsSsMixed!(subocn.worker_ocn)

            syncToMaster!(
                subocn;
                vars2 = sync_to_master_vars2,
                vars3 = sync_to_master_vars3,
            )
        end
    end

    calDiagnostics!(ocn)

    println(format("### Cost: prep={:.1f}s , main={:.1f}s, final={:.1f}s. ###", cost_prep, cost_main, cost_final))

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
export NKLandscape

# Construct an NK lanscape, see:
# Kauffman, S. A. (1993). The Origins of Order: Self-Organization and Selection
# in Evolution. New York, New York, USA: Oxford University Press.
type NKLandscape <: Landscape
  n::Int64
  k::Int64
  a::Int64
  links::AlleleLinks
  contribs::AlleleContribs
end

function NKLandscape(n::Int64, k::Int64; a::Int64=2, near::Bool=false)
  if a != 2
    # TODO:  implement a values other than 2
    error("only a == 2 is implemented at this time")
  end
  links = makelinks(n, k, near)
  contribs = newAlleleContribs(n)
  return NKLandscape(n, k, a, links, contribs)
end

show(io::Base.IO, l::NKLandscape) = print(io, "NKLandscape($(l.n), $(l.k))")


text:
__precompile__(true)

module Randomshit

function l(x::Int, y::Int, n::Int)::Vector{Tuple{Int, Int}}
	mov::Vector{Tuple{Int, Int}} = [
		(x + 1, y + 2), (x - 1, y + 2), (x - 2, y + 1), (x - 2, y - 1),
		(x - 1, y - 2), (x + 1, y - 2), (x + 2, y - 1), (x + 2, y + 1)
	]

	return [i for i in mov if i[1] >= 1 && i[1] <= n && i[2] >= 1 && i[2] <= n]
end

function chess(n::Int, k::Int)::Array{Float64}
	p::Array{Float64, 3} = Array{Float64, 3}(undef, (n, n, k + 1))
	p[:, :, 1] .= 1

	for i in 2 : k + 1
		for x in 1 : n
			for y in 1 : n
				p[x, y, i] = (1 / 8) * sum([p[x_, y_, i - 1]
					for (x_, y_) in l(x, y, n)])
			end
		end
	end

	# return p[:, :, k + 1]
	return p[1, 1, :]
end

testsum() = 0

testsum(a, rest...) = a + testsum(rest...)

end


text:
using TermWin
if VERSION < v"0.4-"
    using Dates
else
    using Base.Dates
end

#TermWin.logstart()
TermWin.initsession()
#s = "a test"
v = newTwList( rootTwScreen; height=25, width=80, posy = :random, posx = :random, box=true, horizontal=true )
v1 = newTwList( v; height=25,width=80)
for i = 1:12
    a = newTwEntry( v1, Float64; width=30, title=@sprintf("float %2d", i), box=true)
    a.data.inputText = "0"
    a.value = 0
end

v2 = newTwList( v; height=25,width=80)
for i = 1:7
    a = newTwEntry( v2, UTF8String; width=30, title=@sprintf("string %2d",i), box=true)
    a.data.inputText = "abcdefghijk"
    a.value = "abcdefghijk"
end

v3 = newTwList( v )
for i=1:10
    a = newTwEntry( v3, Int; width=30, title=@sprintf("integer %2d",i), box=true)
    a.data.inputText = "1"
    a.value = 1
end
TermWin.update_list_canvas( v )

tshow( v )
#activateTwObj( rootTwScreen )
#unregisterTwObj( rootTwScreen, v )
#TermWin.endsession()
println( v )


text:
export pdf2xml
using JSON

mutable struct DeepFigure
end

function read_deepfigure(path::String)
end

function extract(texts::Vector{PDText}, page::Int, frect::Rectangle, crect::Rectangle)
    range = 0:0
    offset = 1
    while true
        i = findnext(texts,offset) do t
            t.page == page && contains(frect,t.fcoord) && !contains(crect,t.fcoord)
        end
        i == nothing && break
        j = findnext(texts, i) do t
            t.page == page && (!contains(frect,t.fcoord) || contains(crect,t.fcoord))
        end
        j -= 1
        length(i:j) > length(range) && (range = i:j)
        offset = j + 1
    end
    texts[range]
end

function pdf2xml(pdfpath::String)
    texts = readpdf(pdfpath)

    jsonpath = pdfpath[1:end-4] * "deepfigures-results.json"
    dict = JSON.parsefile(jsonpath)
    tables = filter(x -> x["figure_type"] == "Table", dict["figures"])
    pdtables = PDTable[]
    for t in tables
        page = t["page"] + 1
        caption = t["caption_text"]
        c = 1.38888888889
        b = t["figure_boundary"]
        x1, x2, y1, y2 = b["x1"], b["x2"], b["y1"], b["y2"]
        frect = Rectangle(x1/c, y1/c, (x2-x1)/c, (y2-y1)/c)

        b = t["caption_boundary"]
        x1, x2, y1, y2 = b["x1"], b["x2"], b["y1"], b["y2"]
        crect = Rectangle(x1/c, y1/c, (x2-x1)/c, (y2-y1)/c)

        tabletexts = extract(texts, page, frect, crect)
        table = PDTable(tabletexts, caption)
        push!(pdtables, table)
    end
    xhtml = toxhtml(pdtables)
    open("$pdfpath.xhtml","w") do io
        println(io, xhtml)
    end

    #table = PDTable(texts)
    #xhtml = toxhtml(table)

    #r = Rectangle
    #i1 = findfirst(t -> contains(r.fcoord,t.fcoord), texts)
    #i2 = findfirst(t)
    #tablerange = texts[i1:i2]
end