Fix #464: crash of Pauli-flow finding algorithm due to incorrect numba signatures.

CHANGELOG.md

@@ -34,6 +34,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - #458: The type for the `input_state` parameter in the simulator now allows `None` to indicate that the input qubits have already been specified in the backend.
 
+- #465: Fix #464. Functions in `graphix._linalg.py` convert arrays to c-contiguous form before passing them to numba-jitted functions.
+
 ### Changed
 
 - #181, #423: Structural separation of Pauli measurements

graphix/_linalg.py

@@ -31,7 +31,9 @@ def __new__(cls, data: npt.ArrayLike, copy: bool = True) -> Self:
             MatGF2
         """
         arr = np.array(data, dtype=np.uint8, copy=copy)
-        return super().__new__(cls, shape=arr.shape, dtype=arr.dtype, buffer=arr)
+        return super().__new__(
+            cls, shape=arr.shape, dtype=arr.dtype, buffer=arr, strides=arr.strides
+        )  # `strides=arr.strides` allows to preserve the memory layout of the original data.
 
     def mat_mul(self, other: MatGF2 | npt.NDArray[np.uint8]) -> MatGF2:
         r"""Multiply two matrices.
@@ -60,7 +62,7 @@ def mat_mul(self, other: MatGF2 | npt.NDArray[np.uint8]) -> MatGF2:
                 f"Dimension mismatch. Attempted to multiply `self` with shape {self.shape} and  `other` with shape {other.shape}"
             )
 
-        return MatGF2(_mat_mul_jit(self, other), copy=False)
+        return MatGF2(_mat_mul_jit(np.ascontiguousarray(self), np.ascontiguousarray(other)), copy=False)
 
     def compute_rank(self) -> np.intp:
         """Get the rank of the matrix.
@@ -149,7 +151,7 @@ def gauss_elimination(self, ncols: int | None = None, copy: bool = True) -> MatG
         ncols_value = self.shape[1] if ncols is None else ncols
         mat_ref = MatGF2(self) if copy else self
 
-        return MatGF2(_elimination_jit(mat_ref, ncols=ncols_value, full_reduce=False), copy=False)
+        return MatGF2(_elimination_jit(np.ascontiguousarray(mat_ref), ncols=ncols_value, full_reduce=False), copy=False)
 
     def row_reduction(self, ncols: int | None = None, copy: bool = True) -> MatGF2:
         """Return row-reduced echelon form (RREF) by performing Gaussian elimination.
@@ -170,7 +172,7 @@ def row_reduction(self, ncols: int | None = None, copy: bool = True) -> MatGF2:
         ncols_value = self.shape[1] if ncols is None else ncols
         mat_ref = self.copy() if copy else self
 
-        return MatGF2(_elimination_jit(mat_ref, ncols=ncols_value, full_reduce=True), copy=False)
+        return MatGF2(_elimination_jit(np.ascontiguousarray(mat_ref), ncols=ncols_value, full_reduce=True), copy=False)
 
 
 def solve_f2_linear_system(mat: MatGF2, b: MatGF2) -> MatGF2:
@@ -192,14 +194,17 @@ def solve_f2_linear_system(mat: MatGF2, b: MatGF2) -> MatGF2:
     -----
     This function is not integrated in `:class: graphix.linalg.MatGF2` because it does not perform any checks on the form of `mat` to ensure that it is in REF or that the system is solvable.
     """
-    return MatGF2(_solve_f2_linear_system_jit(mat, b), copy=False)
+    return MatGF2(_solve_f2_linear_system_jit(np.ascontiguousarray(mat), np.ascontiguousarray(b)), copy=False)
 
 
 @nb.njit("uint8[::1](uint8[:,::1], uint8[::1])")
 def _solve_f2_linear_system_jit(
     mat_data: npt.NDArray[np.uint8], b_data: npt.NDArray[np.uint8]
 ) -> npt.NDArray[np.uint8]:
-    """See docstring of `:func:solve_f2_linear_system` for details."""
+    """See docstring of `:func:solve_f2_linear_system` for details.
+
+    The signature of the numba decorator requires the input arrays to be C_CONTIGUOUS.
+    """
     m, n = mat_data.shape
     x = np.zeros(n, dtype=np.uint8)
 
@@ -238,6 +243,8 @@ def _solve_f2_linear_system_jit(
 def _elimination_jit(mat_data: npt.NDArray[np.uint8], ncols: int, full_reduce: bool) -> npt.NDArray[np.uint8]:
     r"""Return row echelon form (REF) or row-reduced echelon form (RREF) by performing Gaussian elimination.
 
+    The signature of the numba decorator requires the input arrays to be C_CONTIGUOUS.
+
     Parameters
     ----------
     mat_data : npt.NDArray[np.uint8]
@@ -302,7 +309,10 @@ def _elimination_jit(mat_data: npt.NDArray[np.uint8], ncols: int, full_reduce: b
 
 @nb.njit("uint8[:,::1](uint8[:,::1], uint8[:,::1])", parallel=True)
 def _mat_mul_jit(m1: npt.NDArray[np.uint8], m2: npt.NDArray[np.uint8]) -> npt.NDArray[np.uint8]:
-    """See docstring of `:func:MatGF2.__matmul__` for details."""
+    """See docstring of `:func:MatGF2.__matmul__` for details.
+
+    The signature of the numba decorator requires the input arrays to be C_CONTIGUOUS.
+    """
     m, l = m1.shape
     _, n = m2.shape
 

tests/test_linalg.py

@@ -8,6 +8,8 @@
 from graphix._linalg import MatGF2, solve_f2_linear_system
 
 if TYPE_CHECKING:
+    from typing import Literal
+
     from numpy.random import Generator
     from pytest_benchmark import BenchmarkFixture
 
@@ -82,6 +84,13 @@ def prepare_test_matrix() -> list[LinalgTestCase]:
             0,
             False,
         ),
+        # same as before but F-contiguous matrix
+        LinalgTestCase(
+            MatGF2(np.array([[1, 0], [0, 1], [1, 0]], dtype=np.uint8, order="F")),
+            2,
+            0,
+            False,
+        ),
     ]
 
 
@@ -110,6 +119,11 @@ def prepare_test_f2_linear_system() -> list[LSF2TestCase]:
                 mat=MatGF2([[1, 0, 1], [0, 1, 0], [0, 0, 1]]),
                 b=MatGF2([1, 1, 1]),
             ),
+            # Same as previous one but F-contiguous
+            LSF2TestCase(
+                mat=MatGF2(np.array([[1, 0, 1], [0, 1, 0], [0, 0, 1]], dtype=np.uint8, order="F")),
+                b=MatGF2([1, 1, 1]),
+            ),
         )
     )
 
@@ -210,11 +224,29 @@ def test_solve_f2_linear_system(self, benchmark: BenchmarkFixture, test_case: LS
 
         assert np.all((mat @ x) % 2 == b)  # Test with numpy matrix product.
 
-    def test_row_reduction(self, fx_rng: Generator) -> None:
-        sizes = [(10, 10), (3, 7), (6, 2)]
-        ncols = [4, 5, 2]
+    @pytest.mark.parametrize(
+        ("size", "ncol", "order"),
+        [
+            ((10, 10), 4, "K"),
+            ((3, 7), 5, "C"),
+            ((6, 2), 2, "F"),
+        ],
+    )
+    def test_row_reduction(
+        self, fx_rng: Generator, size: tuple[int, int], ncol: int, order: Literal["K", "C", "F"]
+    ) -> None:
+        mat = MatGF2(np.asarray(fx_rng.integers(size=size, low=0, high=2, dtype=np.uint8), order=order))
+        mat_red = mat.row_reduction(ncols=ncol, copy=True)
+        verify_elimination(mat, mat_red, ncol, full_reduce=True)
+
+    def test_initialization(self) -> None:
+        mat_c = MatGF2(np.array([[1, 0], [0, 1], [1, 0]], dtype=np.uint8))
+        mat_f = MatGF2(np.asfortranarray(np.array([[1, 0], [0, 1], [1, 0]], dtype=np.uint8)))
+
+        assert mat_c.flags.c_contiguous
+        assert not mat_c.flags.f_contiguous
+
+        assert not mat_f.flags.c_contiguous
+        assert mat_f.flags.f_contiguous
 
-        for size, ncol in zip(sizes, ncols, strict=True):
-            mat = MatGF2(fx_rng.integers(size=size, low=0, high=2, dtype=np.uint8))
-            mat_red = mat.row_reduction(ncols=ncol, copy=True)
-            verify_elimination(mat, mat_red, ncol, full_reduce=True)
+        assert np.all(mat_c == mat_f)

tests/test_opengraph.py

@@ -557,6 +557,33 @@ def _og_19() -> OpenGraphFlowTestCase:
     return OpenGraphFlowTestCase(og, has_cflow=False, has_gflow=False, has_pflow=True)
 
 
+@register_open_graph_flow_test_case
+def _og_20() -> OpenGraphFlowTestCase:
+    r"""Generate open graph.
+
+    Structure:
+
+    0
+    [(1)]-[(2)]
+
+    Notes
+    -----
+    This opengraph triggered issue #464.
+    https://github.com/TeamGraphix/graphix/issues/464
+    """
+    graph: nx.Graph[int] = nx.Graph([(1, 2)])
+    graph.add_node(0)
+    og = OpenGraph(
+        graph=graph,
+        input_nodes=[],
+        output_nodes=[1, 2],
+        measurements={
+            0: Measurement.YZ(angle=0),
+        },
+    )
+    return OpenGraphFlowTestCase(og, has_cflow=False, has_gflow=True, has_pflow=True)
+
+
 class OpenGraphComposeTestCase(NamedTuple):
     og1: OpenGraph[AbstractMeasurement]
     og2: OpenGraph[AbstractMeasurement]