Allow the boundary operators of a ChainComplex to be abstract.Protopraphs

qldpc/abstract.py

@@ -152,6 +152,9 @@ def __eq__(self, other: object) -> bool:
             isinstance(other, Group) and self._field == other._field and self._group == other._group
         )
 
+    def __hash__(self) -> int:
+        return hash(self._group)
+
     def __contains__(self, member: GroupMember) -> bool:
         return comb.Permutation(member.array_form) in self._group
 
@@ -421,6 +424,9 @@ def __eq__(self, other: object) -> bool:
             and all(self._vec[member] == other._vec[member] for member in other._vec)
         )
 
+    def __bool__(self) -> bool:
+        return any(self._vec.values())
+
     def __iter__(self) -> Iterator[tuple[GroupMember, galois.FieldArray]]:
         yield from self._vec.items()
 
@@ -592,7 +598,9 @@ def T(self) -> Protograph:
 
     @classmethod
     def build(
-        cls, group: Group, array: npt.NDArray[np.object_] | Sequence[Sequence[object]]
+        cls,
+        group: Group,
+        array: npt.NDArray[np.object_ | np.int_] | Sequence[Sequence[object | int]],
     ) -> Protograph:
         """Construct a protograph.
 
@@ -603,7 +611,8 @@ def build(
         """
         array = np.array(array, dtype=object)
         vals = [Element(group, member) if member else Element(group) for member in array.ravel()]
-        return Protograph(np.array(vals, dtype=object).reshape(array.shape))
+        vals_array = np.array(vals, dtype=object).reshape(array.shape)
+        return Protograph(vals_array, group)
 
 
 ################################################################################

qldpc/abstract_test.py

@@ -35,6 +35,7 @@ def test_permutation_group() -> None:
     assert group.random() in group
     assert group.random(seed=0) == group.random(seed=0)
     assert group.to_sympy() == group._group
+    assert hash(group) == hash(group.to_sympy())
 
     assert abstract.Group.from_generating_mats([[1]]) == abstract.CyclicGroup(1)
 
@@ -106,6 +107,7 @@ def test_algebra() -> None:
     group = abstract.TrivialGroup(field=3)
     zero = abstract.Element(group)
     one = abstract.Element(group).one()
+    assert bool(one) and not bool(zero)
     assert zero.group == group
     assert one + 2 == group.identity + 2 * one == -one + 1 == one - 1 == zero
     assert group.identity * one == one * group.identity == one**2 == one

qldpc/codes.py

@@ -2222,4 +2222,6 @@ def __init__(
             chain = ChainComplex(*chain.ops[:2])
 
         matrix_x, matrix_z = chain.op(1), chain.op(2).T
+        assert not isinstance(matrix_x, abstract.Protograph)
+        assert not isinstance(matrix_z, abstract.Protograph)
         CSSCode.__init__(self, matrix_x, matrix_z, field, conjugate=conjugate, skip_validation=True)

qldpc/codes_test.py

@@ -222,7 +222,7 @@ def test_quantum_distance(field: int = 2) -> None:
 
 
 @pytest.mark.parametrize("field", [2, 3])
-def test_graph_product(
+def test_hypergraph_products(
     field: int,
     bits_checks_a: tuple[int, int] = (5, 3),
     bits_checks_b: tuple[int, int] = (3, 2),
@@ -257,10 +257,18 @@ def test_trivial_lift(
     protograph_b = abstract.TrivialGroup.to_protograph(code_b.matrix, field)
     code_LP = codes.LPCode(protograph_a, protograph_b)
 
-    assert np.array_equal(code_HGP.matrix, code_LP.matrix)
+    assert np.array_equal(code_HGP.matrix_x, code_LP.matrix_x)
+    assert np.array_equal(code_HGP.matrix_z, code_LP.matrix_z)
     assert nx.utils.graphs_equal(code_HGP.graph, code_LP.graph)
     assert np.array_equal(code_HGP.sector_size, code_LP.sector_size)
 
+    chain = objects.ChainComplex.tensor_product(protograph_a, protograph_b.T)
+    matrix_x, matrix_z = chain.op(1), chain.op(2).T
+    assert isinstance(matrix_x, abstract.Protograph)
+    assert isinstance(matrix_z, abstract.Protograph)
+    assert np.array_equal(matrix_x.lift(), code_HGP.matrix_x)
+    assert np.array_equal(matrix_z.lift(), code_HGP.matrix_z)
+
 
 def test_lift() -> None:
     """Verify lifting in Eqs. (8) and (10) of arXiv:2202.01702v3."""

qldpc/objects.py

@@ -320,7 +320,6 @@ def satisfies_total_no_conjugacy(
         )
 
 
-# TODO: investigate the lifted product for chain complexes: https://arxiv.org/pdf/2012.04068.pdf
 class ChainComplex:
     """Chain complex: a sequence modules with "boundary operators" that map between them.
 
@@ -342,19 +341,50 @@ class ChainComplex:
     """
 
     _field: type[galois.FieldArray]
-    _ops: tuple[galois.FieldArray, ...]
+    _ops: tuple[npt.NDArray[np.int_] | abstract.Protograph, ...]
 
-    def __init__(self, *ops: npt.NDArray[np.int_], field: int | None = None) -> None:
-        fields = set(type(op) for op in ops if isinstance(op, galois.FieldArray))
-        fields |= set([galois.GF(field)]) if field is not None else set()
-        if len(fields) > 1:
-            raise ValueError("Inconsistent base fields provided for chain complex")
+    # if boundary operators are defined over a group algebra, keep track of their base group
+    _group: abstract.Group | None
+
+    def __init__(
+        self,
+        *ops: npt.NDArray[np.int_] | abstract.Protograph,
+        field: int | None = None,
+        skip_validation: bool = False,
+    ) -> None:
+        # check that either all or none of the operators are defined over a group algebra
+        if not (
+            all(isinstance(op, abstract.Protograph) for op in ops)
+            or not any(isinstance(op, abstract.Protograph) for op in ops)
+        ):
+            raise ValueError("Invalid or inconsistent operator types provided for a ChainComplex")
+
+        # identify the base field and group for the boundary operators of this chain complex
+        fields = set([galois.GF(field)]) if field is not None else set()
+        groups = set()
+        for op in ops:
+            if isinstance(op, abstract.Protograph):
+                fields.add(op.field)
+                groups.add(op.group)
+            elif isinstance(op, galois.FieldArray):
+                fields.add(type(op))
+        if len(fields) > 1 or len(groups) > 1:
+            raise ValueError("Inconsistent base fields (or groups) provided for chain complex")
         self._field = fields.pop() if fields else galois.GF(DEFAULT_FIELD_ORDER)
+        self._group = groups.pop() if groups else None
+
+        # identify the boundary operators of this chain complex
+        if self._group is None:
+            self._ops = tuple(self.field(op) for op in ops)
+        else:
+            self._ops = ops
 
-        self._ops = tuple(self.field(op) for op in ops)
-        for degree in range(1, self.num_links):
-            op_a = self.op(degree)
-            op_b = self.op(degree + 1)
+        if not skip_validation:
+            self._validate_ops()
+
+    def _validate_ops(self) -> None:
+        """Validate the consistency of this the boundary operators in this chain complex."""
+        for op_a, op_b in zip(self.ops, self.ops[1:]):
             if op_a.shape[1] != op_b.shape[0] or np.any(op_a @ op_b):
                 raise ValueError(
                     "Condition for a chain complex not satisfied:\n"
@@ -367,20 +397,31 @@ def field(self) -> type[galois.FieldArray]:
         return self._field
 
     @property
-    def ops(self) -> tuple[npt.NDArray[np.int_], ...]:
-        """The boundary operators of this chain complex."""
-        return self._ops
+    def group(self) -> abstract.Group | None:
+        """The base group of this chain complex."""
+        return self._group
 
     @property
     def num_links(self) -> int:
         """The number of "internal" links in this chain complex."""
         return len(self.ops)
 
+    @property
+    def ops(self) -> tuple[npt.NDArray[np.int_] | abstract.Protograph, ...]:
+        """The boundary operators of this chain complex."""
+        return self._ops
+
     def dim(self, degree: int) -> int:
         """The dimension of the module of the given degree."""
         return self.op(degree).shape[1]
 
-    def op(self, degree: int) -> npt.NDArray[np.int_]:
+    @property
+    def T(self) -> ChainComplex:
+        """Transpose and reverse the order of the boundary operators in this chain complex."""
+        dual_ops = [op.T for op in self.ops[::-1]]
+        return ChainComplex(*dual_ops, skip_validation=True)
+
+    def op(self, degree: int) -> npt.NDArray[np.int_] | abstract.Protograph:
         """The boundary operator of this chain complex that acts on the module of a given degree."""
         assert 0 <= degree <= self.num_links + 1
         if degree == 0:
@@ -389,17 +430,11 @@ def op(self, degree: int) -> npt.NDArray[np.int_]:
             return self.field.Zeros((self.ops[-1].shape[1], 0))
         return self.ops[degree - 1]
 
-    @property
-    def T(self) -> ChainComplex:
-        """Transpose and reverse the order of the boundary operators in this chain complex."""
-        dual_ops = [op.T for op in self.ops[::-1]]
-        return ChainComplex(*dual_ops)
-
     @classmethod
     def tensor_product(  # noqa: C901 ignore complexity check
         cls,
-        chain_a: ChainComplex | galois.FieldArray | npt.NDArray[np.int_],
-        chain_b: ChainComplex | galois.FieldArray | npt.NDArray[np.int_],
+        chain_a: ChainComplex | npt.NDArray[np.int_] | galois.FieldArray | abstract.Protograph,
+        chain_b: ChainComplex | npt.NDArray[np.int_] | galois.FieldArray | abstract.Protograph,
         field: int | None = None,
     ) -> ChainComplex:
         """Tensor product of two chain complexes.
@@ -429,8 +464,8 @@ def tensor_product(  # noqa: C901 ignore complexity check
             chain_a = ChainComplex(chain_a, field=field)
         if not isinstance(chain_b, ChainComplex):
             chain_b = ChainComplex(chain_b, field=field)
-        if chain_a.field is not chain_b.field:
-            raise ValueError("Cannot take tensor product of chain complexes over different fields")
+        if chain_a.field is not chain_b.field or chain_a.group != chain_b.group:
+            raise ValueError("Incompatible chain complexes: different base fields or groups")
         chain_field = chain_a.field
 
         def get_degree_pairs(degree: int) -> Iterator[tuple[int, int]]:
@@ -455,7 +490,7 @@ def get_zero_block(
             cols = chain_a.dim(col_deg_a) * chain_b.dim(col_deg_b)
             return chain_field.Zeros((rows, cols))
 
-        ops: list[npt.NDArray[np.int_]] = []
+        ops = []
         for degree in range(1, chain_a.num_links + chain_b.num_links + 1):
             # fill in zero blocks of the total boundary operator
             blocks = [
@@ -470,12 +505,18 @@ def get_zero_block(
                 if deg_a:
                     row = get_block_index(deg_a - 1, deg_b)
                     iden_b = np.identity(op_b.shape[1], dtype=op_b.dtype)
-                    blocks[row][col] = np.kron(op_a, iden_b)
+                    blocks[row][col] = np.kron(op_a, iden_b)  # type:ignore[assignment,arg-type]
                 if deg_b:
                     row = get_block_index(deg_a, deg_b - 1)
                     iden_a = np.identity(op_a.shape[1], dtype=op_a.dtype)
-                    blocks[row][col] = np.kron(iden_a, op_b) * (-1) ** deg_a
+                    blocks[row][col] = (
+                        np.kron(iden_a, op_b) * (-1) ** deg_a  # type:ignore[arg-type]
+                    )
 
             ops.append(np.block(blocks))
 
-        return ChainComplex(*ops, field=chain_field.order)
+        if chain_a.group is None:
+            ops = [chain_field(op) for op in ops]
+        else:
+            ops = [abstract.Protograph(op) for op in ops]
+        return ChainComplex(*ops, skip_validation=True)

qldpc/objects_test.py

@@ -125,12 +125,21 @@ def test_chain_complex(field: int = 3) -> None:
     assert not np.any(two_chain.op(two_chain.num_links + 1))
 
     # tensor product of a two-complex and its dual
-    objects.ChainComplex.tensor_product(two_chain, two_chain.T, field)
+    four_chain = objects.ChainComplex.tensor_product(two_chain, two_chain.T, field)
+    four_chain._validate_ops()
+
+    # tensor product of one-complexes over a group algebra
+    protograph = abstract.Protograph.build(abstract.TrivialGroup(field), mat)
+    two_chain = objects.ChainComplex.tensor_product(protograph, protograph, field)
+    assert not np.any(two_chain.op(0))
+    assert not np.any(two_chain.op(two_chain.num_links + 1))
 
     # invalid chain complex constructions
+    with pytest.raises(ValueError, match="inconsistent operator types"):
+        objects.ChainComplex(mat, abstract.TrivialGroup.to_protograph([[0]]))
     with pytest.raises(ValueError, match="Inconsistent base fields"):
         objects.ChainComplex(galois.GF(field)(mat), field=field**2)
     with pytest.raises(ValueError, match="boundary operators .* must compose to zero"):
         objects.ChainComplex(mat, mat, field=field)
-    with pytest.raises(ValueError, match="different fields"):
+    with pytest.raises(ValueError, match="different base fields"):
         objects.ChainComplex.tensor_product(galois.GF(field)(mat), galois.GF(field**2)(mat))