diff --git a/CHANGELOG.md b/CHANGELOG.md index fc8da115..53772775 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/graphix/_linalg.py b/graphix/_linalg.py index 3667b597..09456069 100644 --- a/graphix/_linalg.py +++ b/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 diff --git a/tests/test_linalg.py b/tests/test_linalg.py index 2ae07363..a1ad1ce2 100644 --- a/tests/test_linalg.py +++ b/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) diff --git a/tests/test_opengraph.py b/tests/test_opengraph.py index 0cf48c1a..4f6347de 100644 --- a/tests/test_opengraph.py +++ b/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]