Reimplement flow matching using Gaussian elimination

.devcontainer/devcontainer.json

@@ -19,5 +19,5 @@
         }
     },
     "remoteUser": "vscode",
-    "postCreateCommand": "python -m pip install -r ./requirements.txt && python -m pip install -e ./[all] && pre-commit install && pre-commit run -a"
+    "postCreateCommand": "python -m pip install -e ./[all] && pre-commit install && pre-commit run -a"
 }

.github/workflows/ci.yml

@@ -16,7 +16,7 @@ jobs:
       - name: Set up Python
         uses: actions/setup-python@v5
         with:
-          python-version: "3.12"
+          python-version: "3.13"
           cache: "pip"
       - name: Install dependencies
         run: |
@@ -29,25 +29,23 @@ jobs:
     needs: pre-commit-ci
     strategy:
       matrix:
-        python-version: ["3.9", "3.10", "3.11", "3.12"]
+        python-version: ["3.10", "3.11", "3.12", "3.13", "3.14"]
     steps:
       - uses: actions/checkout@v4
       - name: Set up Python
         uses: actions/setup-python@v5
         with:
           python-version: ${{ matrix.python-version }}
           cache: "pip"
-      - name: Install dependencies
+      - name: Build
         run: |
           python -m pip install --upgrade pip
-          python -m pip install -r requirements.txt
+          python -m pip install -e .
       - name: Analyze code with pylint
         run: |
           python -m pip install pylint
           pylint $(git ls-files '*.py')
         continue-on-error: true
-      - name: Build
-        run: python -m pip install -e .
       - name: Test and coverage
         run: |
           python -m pip install pytest pytest-cov
@@ -69,7 +67,7 @@ jobs:
       - name: Set up Python
         uses: actions/setup-python@v5
         with:
-          python-version: "3.12"
+          python-version: "3.13"
           cache: "pip"
       - name: Install dependencies and tqec package
         run: |

.github/workflows/exotic-oses.yml

@@ -21,21 +21,19 @@ jobs:
         uses: actions/setup-python@v5
         with:
           # Use smallest supported Python version as the common denominator.
-          # In theory, if everything succeed in 3.9, everything should succeed in
-          # 3.10, 3.11 and follow-up version.
-          python-version: "3.9"
+          # In theory, if everything succeed in 3.10, everything should succeed in
+          # 3.11 and follow-up versions.
+          python-version: "3.10"
           cache: "pip"
-      - name: Install dependencies
+      - name: Build
         run: |
           python -m pip install --upgrade pip
-          python -m pip install -r requirements.txt
+          python -m pip install '.[all]'
       - name: Analyze code with pylint
         run: |
           python -m pip install pylint
           pylint $(git ls-files '*.py')
         continue-on-error: true
-      - name: Build
-        run: python -m pip install '.[all]'
       - name: Test and coverage
         run: |
           python -m pip install pytest pytest-cov

.github/workflows/gh-pages.yml

@@ -26,7 +26,6 @@ jobs:
           sudo apt-get update
           sudo apt-get install -y --no-install-recommends pandoc
           python -m pip install --upgrade pip
-          python -m pip install -r requirements.txt
           python -m pip install '.[all]'
       - name: Generate the documentation
         run: |

README.md

@@ -16,8 +16,6 @@ Currently, `tqecd` needs to be installed from source using
 python -m pip install git+https://github.com/tqec/tqecd.git
 ```
 
-The `tqecd` package has some dependencies that might be harder to install than a simple `pip install`. If you have any issues with the simple installation method above, please look at the [full installation page](https://tqec.github.io/tqecd/user_guide/installation.html).
-
 ## Basic usage
 
 ```py

docs/user_guide/installation.rst

@@ -10,21 +10,6 @@ Python version
 The ``tqecd`` package only supports Python 3.10 and onward. If you have Python 3.9 or below,
 please update your Python installation.
 
-Additional toolchains
-~~~~~~~~~~~~~~~~~~~~~
-
-Some of the dependencies of ``tqecd`` are implemented using compiled languages. This is for
-example the case of the `pycryptosat <https://pypi.org/project/pycryptosat/>`_ dependency.
-Pre-compiled Python packages that should be compatible with any GNU/Linux are provided
-by the author, but no pre-compiled package exist for Windows or MacOS.
-
-This means that, if you try to install the ``tqecd`` package on Windows or MacOS, a working
-C++ toolchain should also be installed on your system.
-
-Here is a list of potential issues you might encounter and how to solve them:
-
-- `Failed building wheel for pycryptosat <https://github.com/tqec/tqec/issues/311>`_
-
 Installation procedure
 ----------------------
 

pyproject.toml

@@ -4,15 +4,15 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "tqecd"
-version = "0.1.3"
+version = "0.1.4"
 authors = [
     { name = "TQEC community", email = "tqec-design-automation@googlegroups.com" },
 ]
 description = "Automatically find detectors in a topologically quantum error corrected computation"
 readme = "README.md"
 license = { file = "LICENSE" }
 keywords = ["topological quantum error correction", "qec", "detectors"]
-requires-python = ">= 3.9"
+requires-python = ">= 3.10"
 classifiers = [
     "Development Status :: 4 - Beta",
     "Intended Audience :: Developers",
@@ -28,12 +28,19 @@ classifiers = [
     "Programming Language :: Python :: 3.11",
     "Programming Language :: Python :: 3.12",
     "Programming Language :: Python :: 3.13",
+    "Programming Language :: Python :: 3.14",
     "Topic :: Scientific/Engineering",
     "Topic :: Software Development :: Libraries",
     "Topic :: Utilities",
     "Typing :: Typed",
 ]
-dynamic = ["dependencies"]
+dependencies = [
+    "numpy>=1.22,<2.3",                      # Upper bound: https://github.com/tqec/tqec/pull/659
+    "numpy>=1.24; python_version >= '3.11'",
+    "numpy>=1.26; python_version >= '3.12'",
+    "numpy>=2.1; python_version >= '3.13'",
+    "stim>=1.15.0",
+]
 
 [project.urls]
 Documentation = "https://tqec.github.io/tqecd/"
@@ -64,10 +71,6 @@ exclude = ["*test.py"]
 [tool.setuptools.package-data]
 tqecd = ["py.typed"]
 
-[tool.setuptools.dynamic]
-dependencies = { file = "requirements.txt" }
-
-
 [tool.mypy]
 # See https://numpy.org/devdocs/reference/typing.html
 plugins = "numpy.typing.mypy_plugin"
@@ -106,7 +109,7 @@ warn_return_any = true
 
 # See https://mypy.readthedocs.io/en/stable/config_file.html#using-a-pyproject-toml
 [[tool.mypy.overrides]]
-module = ["stim", "pysat", "pysat.solvers"]
+module = ["stim"]
 ignore_missing_imports = true
 
 [tool.coverage.run]

src/tqecd/cover.py

@@ -0,0 +1,150 @@
+"""Provides utility functions to find "covers" of Pauli strings."""
+
+from __future__ import annotations
+
+from collections.abc import Callable, Iterable
+from typing import Any
+
+from tqecd.pauli import PauliString
+
+
+def _find_cover(
+    target: PauliString, sources: list[PauliString], on_qubits: frozenset[int]
+) -> list[int] | None:
+    """Try to find a set of boundary stabilizers from ``sources`` that generate
+    target on qubits ``on_qubits``.
+
+    If multiple valid covers exist, the covers involving the lowest number of
+    :class:`~tqecd.pauli.PauliString` instances from ``sources`` are listed, and
+    a random cover is picked from that list.
+
+    Args:
+        target: the stabilizers to cover with stabilizers from ``sources``.
+        sources: stabilizers that can be used to cover `target`.
+        on_qubits: qubits to consider when trying to cover ``target`` with
+            ``sources``.
+
+    Returns:
+        Either a list of indices over ``sources`` that, when combined, cover
+        exactly the provided ``target`` on all the qubits provided in
+        ``on_qubits``, or ``None`` if such a list could not be found.
+    """
+    return _solve_linear_system(
+        _construct_basis({}, sources, lambda s: s.to_int(on_qubits)),
+        target.to_int(on_qubits),
+    )
+
+
+def find_exact_cover(
+    target: PauliString, sources: list[PauliString]
+) -> list[int] | None:
+    """Try to find a set of pauli strings from ``sources`` that generate exactly
+    ``target``.
+
+    The Pauli strings returned (via indices over the provided ``sources``), once
+    multiplied together, should be exactly equal to ``target``. In particular, the
+    following post-condition should hold:
+
+    .. code-block:: python
+
+        target = None     # to replace
+        sources = [None]  # to replace
+        cover_indices = find_exact_cover(target, sources)
+        resulting_pauli_string = PauliString({})
+        for i in cover_indices:
+            resulting_pauli_string = resulting_pauli_string * sources[i]
+        assert resulting_pauli_string == target, "Should hold"
+
+    Args:
+        target: the stabilizers to cover with stabilizers from ``sources``.
+        sources: stabilizers that can be used to cover ``target``.
+
+    Returns:
+        Either a list of indices over ``sources`` that, when combined, cover
+        exactly the provided ``target``, or ``None`` if such a list could not be
+        found.
+    """
+    # If target is the identity, pick no Pauli string.
+    # Note: we might want to disallow an empty return in the future.
+    if target.non_trivial_pauli_count == 0:
+        return []
+
+    # Else, if there are no sources, we cannot find a solution.
+    if not sources:
+        return None
+
+    # We want an exact (i.e., equality) cover on all qubits, to be sure that
+    # the post-condition in the docstring holds. For that, it is sufficient to
+    # only consider all the qubits where either `target` or at least one of the
+    # items of `sources` acts non-trivially (i.e., something else than the identity).
+    involved_qubits = frozenset(target.qubits)
+    for source in sources:
+        involved_qubits |= frozenset(source.qubits)
+
+    return _find_cover(target, sources, involved_qubits)
+
+
+def find_commuting_cover_on_target_qubits(
+    target: PauliString, sources: list[PauliString]
+) -> list[int] | None:
+    """Try to find a set of boundary stabilizers from ``sources`` that generate a
+    superset of ``target``.
+
+    This function try to find a set of Pauli strings from ``sources`` that
+    includes ``target`` (i.e., on every qubit where ``target`` is non-trivial,
+    the product of each of the returned Pauli strings should commute with
+    ``target``).
+
+    The differences with :func:`find_cover` are:
+
+    1. this function does not restrict the output of the product of each of
+       the returned Pauli string on qubits where ``target`` acts trivially (i.e.
+       "I"). So in practice, on qubits where ``target[qubit] == "I"``, the value
+       of the returned Pauli string can be anything.
+    2. this function does not restrict the output of the product of each of
+       the returned Pauli string to exactly match ``target`` on qubits where
+       it acts non-trivially, but rather requires the output to commute with
+       ``target`` on those qubits.
+
+    Args:
+        target: the stabilizers to cover with stabilizers from ``sources``.
+        sources: stabilizers that can be used to cover ``target``.
+
+    Returns:
+        Either a list of a stabilizers that, when combined, commute with
+        the provided ``target``, or ``None`` if such a list could not be found.
+    """
+    if not sources:
+        return None
+    return _find_cover(target, sources, frozenset(target.qubits))
+
+
+def _solve_linear_system(
+    basis: dict[int, tuple[int, int]], x: int, update_basis: bool = True
+) -> list[int] | None:
+    """Gaussian elimination over GF(2)."""
+    mask = 1 << len(basis)
+    while x:
+        highest_bit = x.bit_length() - 1
+        if highest_bit not in basis:
+            if update_basis:
+                basis[highest_bit] = (x, mask)
+            return None
+        pivot, pivot_mask = basis[highest_bit]
+        x ^= pivot
+        mask ^= pivot_mask
+    return _int_to_bit_indices(mask)[:-1]
+
+
+def _int_to_bit_indices(x: int) -> list[int]:
+    """Convert an integer to a list of indices where the bits are set."""
+    return [i for i in range(x.bit_length()) if (x >> i) & 1]
+
+
+def _construct_basis(
+    basis: dict[int, tuple[int, int]], items: Iterable[Any], func: Callable[[Any], int]
+) -> dict[int, tuple[int, int]]:
+    """Construct a linear basis from the given items using the provided function."""
+    for item in items:
+        _solve_linear_system(basis, func(item))
+    return basis