Feature: Added SHGP and SHYPS Subsystem Quantum Codes

docs/src/references.bib

@@ -1220,3 +1220,28 @@ @article{bravyi2019simulation
   year={2019},
   doi={10.22331/q-2019-09-02-181}
 }
+
+@article{bravyi2011subsystem,
+  title={Subsystem codes with spatially local generators},
+  author={Bravyi, Sergey},
+  journal={Physical Review A},
+  volume={83},
+  number={1},
+  pages={012320},
+  year={2011},
+  publisher={APS}
+}
+
+@article{li2020numerical,
+  title={A Numerical Study of Bravyi-Bacon-Shor and Subsystem Hypergraph Product Codes},
+  author={Li, Muyuan and Yoder, Theodore J},
+  journal={arXiv preprint arXiv:2002.06257},
+  year={2020}
+}
+
+@article{malcolm2025computing,
+  title={Computing Efficiently in QLDPC Codes},
+  author={Malcolm, Alexander J. and Glaudell, Andrew N. and Fuentes, Patricio and Chandra, Daryus and Schotte, Alexis and DeLisle, Colby and Haenel, Rafael and Ebrahimi, Amir and Roffe, Joschka and Quintavalle, Armanda O. and Beale, Stefanie J. and Lee-Hone, Nicholas R. and Simmons, Stephanie},
+  journal={arXiv preprint arXiv:2502.07150},
+  year={2025}
+}

lib/QECCore/ext/QECCoreNemoExt/QECCoreNemoExt.jl

@@ -1,13 +1,15 @@
-module QECCoreNemoExt
+module QECCoreNemoExt
 
 using QECCore
 using DocStringExtensions
+using SparseArrays: sparse
 
 import Nemo
 import Nemo: GF, gen, matrix, rank, transpose, polynomial_ring, evaluate, FqFieldElem,
-    FqPolyRingElem, degree, is_irreducible, gcd, derivative, inv, coeff, is_monic, one
+    FqPolyRingElem, degree, is_irreducible, gcd, derivative, inv, coeff, is_monic, one,
+    ZZ, residue_ring, matrix_space, nullspace
 
-import QECCore: code_k, parity_matrix_x, parity_matrix_z, parity_matrix, generator_polynomial
+import QECCore: code_k, code_n, parity_matrix_x, parity_matrix_z, parity_matrix, generator_polynomial, distance, dual
 
 import Random
 import Random: MersenneTwister, GLOBAL_RNG, AbstractRNG, rand
@@ -31,5 +33,6 @@ function QECCore.code_k(c::AbstractCECC)
 end
 
 include("goppa.jl")
+include("simplex.jl")
 
 end # module

lib/QECCore/ext/QECCoreNemoExt/simplex.jl

@@ -0,0 +1,66 @@
+"""
+    $TYPEDEF
+
+The `[2ʳ - 1, r, 2ʳ⁻¹]` simplex code family, dual to the binary Hamming codes.
+
+`C(r)` is the dual of `Hamming(r)`. Its codewords are the rows of the Hamming
+parity check matrix, and every nonzero codeword has weight `2ʳ⁻¹`.
+
+Used as the seed code in the SHYPS construction [malcolm2025computing](@cite).
+ECC Zoo: [Simplex code family](https://errorcorrectionzoo.org/c/simplex).
+
+### Fields
+    \$TYPEDFIELDS
+"""
+struct Simplex <: AbstractCECC
+    """Parameter `r` (must be ≥ 2). Determines code length `n = 2ʳ - 1`."""
+    r::Int
+    function Simplex(r)
+        if r < 2
+            throw(ArgumentError("Invalid parameters: `r` must be ≥ 2 to obtain a valid code."))
+        end
+        new(r)
+    end
+end
+
+"""
+    dual(H)
+
+Compute the dual code of a binary parity check matrix `H` using Nemo's nullspace.
+Returns the parity check matrix of the dual code as a transposed nullspace matrix.
+
+This is a general-purpose utility: for any binary matrix `H`, the dual code's
+generator matrix `G` satisfies `H * Gᵀ = 0` over GF(2).
+"""
+function QECCore.dual(H)
+    H_nemo = matrix(GF(2), H)
+    null = Nemo.nullspace(H_nemo)[2]
+    @assert all(iszero, H_nemo * null)
+    return Nemo.transpose(null)
+end
+
+function QECCore.parity_matrix(c::Simplex)
+    r = c.r
+    n = 2^r - 1
+    # Building the Hamming parity check matrix -- column j is the number j in r-bit binary
+    H_hamming = zeros(Int, r, n)
+    for j in 1:n, i in 1:r
+        H_hamming[i, j] = (j >> (r - i)) & 1
+    end
+    # The dual of the Hamming code is the Simplex code
+    dual_mat = QECCore.dual(H_hamming)
+    # Converting Nemo matrix back to sparse Int matrix
+    nr, nc = size(dual_mat)
+    rows_idx = Int[]
+    cols_idx = Int[]
+    for i in 1:nr, j in 1:nc
+        if !iszero(dual_mat[i, j])
+            push!(rows_idx, i)
+            push!(cols_idx, j)
+        end
+    end
+    return sparse(rows_idx, cols_idx, ones(Int, length(rows_idx)), nr, n)
+end
+
+QECCore.code_n(c::Simplex) = 2^c.r - 1
+QECCore.distance(c::Simplex) = 2^(c.r - 1)

lib/QECCore/src/QECCore.jl

@@ -23,16 +23,17 @@ DelfosseReichardt, DelfosseReichardtRep, DelfosseReichardt823, QuantumTannerGrap
 TillichZemor, random_TillichZemor_code, BivariateBicycleViaCirculantMat
 
 # Classical Codes
-export RepCode, ReedMuller, RecursiveReedMuller, Golay, Hamming, random_Gallager_ldpc, Goppa, random_Goppa_code
+export RepCode, ReedMuller, RecursiveReedMuller, Golay, Hamming, Simplex, random_Gallager_ldpc, Goppa, random_Goppa_code
 
 # utilities
-export search_self_orthogonal_rm_code, hgp
+export search_self_orthogonal_rm_code, hgp, dual
 
 include("interface.jl")
 include("codes/util.jl")
 
 # Classical Codes
 include("codes/classical/hamming.jl")
+include("codes/classical/simplex.jl")
 include("codes/classical/repetition.jl")
 include("codes/classical/golay.jl")
 include("codes/classical/gallager.jl")

lib/QECCore/src/codes/classical/simplex.jl

@@ -0,0 +1,19 @@
+"""
+    dual(H)
+
+Compute the dual code of a binary parity check matrix `H`.
+Returns the parity check matrix of the dual code.
+
+Requires `Nemo` to be loaded.
+"""
+function dual end
+
+"""Simplex code, the dual of the Hamming code [malcolm2025computing](@cite)."""
+function Simplex(args...; kwargs...)
+    ext = Base.get_extension(QECCore, :QECCoreNemoExt)
+    if isnothing(ext)
+        throw("The `Simplex` code depends on the package `Nemo` but you have not installed or imported it yet. Immediately after you import `Nemo`, the `Simplex` will be available.")
+    end
+    return ext.Simplex(args...; kwargs...)
+end
+

src/ecc/ECC.jl

@@ -4,19 +4,20 @@ using QECCore
 import QECCore: code_n, code_s, code_k, rate, distance, parity_matrix_x, parity_matrix_z, parity_matrix,
 metacheck_matrix_x, metacheck_matrix_z, metacheck_matrix, hgp, generator_polynomial, hasmetachecks
 using QuantumClifford: QuantumClifford, AbstractOperation, AbstractStabilizer,
-    AbstractTwoQubitOperator, Stabilizer, PauliOperator,
+    AbstractTwoQubitOperator, Stabilizer, Tableau, PauliOperator,
     random_brickwork_clifford_circuit, random_all_to_all_clifford_circuit,
     canonicalize!, canonicalize_gott!,
     logicalxview, logicalzview, stabilizerview, destabilizerview, tab, phases,
     sCNOT, sSWAP, sHadamard, sPhase, sInvPhase,
     sZCX, sZCY, sZCZ, sXCX, sXCY, sXCZ, sYCX, sYCY, sYCZ, sZ, sX, sY, sMRZ, sMRX,
-    single_x, single_y, single_z, random_pauli!, PauliError,
+    single_x, single_y, single_z, random_pauli!, PauliError, xbit, zbit,
     apply!, comm, comm!, stab_to_gf2, embed, @S_str, affectedqubits, affectedbits,
     pftrajectories, measurements, mctrajectories
 import QuantumClifford: Stabilizer, MixedDestabilizer, nqubits
 
 using Combinatorics: combinations
 using LinearAlgebra: LinearAlgebra, I, rank, tr
+import Nemo
 using Nemo: ZZ, residue_ring, matrix, finite_field, GF, minpoly, coeff, lcm, FqPolyRingElem, FqFieldElem, is_zero, degree, defining_polynomial, is_irreducible, echelon_form
 using SparseArrays: sparse
 using Statistics: std
@@ -44,7 +45,9 @@ export parity_checks, parity_matrix_x, parity_matrix_z, iscss,
     GeneralizedBicycle, ExtendedGeneralizedBicycle,
     HomologicalProduct, DoubleHomologicalProduct,
     GeneralizedToric, TrivariateTricycle, BivariateBicycleViaPoly,
-    MultivariateMulticycle,
+    MultivariateMulticycle, SubsystemHypergraphProduct,
+    SubsystemHypergraphProductSimplex,
+    gauge_generators, code_g,
     evaluate_decoder,
     CommutationCheckECCSetup, NaiveSyndromeECCSetup, ShorSyndromeECCSetup,
     TableDecoder, CSSTableDecoder,
@@ -403,6 +406,59 @@ include("codes/concat.jl")
 include("codes/random_circuit.jl")
 include("codes/classical/bch.jl")
 
+"""Gauge group generators of a subsystem code, returned as a `Tableau`."""
+function gauge_generators end
+
+"""
+The number of gauge qubits in a subsystem code.
+
+For a subsystem code with `n` physical qubits, `k` logical qubits, and `s`
+independent stabilizer generators: `code_g(c) = n - k - s`.
+
+See also: [`code_n`](@ref), [`code_k`](@ref), [`gauge_generators`](@ref)
+"""
+function code_g end
+
+"""Compute the rank (number of independent generators) of a code's stabilizer group."""
+function _stabilizer_rank(c)
+    stab = parity_checks(c)
+    _, _, r = canonicalize!(Base.copy(stab), ranks=true)
+    return r
+end
+
+"""Extract X-only stabilizer rows from a mixed stabilizer tableau as a binary matrix."""
+function _stab_to_parity_x(stab::Stabilizer)
+    n = nqubits(stab)
+    rows = BitVector[]
+    for p in stab
+        xb, zb = xbit(p), zbit(p)
+        if any(xb) && !any(zb)
+            push!(rows, xb)
+        end
+    end
+    isempty(rows) && return falses(0, n)
+    return reduce(vcat, transpose.(rows))
+end
+
+"""Extract Z-only stabilizer rows from a mixed stabilizer tableau as a binary matrix."""
+function _stab_to_parity_z(stab::Stabilizer)
+    n = nqubits(stab)
+    rows = BitVector[]
+    for p in stab
+        xb, zb = xbit(p), zbit(p)
+        if !any(xb) && any(zb)
+            push!(rows, zb)
+        end
+    end
+    isempty(rows) && return falses(0, n)
+    return reduce(vcat, transpose.(rows))
+end
+
+# Subsystem codes
+
+include("codes/subsystem_hypergraph_product.jl")
+include("codes/subsystem_shyps.jl")
+
 # qLDPC
 include("codes/classical/lifted.jl")
 include("codes/qeccs_from_extensions.jl")

src/ecc/codes/subsystem_hypergraph_product.jl

@@ -0,0 +1,135 @@
+"""
+    $TYPEDEF
+
+A Subsystem Hypergraph Product (SHP) code built from two classical parity
+check matrices `H1` (m₁ × n₁) and `H2` (m₂ × n₂).
+
+Physical qubits: `n = n₁ · n₂`. Logical qubits: `k = nullity(H1) · nullity(H2)`.
+X-type gauge generators are `H1 ⊗ I`; Z-type are `I ⊗ H2`. Stabilizers:
+`S_X = H1 ⊗ ker(H2)` and `S_Z = ker(H1) ⊗ H2`.
+
+Based on [li2020numerical](@cite).
+ECC Zoo: [Subsystem product code family](https://errorcorrectionzoo.org/c/subsystem_product).
+
+See also: [`SubsystemHypergraphProductSimplex`](@ref)
+
+### Fields
+    $TYPEDFIELDS
+"""
+struct SubsystemHypergraphProduct <: AbstractCSSCode
+    H1::Matrix{Int}
+    H2::Matrix{Int}
+    gauge_generators::Tableau
+    stabilizer::Stabilizer
+end
+
+function SubsystemHypergraphProduct(H1::AbstractMatrix, H2::AbstractMatrix)
+    m1, n1 = size(H1)
+    m2, n2 = size(H2)
+    N = n1 * n2
+
+    gauge_ops = PauliOperator[]
+
+    # X gauge generators = H1 ⊗ I -- each row of H1 gives one X-type generator acting on n2 qubits
+    I_n2 = Matrix{Int}(I, n2, n2)
+    GX = kron(H1, I_n2) .% 2
+    for r in 1:size(GX, 1)
+        p = PauliOperator(0x0, falses(N), falses(N))
+        empty = true
+        for c in 1:N
+            if GX[r, c] == 1
+                p[c] = (true, false)
+                empty = false
+            end
+        end
+        if !empty push!(gauge_ops, p) end
+    end
+
+    # Z gauge generators = I ⊗ H2 -- same idea, each row of H2 gives one Z-type generator
+    I_n1 = Matrix{Int}(I, n1, n1)
+    GZ = kron(I_n1, H2) .% 2
+    for r in 1:size(GZ, 1)
+        p = PauliOperator(0x0, falses(N), falses(N))
+        empty = true
+        for c in 1:N
+            if GZ[r, c] == 1
+                p[c] = (false, true)
+                empty = false
+            end
+        end
+        if !empty push!(gauge_ops, p) end
+    end
+
+    gauge_generators = Tableau(gauge_ops)
+
+    F = GF(2)
+
+    # G1 = right nullspace of H1, i.e. vectors v where H1*v = 0
+    # Nemo.kernel gives left kernel by default (K*M=0), so we transpose H1 to get what we want
+    K1 = Nemo.kernel(transpose(matrix(F, H1)))
+    G1 = zeros(Int, Nemo.nrows(K1), n1)
+    for i in 1:Nemo.nrows(K1), j in 1:n1
+        G1[i,j] = K1[i,j] == F(1) ? 1 : 0
+    end
+
+    # G2 = right nullspace of H2 -- same thing as above but for H2
+    K2 = Nemo.kernel(transpose(matrix(F, H2)))
+    G2 = zeros(Int, Nemo.nrows(K2), n2)
+    for i in 1:Nemo.nrows(K2), j in 1:n2
+        G2[i,j] = K2[i,j] == F(1) ? 1 : 0
+    end
+
+    stabs = PauliOperator[]
+
+    # X stabilizers = H1 ⊗ G2 -- tensor each row of H1 with each row of G2 (nullspace of H2)
+    SX = kron(H1, G2) .% 2
+    for r in 1:size(SX, 1)
+        p = PauliOperator(0x0, falses(N), falses(N))
+        empty = true
+        for c in 1:N
+            if SX[r, c] == 1
+                p[c] = (true, false)
+                empty = false
+            end
+        end
+        if !empty push!(stabs, p) end
+    end
+
+    # Z stabilizers = G1 ⊗ H2 -- same idea, tensor nullspace of H1 with each row of H2
+    SZ = kron(G1, H2) .% 2
+    for r in 1:size(SZ, 1)
+        p = PauliOperator(0x0, falses(N), falses(N))
+        empty = true
+        for c in 1:N
+            if SZ[r, c] == 1
+                p[c] = (false, true)
+                empty = false
+            end
+        end
+        if !empty push!(stabs, p) end
+    end
+
+    stabilizer = Stabilizer(stabs)
+
+    return SubsystemHypergraphProduct(Matrix{Int}(H1), Matrix{Int}(H2), gauge_generators, stabilizer)
+end
+
+parity_checks(c::SubsystemHypergraphProduct) = c.stabilizer
+
+gauge_generators(c::SubsystemHypergraphProduct) = c.gauge_generators
+
+function code_k(c::SubsystemHypergraphProduct)
+    F = GF(2)
+    n1 = size(c.H1, 2)
+    n2 = size(c.H2, 2)
+    r1 = Nemo.rank(matrix(F, c.H1))
+    r2 = Nemo.rank(matrix(F, c.H2))
+    return (n1 - r1) * (n2 - r2)
+end
+
+function code_g(c::SubsystemHypergraphProduct)
+    return code_n(c) - code_k(c) - _stabilizer_rank(c)
+end
+
+parity_matrix_x(c::SubsystemHypergraphProduct) = _stab_to_parity_x(c.stabilizer)
+parity_matrix_z(c::SubsystemHypergraphProduct) = _stab_to_parity_z(c.stabilizer)