Add native Clifford map compilation pass

graphix/circ_ext/compilation.py

@@ -5,14 +5,23 @@
 from itertools import chain, pairwise
 from typing import TYPE_CHECKING
 
+import numpy as np
+
 from graphix.fundamentals import ANGLE_PI, Axis
+from graphix.instruction import CNOT, SWAP, H, S, X, Y, Z
 from graphix.sim.base_backend import NodeIndex
 from graphix.transpiler import Circuit
 
 if TYPE_CHECKING:
     from collections.abc import Callable
+    from typing import TypeAlias
 
+    from graphix._linalg import MatGF2
     from graphix.circ_ext.extraction import CliffordMap, ExtractionResult, PauliExponential, PauliExponentialDAG
+    from graphix.instruction import Instruction
+
+    # NOTE: This alias could be defined at the level of graphix.instruction, and treat all qubit indices as `Qubit`. This change would affect many files in the codebase, so as a temporary solution `Qubit` is casted to `int` in this module.
+    Qubit: TypeAlias = int | np.int_
 
 
 def er_to_circuit(
@@ -31,7 +40,7 @@ def er_to_circuit(
     pexp_cp: Callable[[PauliExponentialDAG, Circuit], None] | None
         Compilation pass to synthesize a Pauli exponential DAG. If ``None`` (default), :func:`pexp_ladder_pass` is employed.
     cm_cp: Callable[[CliffordMap, Circuit], None] | None
-        Compilation pass to synthesize a Clifford map. If ``None`` (default), a ``ValueError`` is raised since there is still no default pass for Clifford map integrated in Graphix.
+        Compilation pass to synthesize a Clifford map. If ``None`` (default), :func:`cm_berg_pass` is employed. This pass only handles unitaries so far (Clifford maps with the same number of input and ouptut nodes).
 
     Returns
     -------
@@ -51,9 +60,7 @@ def er_to_circuit(
         pexp_cp = pexp_ladder_pass
 
     if cm_cp is None:
-        raise ValueError(
-            "Clifford-map pass is missing: there is still no default pass for Clifford map integrated in Graphix. You may use graphix-stim-compiler plugin."
-        )
+        cm_cp = cm_berg_pass
 
     n_qubits = len(er.pexp_dag.output_nodes)
     circuit = Circuit(n_qubits)
@@ -181,3 +188,167 @@ def add_hy(qubit: int, circuit: Circuit) -> None:
     for node in chain(*reversed(pexp_dag.partial_order_layers[1:])):
         pexp = pexp_dag.pauli_exponentials[node]
         add_pexp(pexp, circuit)
+
+
+def cm_berg_pass(clifford_map: CliffordMap, circuit: Circuit) -> None:
+    r"""Add a Clifford map to a circuit by using an adaptation of van den Berg's sweeping algorithm introduced in Ref. [1].
+
+    The input circuit is modified in-place. This function assumes that the Clifford Map has been remap, i.e., its Pauli strings are defined on qubit indices instead of output nodes. See :meth:`PauliString.remap` for additional information.
+
+    Parameters
+    ----------
+    clifford_map: CliffordMap
+        The Clifford map to be transpiled. Its Pauli strings are assumed to be defined on qubit indices.
+    circuit : Circuit
+        The circuit to which the operation is added. The input circuit is assumed to be compatible with
+        ``CliffordMap.input_nodes`` and ``CliffordMap.output_nodes``.
+
+    Raises
+    ------
+    NotImplementedError
+        If ``len(clifford_map.input_nodes) != len(clifford_map.output_nodes)``.
+    AssertionError
+        If an unexpected pivot position is encountered during Step 4.
+
+    Notes
+    -----
+    This pass only handles unitaries so far (Clifford maps with the same number of input and ouptut nodes).
+
+    Gate set: H, S, CNOT, SWAP, X, Y, Z
+
+    This function converts a ``CliffordMap`` into a sequence of quantum
+    gate instructions by operating on its binary tableau representation.
+    The synthesis proceeds qubit-by-qubit, applying a sequence of local
+    transformations (H, S, CNOT, and SWAP gates) to reduce the
+    tableau into a canonical form. The resulting sequence represents the
+    adjoint of the input Clifford map, therefore it's appended in reverse
+    order (and exchanging S by Sdagger) to the provided ``Circuit``.
+
+    The synthesis applies a series of steps on every qubit subtableau:
+
+    .. math::
+
+    T_q = \begin{pmatrix}
+        XX & XZ \\
+        ZX & ZZ
+    \end{pmatrix}
+
+    1. Clear elements in the XZ-block by applying single-qubit gates (H or S).
+
+    2. Use CNOT gates to reduce the XX block to a single pivot column.
+
+    3. Apply a SWAP gate to bring the pivot to the diagonal if neccesary.
+
+    4. Ensure the ZX and ZZ blocks of the tableau have the correct canonical form
+    by redoing steps 1. and 2.
+
+    After processing all qubits, a final sign correction step applies
+    Pauli gates (X, Y, Z) to fix phase bits in the tableau.
+
+    The generated instructions are accumulated during the forward pass
+    and then appended to the circuit in reverse order to yield the
+    correct overall transformation.
+
+    For the mapping between tableau updates and Clifford gates (H, S, CNOT) see [2].
+
+    References
+    ----------
+    [1] Van Den Berg, 2021. A simple method for sampling random Clifford operators (arxiv:2008.06011).
+    [2] Aaronson, Gottesman, (2004). Improved Simulation of Stabilizer Circuits (arXiv:quant-ph/0406196).
+    """
+    tab = clifford_map.to_tableau()
+    n = len(clifford_map.output_nodes)
+    if len(clifford_map.input_nodes) != n:
+        raise NotImplementedError(
+            ":func:`cm_berg_pass` does not support circuit compilation if the number of input and output nodes is different (isometry)."
+        )
+    instructions: list[Instruction] = []
+
+    def process_qubit(tab: MatGF2, instructions: list[Instruction], q: int) -> None:
+        """Bring to canonical form two tableau rows corresponding to qubit ``q``."""
+        # Step 1
+        do_step_1(tab, instructions, row_idx=q)
+
+        # Step 2
+        pivot = do_step_2(tab, instructions, row_idx=q)
+
+        # Step 3
+        if pivot != q:
+            add_swap(tab, instructions, q, pivot)
+
+        # Step 4
+        col_idx_z = np.flatnonzero(tab[q + n, :-1])  # ZX and ZZ blocks of qubit q, without sign.
+        if not (len(col_idx_z) == 1 and col_idx_z[0] == q + n):
+            # ZX and ZZ blocks don't have the canonical form.
+            add_h(tab, instructions, q)
+            do_step_1(tab, instructions, row_idx=q + n)
+            pivot = do_step_2(tab, instructions, row_idx=q + n)
+            if pivot != q:
+                raise AssertionError(
+                    f"Pivot in block ZZ should be at q = {q}. This error probably means that `CliffordMap` doesn't describe a valid Clifford operation. All Pauli strings must commute, except for `x_map[q]` anticommuting with `z_map[q]` for each q."
+                )
+            add_h(tab, instructions, q)
+
+    def do_step_1(tab: MatGF2, instructions: list[Instruction], row_idx: int) -> None:
+        col_idx_zx = np.flatnonzero(tab[row_idx, n : 2 * n])  # Don't take the sign column
+        for j in col_idx_zx:
+            # Each iteration sets the element `tab[row_idx, n+j]` to 0.
+            add_s(tab, instructions, int(j)) if tab[row_idx, j] else add_h(tab, instructions, int(j))
+
+    def do_step_2(tab: MatGF2, instructions: list[Instruction], row_idx: int) -> int:
+        col_idx_xx = np.flatnonzero(tab[row_idx, :n])
+        while len(col_idx_xx) > 1:
+            for i in range(0, len(col_idx_xx) - 1, 2):  # itertools.batched only avaialble in Python 3.12+
+                # Apply CNOTS to disjoint qubits in parallel
+                add_cnot(tab, instructions, qc=col_idx_xx[i], qt=col_idx_xx[i + 1])
+            col_idx_xx = col_idx_xx[::2]
+
+        return int(col_idx_xx[0])  # Return pivot
+
+    def add_h(tab: MatGF2, instructions: list[Instruction], q: Qubit) -> None:
+        q = int(q)  # Cast to `int` to avoid typing issues
+        tab[:, -1] ^= tab[:, q] & tab[:, q + n]
+        tab[:, [q, q + n]] = tab[:, [q + n, q]]  # The usual tuple assignment `a, b = b, a` does not work here.
+        instructions.append(H(q))
+
+    def add_s(tab: MatGF2, instructions: list[Instruction], q: Qubit) -> None:
+        tab[:, -1] ^= tab[:, q] & tab[:, q + n]
+        tab[:, q + n] = tab[:, q] ^ tab[:, q + n]
+        q = int(q)
+        instructions.extend((S(q), Z(q)))  # We append Sdagger to get C instead of C^dagger
+
+    def add_cnot(tab: MatGF2, instructions: list[Instruction], qc: Qubit, qt: Qubit) -> None:
+        tab[:, -1] ^= tab[:, qc] * tab[:, qt + n] * (tab[:, qt] ^ tab[:, qc + n] ^ 1)
+        tab[:, qt] ^= tab[:, qc]
+        tab[:, qc + n] ^= tab[:, qt + n]
+        instructions.append(CNOT(control=int(qc), target=int(qt)))
+
+    def add_swap(tab: MatGF2, instructions: list[Instruction], q0: Qubit, q1: Qubit) -> None:
+        q0, q1 = int(q0), int(q1)  # Cast to `int` to avoid typing issues
+        for shift in [0, n]:
+            tab[:, [q0 + shift, q1 + shift]] = tab[:, [q1 + shift, q0 + shift]]
+
+        instructions.append(SWAP((q0, q1)))
+
+    def correct_signs(tab: MatGF2, instructions: list[Instruction]) -> None:
+        for q in range(n):
+            sign_xz = tab[q, -1], tab[q + n, -1]
+            match sign_xz:
+                case (0, 1):
+                    instructions.append(X(q))
+                case (1, 1):
+                    instructions.append(Y(q))
+                case (1, 0):
+                    instructions.append(Z(q))
+
+            # The tableau sign column should be set to 0, but we don't need to do it since it's the last step.
+            # tab[q, -1], tab[q + n, -1] = 0, 0
+
+    for q in range(n):
+        process_qubit(tab, instructions, q)
+
+    correct_signs(tab, instructions)
+
+    # Append instructions in reverse order to get C instead of Cdagger
+    for instr in instructions[::-1]:
+        circuit.add(instr)

graphix/circ_ext/extraction.py

@@ -5,8 +5,10 @@
 from dataclasses import dataclass, replace
 from typing import TYPE_CHECKING
 
+import numpy as np
 from typing_extensions import Self  # Self introduced in 3.11
 
+from graphix._linalg import MatGF2
 from graphix.fundamentals import Axis, ParameterizedAngle, Plane, Sign
 from graphix.measurements import BlochMeasurement, Measurement, PauliMeasurement
 
@@ -56,9 +58,9 @@ def to_circuit(
         Parameters
         ----------
         pexp_cp: Callable[[PauliExponentialDAG, Circuit], None] | None
-            Compilation pass to synthesize a Pauli exponential DAG. If ``None`` (default), :func:`pexp_ladder_pass` is employed.
-        cm_cp: Callable[[PauliExponentialDAG, Circuit], None] | None
-            Compilation pass to synthesize a Clifford map. If ``None`` (default), a `ValueError` is raised since there is still no default pass for Clifford map integrated in Graphix.
+            Compilation pass to synthesize a Pauli exponential DAG. If ``None`` (default), :func:`graphix.circ_ext.compilation.pexp_ladder_pass` is employed.
+        cm_cp: Callable[[CliffordMap, Circuit], None] | None
+            Compilation pass to synthesize a Clifford map. If ``None`` (default), :func:`graphix.circ_ext.compilation.cm_berg_pass` is employed. This pass only handles unitaries so far (Clifford maps with the same number of input and ouptut nodes).
 
         Returns
         -------
@@ -369,6 +371,63 @@ def remap(self, inputs_mapping: Callable[[int], int], outputs_mapping: Callable[
         z_map = {inputs_mapping(node): ps.remap(outputs_mapping) for node, ps in self.z_map.items()}
         return replace(self, x_map=x_map, z_map=z_map)
 
+    def to_tableau(self) -> MatGF2:
+        """Convert the CliffordMap into its binary tableau representation.
+
+        The returned tableau is a ``(2n, 2n + 1)`` binary matrix over GF(2),
+        where ``n`` is the number of qubits. The first ``n`` rows correspond
+        to the images of X generators, and the next ``n`` rows correspond to
+        the images of Z generators. Columns encode the X and Z components of
+        the resulting Pauli strings, along with a sign column.
+
+        Each PauliString in ``x_map`` and ``z_map`` is decomposed into its
+        X/Z support:
+        - X contributes a 1 to the X block.
+        - Z contributes a 1 to the Z block.
+        - Y contributes a 1 to both X and Z blocks.
+        The sign of the Pauli string is stored in the final column
+        (0 for ``Sign.PLUS`` and 1 for ``Sign.MINUS``).
+
+        Parameters
+        ----------
+        None
+            This method operates on the current CliffordMap instance. It
+            assumes the map has already been remapped such that input and
+            output nodes correspond to qubit indices.
+
+        Returns
+        -------
+        MatGF2
+            A binary matrix of shape ``(2n, 2n + 1)`` representing the
+            Clifford tableau.
+
+        Raises
+        ------
+        NotImplementedError
+            If the number of input nodes differs from the number of output
+            nodes (i.e., the map is an isometry instead of a square Clifford).
+        """
+        n = len(self.input_nodes)
+        if n != len(self.output_nodes):
+            raise NotImplementedError(
+                f"Isometries are not supported yet: # of inputs ({len(self.input_nodes)}) must be equal to the # of outputs ({len(self.output_nodes)})."
+            )
+
+        tab = MatGF2(np.zeros((2 * n, 2 * n + 1)))
+
+        for mapping, shift in zip((self.x_map, self.z_map), (0, n), strict=True):
+            for i, ps in mapping.items():  # Clifford map has been remap so keys correspond to qubits.
+                for j, ax in ps.axes.items():
+                    if ax in {Axis.X, Axis.Y}:
+                        tab[i + shift, j] = 1
+                    if ax in {Axis.Y, Axis.Z}:
+                        tab[i + shift, j + n] = 1
+
+                if ps.sign is Sign.MINUS:
+                    tab[i + shift, 2 * n] = 1
+
+        return tab
+
 
 def extend_input(og: OpenGraph[Measurement]) -> tuple[OpenGraph[Measurement], dict[int, int]]:
     r"""Extend the inputs of a given open graph.

graphix/opengraph.py

@@ -17,9 +17,11 @@
 if TYPE_CHECKING:
     from collections.abc import Callable, Collection, Iterable, Mapping, Sequence
 
+    from graphix.circ_ext.extraction import CliffordMap, PauliExponentialDAG
     from graphix.flow.core import CausalFlow
     from graphix.parameter import ExpressionOrSupportsFloat, Parameter
     from graphix.pattern import Pattern
+    from graphix.transpiler import Circuit
 
 # TODO: Maybe move these definitions to graphix.fundamentals and graphix.measurements ? Now they are redefined in graphix.flow._find_gpflow, not very elegant.
 _AM_co = TypeVar("_AM_co", bound=AbstractMeasurement, covariant=True)
@@ -541,6 +543,87 @@ def find_pauli_flow(self: OpenGraph[_AM_co], *, stacklevel: int = 1) -> PauliFlo
             correction_matrix
         )  # The constructor returns `None` if the correction matrix is not compatible with any partial order on the open graph.
 
+    def extract_circuit(
+        self: OpenGraph[Measurement],
+        pexp_cp: Callable[[PauliExponentialDAG, Circuit], None] | None = None,
+        cm_cp: Callable[[CliffordMap, Circuit], None] | None = None,
+        *,
+        stacklevel: int = 1,
+    ) -> Circuit:
+        """Extract a unitary in the form of a circuit from an open graph resource state.
+
+        This method acts as a wrapper around the circuit extraction routine, simplifying
+        its usage. It first attempts to extract the Pauli flow of the open graph, then
+        applies the circuit extraction procedure described in Ref. [1], and finally compiles
+        the resulting circuit using the provided passes.
+        To obtain the open graph's unitary in the form of a Pauli exponential DAG along with
+        a Clifford transformation, as presented in Ref. [1], one should instead operate
+        directly on the flow object using :meth:`PauliFlow.extract_circuit`.
+
+        Parameters
+        ----------
+        pexp_cp: Callable[[PauliExponentialDAG, Circuit], None] | None
+            Compilation pass to synthesize a Pauli exponential DAG.
+            If ``None`` (default), :func:`graphix.circ_ext.compilation.pexp_ladder_pass` is
+            employed.
+        cm_cp: Callable[[CliffordMap, Circuit], None] | None
+            Compilation pass to synthesize a Clifford map. If ``None`` (default),
+            :func:`graphix.circ_ext.compilation.cm_berg_pass` is employed. This pass
+            only handles unitaries so far (Clifford maps with the same number of input
+            and ouptut nodes).
+        stacklevel : int, optional
+            Stack level to use for warnings. Defaults to 1, meaning that warnings
+            are reported at this function's call site.
+
+        Returns
+        -------
+        Circuit
+            Quantum circuit represented as a set of instructions.
+
+        Notes
+        -----
+        - The open graph instance must be of parametric type ``Measurement`` to allow
+        for a circuit extraction, otherwise it does not contain information about the
+        measurement angles.
+
+        - This wrapper extracts a Pauli flow rather than a gflow, as the former is more
+        general while the underlying extraction algorithms have the same computational
+        complexity in both cases. The resulting unitary is identical whether it is
+        obtained from a ``GFlow`` or from a ``PauliFlow`` with inferred Pauli measurements.
+        However, compilation passes that simultaneously diagonalize Pauli exponentials
+        within the same layer of the Pauli exponential DAG may benefit from flows of
+        lower depth, which is often the case for Pauli flow.  The pass
+        :func:`graphix.circ_ext.compilation.pexp_ladder_pass` does not take into account
+        the flow's depth.
+
+        References
+        ----------
+        [1] Simmons, 2021 (arXiv:2109.05654).
+
+        Examples
+        --------
+        >>> import networkx as nx
+        >>> from graphix.opengraph import OpenGraph
+        >>> from graphix.measurements import Measurement
+        >>> og = OpenGraph(
+        ...     graph=nx.Graph([(0, 1), (1, 2), (3, 4), (4, 5), (6, 7), (7, 8), (1, 3), (4, 6)]),
+        ...     input_nodes=(0, 3, 6),
+        ...     output_nodes=(2, 5, 8),
+        ...     measurements=dict.fromkeys((0, 1, 3, 4, 6, 7), Measurement.XY(angle=0)),
+        ... )
+        >>> og.extract_circuit()
+        Circuit(width=3, instr=[H(2), H(1), CNOT(2, 1), H(1), H(1), H(0), CNOT(2, 0), CNOT(1, 0), H(2), H(1), H(0)])
+        >>> # The default compilation passes do not exploit the lower depth of the Pauli flow
+        >>> # compared the gflow.
+        >>> og.infer_pauli_measurements().extract_circuit()
+        Circuit(width=3, instr=[H(2), H(1), CNOT(2, 1), H(1), H(1), H(0), CNOT(2, 0), CNOT(1, 0), H(2), H(1), H(0)])
+        """
+        return (
+            self.extract_pauli_flow(stacklevel=stacklevel + 1)
+            .extract_circuit()
+            .to_circuit(pexp_cp=pexp_cp, cm_cp=cm_cp)
+        )
+
     def compose(self, other: OpenGraph[_AM_co], mapping: Mapping[int, int]) -> tuple[OpenGraph[_AM_co], dict[int, int]]:
         r"""Compose two open graphs by merging subsets of nodes from ``self`` and ``other``, and relabeling the nodes of ``other`` that were not merged.