Implement balanced product codes

qldpc/abstract.py

@@ -88,17 +88,25 @@ def __matmul__(self, other: GroupMember) -> GroupMember:
         """
         return GroupMember(self.array_form + [val + self.size for val in other.array_form])
 
-    def to_matrix(self) -> npt.NDArray[np.int_]:
+    def to_matrix(self, size: int | None = None) -> npt.NDArray[np.int_]:
         """Lift this permutation object to a permutation matrix.
 
         For consistency with how SymPy composes permutations, the permutation matrix constructed
         here is right-acting, meaning that it acts on a vector v as v --> v @ p.to_matrix().  This
         convension ensures that this lift is a homomorphism on SymPy Permutation objects, which is
         to say that (p * q).to_matrix() = p.to_matrix() @ q.to_matrix().
+
+        Size specifies to final size of the output matrix. Must be strictly larger than original
+        size of the permutation. Each additional element is treated as a fixed point
         """
         matrix = np.zeros([self.size] * 2, dtype=int)
         for ii in range(self.size):
             matrix[ii, self.apply(ii)] = 1
+
+        if size is not None:
+            if size > len(matrix):
+                matrix = scipy.linalg.block_diag(matrix, np.eye(size - len(matrix), dtype=int))
+
         return matrix
 
     def to_gap_cycles(self) -> str:

qldpc/abstract_test.py

@@ -25,7 +25,9 @@
 import numpy as np
 import numpy.typing as npt
 import pytest
+from sympy.combinatorics import Permutation
 
+import qldpc
 from qldpc import abstract
 
 
@@ -383,3 +385,12 @@ def test_small_group() -> None:
         group = abstract.SmallGroup(1, 1)
     assert group == abstract.TrivialGroup()
     assert group.random() == group.identity
+
+
+def test_to_matrix() -> None:
+    perm = qldpc.abstract.GroupMember.from_sympy(Permutation(0, 2, 1))
+    print("\n", perm.to_matrix())
+    perm = qldpc.abstract.GroupMember.from_sympy(Permutation(3)(0, 2, 1))
+    print("\n", perm.to_matrix())
+
+    ...

qldpc/codes/__init__.py

@@ -22,8 +22,10 @@
 from .quantum import (
     BaconShorCode,
     BBCode,
+    BPCode,
     C4Code,
     C6Code,
+    CCode,
     FiveQubitCode,
     GeneralizedSurfaceCode,
     HGPCode,
@@ -56,8 +58,10 @@
     "QuditCode",
     "BaconShorCode",
     "BBCode",
+    "BPCode",
     "C4Code",
     "C6Code",
+    "CCode",
     "FiveQubitCode",
     "GeneralizedSurfaceCode",
     "HGPCode",

qldpc/codes/quantum.py

@@ -31,6 +31,7 @@
 import scipy
 import sympy
 
+import qldpc.external.groups
 from qldpc import abstract
 from qldpc.abstract import DEFAULT_FIELD_ORDER
 from qldpc.math import first_nonzero_cols
@@ -1564,6 +1565,119 @@ def __init__(
         CSSCode.__init__(self, matrix_x, matrix_z, field)
 
 
+class BPCode(CSSCode):
+    """Code created from the product of classical codes. Similar to hypergraph codes.
+
+    Binary matrix (M) must have 2 permutations, r and c. Where r*M = M*C^T and r^l = I and c^l for some l.
+
+    qldpc.external.groups.get_balanced_permutations_of_matrix can be used to find these permutations.
+
+    I.e.
+    matrix = np.array(
+        [
+            [1, 1, 0, 0, 0, 0],
+            [0, 1, 1, 0, 0, 0],
+            [0, 0, 1, 1, 0, 0],
+            [0, 0, 0, 1, 1, 0],
+            [0, 0, 0, 0, 1, 1],
+            [1, 0, 0, 0, 0, 1],
+        ]
+    )
+    R,C = ldpc.external.groups.get_balanced_permutations_of_matrix(matrix,3)
+    balanced_code = codes.BalancedProductCode(matrix,R,C)
+
+    References:
+    - https://errorcorrectionzoo.org/c/balanced_product
+    - https://arxiv.org/pdf/2505.13679
+
+    """
+
+    def __init__(
+        self,
+        seed_matrix: np.typing.NDArray[np.int_],
+        R: qldpc.abstract.GroupMember,
+        C: qldpc.abstract.GroupMember,
+        field: int | None = None,
+    ) -> None:
+        if not np.all((seed_matrix == 0) | seed_matrix == 1):
+            raise ValueError("Binary matrix must only have zeros or ones")
+
+        if seed_matrix.shape[0] != seed_matrix.shape[1]:
+            raise ValueError("Only square binary matrix is supported at this time")
+        if not np.array_equal(
+            R.to_matrix(seed_matrix.shape[0]) @ seed_matrix,
+            seed_matrix @ C.to_matrix(seed_matrix.shape[1]).T,
+        ):
+            raise ValueError("Invalid permutations provided")
+
+        self._r = R
+        self._c = C
+        self._seed_matrix = seed_matrix
+
+        check_x = np.hstack(
+            (seed_matrix.T, np.eye(seed_matrix.shape[0], dtype=int) + self._c.to_matrix())
+        )
+        check_z = np.hstack(
+            (np.eye(seed_matrix.shape[1], dtype=int) + self._r.to_matrix(), seed_matrix)
+        )
+        CSSCode.__init__(self, check_x, check_z, field, is_subsystem_code=False)
+
+    @classmethod
+    def from_codes(
+        cls, code_q: CSSCode, code_c: ClassicalCode, distance_balancing: bool = False
+    ) -> qldpc.codes.BPCode:
+        """
+        Allows for the creation of balanced codes from a classical code and a CSS code.
+
+        https://arxiv.org/pdf/2505.13679v1
+
+        If 'distance_balancing', function will determine if x and z errors have the same distance
+        and if not, will construct the code to maximize the final distance. This involve exact distance calculations
+        which will only work on relatively small codes
+
+        """
+
+        # i.e. rows of matrix, number qubits of matrix
+        r_C, n_C = code_c.matrix.shape
+
+        high_distance = code_q.code_z
+        low_distance = code_q.code_x
+        if distance_balancing:
+            if code_q.code_x.get_distance_exact() > code_q.code_z.get_distance_exact():
+                temp = high_distance
+                high_distance = low_distance
+                low_distance = temp
+
+        r_low, n_low = low_distance.matrix.shape
+
+        # Make new H_x
+        Hx_tensor = np.kron(low_distance.matrix, np.eye(n_C, dtype=int))
+        classical_T_tensor = np.kron(
+            np.eye(r_low, dtype=int), code_c.matrix.T
+        )  # shape: (r_low * n_C, r_low * r_C)
+        H_new_X = np.hstack((Hx_tensor, classical_T_tensor))
+
+        # Make new H_z
+        Hz_tensor = np.kron(high_distance.matrix, np.eye(n_C, dtype=int))
+
+        bottom_left = np.kron(np.eye(n_low, dtype=int), code_c.matrix)
+        bottom_right = np.kron(low_distance.matrix.T, np.eye(r_C, dtype=int))
+
+        top_right_zeros = np.zeros(
+            (Hz_tensor.shape[0], bottom_left.shape[1] + bottom_right.shape[1] - Hz_tensor.shape[1]),
+            dtype=int,
+        )
+        top_row = np.hstack((Hz_tensor, top_right_zeros))
+
+        bottom_row = np.hstack((bottom_left, bottom_right))
+
+        H_new_Z = np.vstack((top_row, bottom_row))
+
+        new_code = cls.__new__(cls)
+        CSSCode.__init__(new_code, H_new_X, H_new_Z)
+        return new_code
+
+
 class BaconShorCode(SHPCode):
     """Bacon-Shor code on a square grid, implemented as a subsystem hypergraph product code.
 
@@ -1614,3 +1728,86 @@ def __init__(
         self._dimension = dim_x * dim_z
         self._distance_x = code_z.get_distance()  # X errors are witnessed by the Z code
         self._distance_z = code_x.get_distance()  # Z errors are witnessed by the X code
+
+class CCode(CSSCode):
+    def __init__(
+            self,
+            colors: list[int],
+            plaquettes: list[list[list[int]]],
+            dimension : int
+    ) -> None:
+        if len(colors) != len(plaquettes):
+            raise ValueError("Every plaquette needs a color")
+        unique_colors = set()
+        for color in colors:
+            unique_colors.add(color)
+        if len(unique_colors) <= dimension:
+            raise ValueError("N-colorable, may only have n colors")
+        unique_qubits =set() # do i really need this
+        for index_f,f in enumerate(plaquettes):
+            for index_g,g in enumerate(plaquettes):
+                for edge_f,edge_g in itertools.product(f, g):
+                    if edge_f == edge_g and colors[index_f] == colors[index_g] and index_f != index_g:
+                        raise ValueError("Invalid coloring")
+                    for qubit in edge_f + edge_g:
+                        unique_qubits.add(qubit)
+        #each generator, loop through qubit add identity if in, otherwise not, find ordering of
+        ##add
+        code_x = []
+        # code_z = []
+        for plaquette in plaquettes:
+            qubits = set()
+            for edge in plaquette:
+                for qubit in edge:
+                    qubits.add(qubit)
+            stabilizer = []
+            for bit in range(0,len(unique_qubits)):
+                if bit in qubits:
+                    stabilizer.append(1)
+                else:
+                    stabilizer.append(0)
+                code_x.append(stabilizer)
+
+        CSSCode.__init__(self,code_x,code_x)
+
+#
+# ###
+#   r_C, n_C = code_c.matrix.shape
+#
+#         high_distance = code_q.code_z
+#         low_distance = code_q.code_x
+#         if distance_balancing:
+#             if code_q.code_x.get_distance_exact() > code_q.code_z.get_distance_exact():
+#                 temp = high_distance
+#                 high_distance = low_distance
+#                 low_distance = temp
+#
+#         r_low, n_low = low_distance.matrix.shape
+#
+#         # Make new H_x
+#         Hx_tensor = np.kron(low_distance.matrix, np.eye(n_C, dtype=int))
+#         classical_T_tensor = np.kron(
+#             np.eye(r_low, dtype=int), code_c.matrix.T
+#         )  # shape: (r_low * n_C, r_low * r_C)
+#         H_new_X = np.hstack((Hx_tensor, classical_T_tensor))
+#
+#         # Make new H_z
+#         Hz_tensor = np.kron(high_distance.matrix, np.eye(n_C, dtype=int))
+#
+#         bottom_left = np.kron(np.eye(n_low, dtype=int), code_c.matrix)
+#         bottom_right = np.kron(low_distance.matrix.T, np.eye(r_C, dtype=int))
+#
+#         top_right_zeros = np.zeros(
+#             (Hz_tensor.shape[0], bottom_left.shape[1] + bottom_right.shape[1] - Hz_tensor.shape[1]),
+#             dtype=int,
+#         )
+#         top_row = np.hstack((Hz_tensor, top_right_zeros))
+#
+#         bottom_row = np.hstack((bottom_left, bottom_right))
+#
+#         H_new_Z = np.vstack((top_row, bottom_row))
+#
+#         new_code = cls.__new__(cls)
+#         CSSCode.__init__(new_code, H_new_X, H_new_Z)
+#         return new_code
+