diff --git a/CHANGELOG.md b/CHANGELOG.md
index 897d20fed..b6e83814e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,10 +9,22 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 
+- #476 Introduced new methods `OpenGraph.extract_circuit`, `CliffordMap.to_tableau` and new function `graphix.circ_ext.compilation.cm_berg_pass`. Circuit extraction can be done natively in Graphix.
+
+### Fixed
+
+### Changed
+
 - #479: Added new methods `OpenGraph.draw`, `PauliFlow.draw` and `XZCorrections.draw`.
 
 - #454, #481: New space minimization API that allows users to select or define custom heuristics.
 
+- #476
+  - Added new field `dim` to `PauliString` to represent the dimension of the Hilbert space.
+  - Define `PauliString` on qubit indices and remove all `remap` methods.
+  - Represent `x_map` and `z_map` attributes of `CliffordMap` as sequences of `PauliString` instead of mappings.
+  - Rename `PauliFlow.pauli_strings` property as `PauliFlow.extraction_pauli_strings`.
+
 ### Fixed
 
 - #454, #481: Ensure `Pattern.minimize_space` only reduces max-space and does not increase it.
diff --git a/README.md b/README.md
index 39603616d..d40169b05 100644
--- a/README.md
+++ b/README.md
@@ -90,13 +90,17 @@ state_out = pattern.simulate_pattern(backend="statevector")
 - For theoretical background, read our quick introduction into [MBQC](https://graphix.readthedocs.io/en/latest/intro.html) and [LC-MBQC](https://graphix.readthedocs.io/en/latest/lc-mbqc.html).
 - Full API docs is [here](https://graphix.readthedocs.io/en/latest/references.html).
 
-## Related packages
+## Graphix plugins
 
 - [graphix-stim-backend](https://github.com/thierry-martinez/graphix-stim-backend): `stim` backend for efficient Clifford pattern simulation
 - [graphix-symbolic](https://github.com/TeamGraphix/graphix-symbolic): parameterized patterns with symbolic simulation
 - [graphix-ibmq](https://github.com/TeamGraphix/graphix-ibmq): pattern transpiler for IBMQ / `qiskit`
 - [graphix-perceval](https://github.com/TeamGraphix/graphix-perceval): pattern transpiler for Quandela's `perceval` simulator and QPU
 - [graphix-qasm-parser](https://github.com/TeamGraphix/graphix-qasm-parser): a plugin for parsing OpenQASM circuit.
+- [`graphix-stim-compiler`](https://github.com/qat-inria/graphix-stim-compiler): `stim` backend for efficient compilation of Clifford maps.
+
+## Related packages
+
 - [swiflow](https://github.com/TeamGraphix/swiflow): rust-based implementation of flow-finding algorithms.
 - [graphqomb](https://github.com/TeamGraphix/graphqomb): modular graph state compiler for fault-tolerant MBQC and more.
 
diff --git a/graphix/circ_ext/compilation.py b/graphix/circ_ext/compilation.py
index 503532910..e8897c9ae 100644
--- a/graphix/circ_ext/compilation.py
+++ b/graphix/circ_ext/compilation.py
@@ -5,14 +5,18 @@
 from itertools import chain, pairwise
 from typing import TYPE_CHECKING
 
+import numpy as np
+
 from graphix.fundamentals import ANGLE_PI, Axis
-from graphix.sim.base_backend import NodeIndex
+from graphix.instruction import CNOT, SWAP, H, S, X, Y, Z
 from graphix.transpiler import Circuit
 
 if TYPE_CHECKING:
     from collections.abc import Callable
 
+    from graphix._linalg import MatGF2
     from graphix.circ_ext.extraction import CliffordMap, ExtractionResult, PauliExponential, PauliExponentialDAG
+    from graphix.instruction import Instruction
 
 
 def er_to_circuit(
@@ -22,7 +26,7 @@ def er_to_circuit(
 ) -> Circuit:
     """Convert a circuit extraction result into a quantum circuit representation.
 
-    This method synthesizes a circuit by sequentially applying the Clifford map and the Pauli exponential DAG (Directed Acyclic Graph) in the extraction result. It performs a validation check to ensure that the output nodes of both components are identical and it maps the output node numbers to qubit indices.
+    This method synthesizes a circuit by sequentially applying the Clifford map and the Pauli exponential DAG (Directed Acyclic Graph) in the extraction result. It performs a validation check to ensure that the output nodes of both components are identical.
 
     Parameters
     ----------
@@ -31,7 +35,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 output nodes).
 
     Returns
     -------
@@ -51,32 +55,25 @@ 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)
-    outputs_mapping = NodeIndex()
-    outputs_mapping.extend(er.pexp_dag.output_nodes)
 
-    inputs_mapping = NodeIndex()
-    inputs_mapping.extend(er.clifford_map.input_nodes)
-
-    cm_cp(er.clifford_map.remap(inputs_mapping.index, outputs_mapping.index), circuit)
-    pexp_cp(er.pexp_dag.remap(outputs_mapping.index), circuit)
+    cm_cp(er.clifford_map, circuit)
+    pexp_cp(er.pexp_dag, circuit)
     return circuit
 
 
 def pexp_ladder_pass(pexp_dag: PauliExponentialDAG, circuit: Circuit) -> None:
     r"""Add a Pauli exponential DAG to a circuit by using a ladder decomposition.
 
-    The input circuit is modified in-place. This function assumes that the Pauli exponential DAG has been remapped, i.e., its Pauli strings are defined on qubit indices instead of output nodes. See :meth:`PauliString.remap` for additional information.
+    The input circuit is modified in-place.
 
     Parameters
     ----------
     pexp_dag: PauliExponentialDAG
-        The Pauli exponential rotation to be added to the circuit. Its Pauli strings are assumed to be defined on qubit indices.
+        The Pauli exponential rotation to be added to the circuit. Its Pauli strings are defined on qubit indices.
     circuit : Circuit
         The circuit to which the operation is added. The input circuit is assumed to be compatible with ``pexp_dag.output_nodes``.
 
@@ -106,15 +103,11 @@ def add_pexp(pexp: PauliExponential, circuit: Circuit) -> None:
             The Pauli exponential to add.
         circuit : Circuit
             The quantum circuit to which the Pauli exponential is added.
-
-        Notes
-        -----
-        It is assumed that the ``x``, ``y``, and ``z`` node sets of the Pauli string in the exponential are well-formed, i.e., contain valid qubit indices and are pairwise disjoint.
         """
         if pexp.angle == 0:  # No rotation
             return
 
-        # We assume that nodes in the Pauli strings have been mapped to qubits.
+        # Pauli strings are defined on qubit indices.
         # The order on which we iterate over the modified qubits does not matter.
         modified_qubits = list(pexp.pauli_string.axes)
         angle = -2 * pexp.angle * pexp.pauli_string.sign
@@ -181,3 +174,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.
+
+    Parameters
+    ----------
+    clifford_map: CliffordMap
+        The Clifford map to be transpiled.
+    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 output 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 necessary.
+
+    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 available 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: int) -> 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: int) -> None:
+        tab[:, -1] ^= tab[:, q] & tab[:, q + n]
+        tab[:, q + n] ^= tab[:, q]
+        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: int, qt: int) -> 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: int, q1: int) -> None:
+        q0, q1 = int(q0), int(q1)  # Cast to `int` to avoid typing issues
+        tab[:, [q0, q1, q0 + n, q1 + n]] = tab[:, [q1, q0, q1 + n, q0 + n]]
+
+        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)
diff --git a/graphix/circ_ext/extraction.py b/graphix/circ_ext/extraction.py
index 88b521295..fd33537ec 100644
--- a/graphix/circ_ext/extraction.py
+++ b/graphix/circ_ext/extraction.py
@@ -5,10 +5,12 @@
 from dataclasses import dataclass, replace
 from typing import TYPE_CHECKING
 
-from typing_extensions import Self  # Self introduced in 3.11
+import numpy as np
 
+from graphix._linalg import MatGF2
 from graphix.fundamentals import Axis, ParameterizedAngle, Plane, Sign
 from graphix.measurements import BlochMeasurement, Measurement, PauliMeasurement
+from graphix.sim.base_backend import NodeIndex
 
 if TYPE_CHECKING:
     from collections.abc import Callable, Mapping, Sequence
@@ -55,10 +57,10 @@ 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.
+        pexp_cp: Callable[[PauliExponentialDAG, Circuit], None] | None, default
+            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 output nodes).
 
         Returns
         -------
@@ -73,20 +75,23 @@ def to_circuit(
 
 @dataclass(frozen=True)
 class PauliString:
-    """Dataclass representing a Pauli string over a set of MBQC nodes.
+    """Dataclass representing a Pauli string.
 
     Attributes
     ----------
+    dim : int
+        Dimension of the Hilbert space on which the Pauli string acts.
     axes : Mapping[int, Axis]
-        Mapping between nodes and the applied Pauli operator.
+        Mapping between qubit indices and the applied Pauli operator.
     sign : Sign
-        Phase of the Pauli string.
+        Sign of the Pauli string.
 
     Notes
     -----
-    The identity operator is omitted in this representation, which means that in general it is not possible to infer the size of the Hilbert space from an instance of ``PauliString`` alone.
+    The identity operators in the Pauli string are omitted in ``axes``, but they can be inferred from the dimension of the Hilbert space ``dim``.
     """
 
+    dim: int
     axes: Mapping[int, Axis]
     sign: Sign = Sign.PLUS
 
@@ -104,7 +109,7 @@ def from_measured_node(flow: PauliFlow[Measurement], node: Node) -> PauliString:
         Returns
         -------
         PauliString
-            Primary extraction string associated to the input measured nodes.
+            Primary extraction string associated to the input measured nodes. The Pauli string is defined over qubit indices corresponding to positions in ``output_nodes``.
 
         Notes
         -----
@@ -115,6 +120,7 @@ def from_measured_node(flow: PauliFlow[Measurement], node: Node) -> PauliString:
         [1] Simmons, 2021 (arXiv:2109.05654).
         """
         og = flow.og
+        dim = len(flow.og.output_nodes)
         c_set = set(flow.correction_function[node])
         odd_c_set = og.odd_neighbors(c_set)
         inter_c_odd_set = c_set & odd_c_set
@@ -141,36 +147,23 @@ def from_measured_node(flow: PauliFlow[Measurement], node: Node) -> PauliString:
         # One phase flip if measured on the YZ plane.
         negative_sign ^= flow.node_measurement_label(node) == Plane.YZ
 
-        nodes = {}
-        # Sets `x_corrections`, `y_corrections` and `z_corrections` are disjoint.
-        for cnode in x_corrections:
-            nodes[cnode] = Axis.X
-        for cnode in y_corrections:
-            nodes[cnode] = Axis.Y
-        for cnode in z_corrections:
-            nodes[cnode] = Axis.Z
-        return PauliString(nodes, Sign.minus_if(negative_sign))
-
-    def remap(self, outputs_mapping: Callable[[int], int]) -> PauliString:
-        """Remap nodes to qubit indices.
+        axes_dict: dict[int, Axis] = {}
+        output_to_qubit_mapping = NodeIndex()
+        output_to_qubit_mapping.extend(og.output_nodes)
 
-        Parameters
-        ----------
-        outputs_mapping: Callable[[int], int]
-            Mapping between node numbers of the original MBQC pattern or open graph and qubit indices of a quantum circuit.
+        # Sets `x_corrections`, `y_corrections` and `z_corrections` are disjoint.
+        corrections = (x_corrections, y_corrections, z_corrections)
+        for correction, axis in zip(corrections, Axis, strict=True):
+            for cnode in correction:
+                qubit = output_to_qubit_mapping.index(cnode)
+                axes_dict[qubit] = axis
 
-        Returns
-        -------
-        PauliString
-            Pauli string defined on qubit indices.
-        """
-        axes = {outputs_mapping(n): axis for n, axis in self.axes.items()}
-        return PauliString(axes, self.sign)
+        return PauliString(dim, axes_dict, Sign.minus_if(negative_sign))
 
 
 @dataclass(frozen=True)
 class PauliExponential:
-    r"""Dataclass representing a Pauli exponential over a set of MBQC nodes.
+    r"""Dataclass representing a Pauli exponential.
 
     A Pauli exponential corresponds to the unitary operator
 
@@ -185,7 +178,7 @@ class PauliExponential:
     angle : ParameterizedAngle
         The Pauli exponential angle :math:`\alpha` in units of :math:`\pi`. When extracted from a corrected node, it corresponds to the node's measurement divided by two.
     pauli_string : PauliString
-        The signed Pauli string :math:`P` specifying the tensor product of Pauli operators acting on the corresponding MBQC nodes.
+        The signed Pauli string :math:`P` specifying the tensor product of Pauli operators acting on qubit indices.
     """
 
     angle: ParameterizedAngle
@@ -215,20 +208,13 @@ def from_measured_node(flow: PauliFlow[Measurement], node: Node) -> PauliExponen
         ----------
         [1] Simmons, 2021 (arXiv:2109.05654).
         """
-        pauli_string = flow.pauli_strings[node]
+        pauli_string = flow.extraction_pauli_strings[node]
         meas = flow.og.measurements[node]
         # We don't extract any rotation from Pauli Measurements. This is equivalent to setting the angle to 0.
         angle = meas.angle / 2 if isinstance(meas, BlochMeasurement) else 0
 
         return PauliExponential(angle, pauli_string)
 
-    def remap(self, outputs_mapping: Callable[[int], int]) -> Self:
-        """Remap nodes to qubit indices.
-
-        See documentation in :meth:`PauliString.remap` for additional information.
-        """
-        return replace(self, pauli_string=self.pauli_string.remap(outputs_mapping))
-
 
 @dataclass(frozen=True)
 class PauliExponentialDAG:
@@ -247,6 +233,9 @@ class PauliExponentialDAG:
     -----
     See Definition 3.3 in Ref. [1].
 
+    The Pauli strings in the Pauli exponentials are defined on qubit indices
+    which correspond to the indices of the sequence ``output_nodes``.
+
     References
     ----------
     [1] Simmons, 2021 (arXiv:2109.05654).
@@ -283,27 +272,38 @@ def from_focused_flow(flow: PauliFlow[Measurement]) -> PauliExponentialDAG:
 
         return PauliExponentialDAG(pauli_strings, flow.partial_order_layers, flow.og.output_nodes)
 
-    def remap(self, outputs_mapping: Callable[[int], int]) -> Self:
-        """Remap nodes to qubit indices.
-
-        See documentation in :meth:`PauliString.remap` for additional information.
-        """
-        pauli_exponentials = {node: pexp.remap(outputs_mapping) for node, pexp in self.pauli_exponentials.items()}
-        return replace(self, pauli_exponentials=pauli_exponentials)
-
 
 @dataclass(frozen=True)
 class CliffordMap:
-    """Dataclass to represent a Clifford map.
+    r"""Dataclass to represent a Clifford map.
+
+    A Clifford map encodes the action of a Clifford operator on Pauli generators.
+    It describes how single-qubit Pauli operators on the input qubits are mapped,
+    under conjugation, to Pauli strings on the output qubits.
 
-    A Clifford map describes a linear transformation between the space of input qubits and the space of output qubits. It is encoded as a map from the Pauli-group generators (X and Z) over the input nodes to Pauli strings over the output nodes.
+    For each input qubit :math:`i`, the map specifies:
+
+    .. math::
+
+        P_i = C X_i C^\dagger
+
+    .. math::
+
+        P_i = C Z_i C^\dagger
+
+    where the resulting operators :math:`P_i` are Pauli strings over the output qubits.
+
+    The sequences ``input_nodes`` and ``output_nodes`` define the correspondence
+    between qubit indices and node labels in the input and output spaces.
 
     Attributes
     ----------
-    x_map: Mapping[int, PauliString]
-        Map for the X generators. ``keys`` correspond to input nodes and ``values`` to their corresponding Pauli string over the outputs nodes.
-    z_map: Mapping[int, PauliString]
-        Map for the Z generators. ``keys`` correspond to input nodes and ``values`` to their corresponding Pauli string over the outputs nodes.
+    x_map: Sequence[PauliString]
+        Images of the :math:`X` generators. The :math:`i`-th element is the Pauli
+        string corresponding to :math:`C X_i C^\dagger`.
+    z_map: Sequence[PauliString]
+        Images of the :math:`Z` generators. The :math:`i`-th element is the Pauli
+        string corresponding to :math:`C Z_i C^\dagger`.
     input_nodes: Sequence[int]
         Sequence of inputs nodes.
     output_nodes: Sequence[int]
@@ -313,13 +313,15 @@ class CliffordMap:
     -----
     See Definition 3.3 in Ref. [1].
 
+    Elements of ``x_map`` and ``z_map`` are in one-to-one correspondance with ``input_nodes``. Each Pauli string is defined over qubit indices corresponding to positions in ``output_nodes``.
+
     References
     ----------
     [1] Simmons, 2021 (arXiv:2109.05654).
     """
 
-    x_map: Mapping[int, PauliString]
-    z_map: Mapping[int, PauliString]
+    x_map: Sequence[PauliString]
+    z_map: Sequence[PauliString]
     input_nodes: Sequence[int]
     output_nodes: Sequence[int]
 
@@ -350,24 +352,55 @@ def from_focused_flow(flow: PauliFlow[Measurement]) -> CliffordMap:
         x_map = clifford_x_map_from_focused_flow(flow)
         return CliffordMap(x_map, z_map, flow.og.input_nodes, flow.og.output_nodes)
 
-    def remap(self, inputs_mapping: Callable[[int], int], outputs_mapping: Callable[[int], int]) -> Self:
-        """Remap nodes to qubit indices.
+    def to_tableau(self) -> MatGF2:
+        """Convert the CliffordMap into its binary tableau representation.
 
-        Parameters
-        ----------
-        inputs_mapping: Callable[[int], int]
-            Mapping between input node numbers of the original MBQC pattern or open graph and qubit indices of a quantum circuit.
-        outputs_mapping: Callable[[int], int]
-            Mapping between output node numbers of the original MBQC pattern or open graph and qubit indices of a quantum circuit.
+        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``).
 
         Returns
         -------
-        CliffordMap
-            Clifford map defined on qubit indices.
+        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).
         """
-        x_map = {inputs_mapping(node): ps.remap(outputs_mapping) for node, ps in self.x_map.items()}
-        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)
+        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 (self.x_map, 0), (self.z_map, n):
+            for i, ps in enumerate(mapping):  # Indices in the Clifford map correspond to qubits (0 to n-1).
+                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]]:
@@ -409,8 +442,8 @@ def extend_input(og: OpenGraph[Measurement]) -> tuple[OpenGraph[Measurement], di
     return replace(og, graph=graph, input_nodes=new_input_nodes[::-1], measurements=measurements), ancillary_inputs_map
 
 
-def clifford_z_map_from_focused_flow(flow: PauliFlow[Measurement]) -> dict[int, PauliString]:
-    """Extract a map between Z over the input nodes and Pauli strings over the output nodes from a focused Pauli flow.
+def clifford_z_map_from_focused_flow(flow: PauliFlow[Measurement]) -> tuple[PauliString, ...]:
+    r"""Extract the images of the Z generators of a Clifford map from a focused Pauli flow.
 
     If the input node is a measured node, the resulting Pauli string is given by the correction set. If the input node is also an output node, the resulting Pauli string is Z (representing the identity map).
 
@@ -421,8 +454,9 @@ def clifford_z_map_from_focused_flow(flow: PauliFlow[Measurement]) -> dict[int,
 
     Returns
     -------
-    dict[int, PauliString]
-        Map between input nodes (``keys``) and Pauli strings over the output nodes (``values``).
+    tuple[PauliString,...]
+        Images of the :math:`Z` generators. The :math:`i`-th element is the Pauli string
+        corresponding to :math:`C Z_i C^\dagger`, where :math:`C` is the Clifford map.
 
     Notes
     -----
@@ -432,15 +466,20 @@ def clifford_z_map_from_focused_flow(flow: PauliFlow[Measurement]) -> dict[int,
     ----------
     [1] Simmons, 2021 (arXiv:2109.05654).
     """
-    # Nodes are either measured or outputs.
-    return {
-        node: flow.pauli_strings[node] if node in flow.og.measurements else PauliString({node: Axis.Z})
+    dim = len(flow.og.output_nodes)
+    output_to_qubit_mapping = NodeIndex()
+    output_to_qubit_mapping.extend(flow.og.output_nodes)
+    # Input nodes are either measured or outputs.
+    return tuple(
+        flow.extraction_pauli_strings[node]
+        if node in flow.og.measurements
+        else PauliString(dim, {output_to_qubit_mapping.index(node): Axis.Z})
         for node in flow.og.input_nodes
-    }
+    )
 
 
-def clifford_x_map_from_focused_flow(flow: PauliFlow[Measurement]) -> Mapping[int, PauliString]:
-    """Extract a map between X over the input nodes and Pauli strings over the output nodes from a focused Pauli flow.
+def clifford_x_map_from_focused_flow(flow: PauliFlow[Measurement]) -> tuple[PauliString, ...]:
+    r"""Extract the images of the X generators of a Clifford map from a focused Pauli flow.
 
     The resulting Pauli string is given by the correction set of a focused flow of the extended open graph.
 
@@ -451,8 +490,9 @@ def clifford_x_map_from_focused_flow(flow: PauliFlow[Measurement]) -> Mapping[in
 
     Returns
     -------
-    dict[int, PauliString]
-        Map between input nodes (``keys``) and Pauli strings over the output nodes (``values``).
+    tuple[PauliString,...]
+        Images of the :math:`X` generators. The :math:`i`-th element is the Pauli string
+        corresponding to :math:`C X_i C^\dagger`, where :math:`C` is the Clifford map.
 
     Notes
     -----
@@ -475,4 +515,4 @@ def clifford_x_map_from_focused_flow(flow: PauliFlow[Measurement]) -> Mapping[in
     # It's better to call the `PauliString` constructor instead of the cached property `flow_extended.pauli_strings` since the latter will compute a `PauliString` for _every_ node in the correction function and we just need it for the input nodes.
     x_map_ancillas = {node: PauliString.from_measured_node(flow_extended, node) for node in og_extended.input_nodes}
 
-    return {input_node: x_map_ancillas[ancillary_inputs_map[input_node]] for input_node in og.input_nodes}
+    return tuple(x_map_ancillas[ancillary_inputs_map[input_node]] for input_node in og.input_nodes)
diff --git a/graphix/flow/core.py b/graphix/flow/core.py
index f2e844ed0..be558dace 100644
--- a/graphix/flow/core.py
+++ b/graphix/flow/core.py
@@ -807,8 +807,8 @@ def is_focused(self) -> bool:
         return True
 
     @cached_property
-    def pauli_strings(self: PauliFlow[Measurement]) -> dict[int, PauliString]:
-        """Compute the Pauli strings associated with each node in the correction function.
+    def extraction_pauli_strings(self: PauliFlow[Measurement]) -> dict[int, PauliString]:
+        """Compute the extraction Pauli strings associated with each node in the correction function.
 
         This property requires the flow to be focused.
 
diff --git a/graphix/opengraph.py b/graphix/opengraph.py
index 717a6d5d2..1fc2e9f87 100644
--- a/graphix/opengraph.py
+++ b/graphix/opengraph.py
@@ -20,9 +20,11 @@
     # Unpack introduced in Python 3.12
     from typing_extensions import Unpack
 
+    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
     from graphix.visualization import DrawKwargs
 
 # TODO: Maybe move these definitions to graphix.fundamentals and graphix.measurements ? Now they are redefined in graphix.flow._find_gpflow, not very elegant.
@@ -559,6 +561,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.
 
diff --git a/noxfile.py b/noxfile.py
index 3418a52db..b9e3a0aae 100644
--- a/noxfile.py
+++ b/noxfile.py
@@ -119,7 +119,7 @@ class ReverseDependency:
             install_target=".[dev]",
         ),
         ReverseDependency("https://github.com/TeamGraphix/graphix-ibmq", doctest_modules=False),
-        ReverseDependency("https://github.com/qat-inria/graphix-stim-compiler"),
+        ReverseDependency("https://github.com/qat-inria/graphix-stim-compiler", branch="ps_dim"),
     ],
 )
 def tests_reverse_dependencies(session: Session, package: ReverseDependency) -> None:
diff --git a/tests/test_circ_extraction.py b/tests/test_circ_extraction.py
index 7ff891bf2..4116860ab 100644
--- a/tests/test_circ_extraction.py
+++ b/tests/test_circ_extraction.py
@@ -3,21 +3,25 @@
 from typing import TYPE_CHECKING, NamedTuple
 
 import networkx as nx
+import numpy as np
 import pytest
+from numpy.random import Generator
 
-from graphix.circ_ext.compilation import pexp_ladder_pass
-from graphix.circ_ext.extraction import PauliExponential, PauliExponentialDAG, PauliString, extend_input
+from graphix._linalg import MatGF2
+from graphix.circ_ext.compilation import cm_berg_pass, pexp_ladder_pass
+from graphix.circ_ext.extraction import CliffordMap, PauliExponential, PauliExponentialDAG, PauliString, extend_input
 from graphix.flow.core import PauliFlow
 from graphix.fundamentals import ANGLE_PI, Axis, Sign
-from graphix.instruction import CNOT, RX, RY, RZ, H
+from graphix.instruction import CNOT, RX, RY, RZ, H, S
 from graphix.measurements import Measurement
 from graphix.opengraph import OpenGraph
-from graphix.sim.base_backend import NodeIndex
+from graphix.parameter import Placeholder
+from graphix.random_objects import rand_circuit
 from graphix.states import BasicStates
 from graphix.transpiler import Circuit
 
 if TYPE_CHECKING:
-    from numpy.random import Generator
+    from numpy.random import PCG64
 
 
 class PauliExpTestCase(NamedTuple):
@@ -25,7 +29,7 @@ class PauliExpTestCase(NamedTuple):
     qc: Circuit
 
 
-class TestPauliExponential:
+class TestPauliExponentialDAG:
     # Angles of Pauli exponentials are in units of pi
     alpha = 0.3 * ANGLE_PI
 
@@ -35,7 +39,7 @@ class TestPauliExponential:
             PauliExpTestCase(
                 PauliExponentialDAG(
                     pauli_exponentials={
-                        0: PauliExponential(alpha / 2, PauliString({1: Axis.Z}, sign=Sign.MINUS)),
+                        0: PauliExponential(alpha / 2, PauliString(dim=1, axes={0: Axis.Z}, sign=Sign.MINUS)),
                     },
                     partial_order_layers=[{1}, {0}],
                     output_nodes=[1],
@@ -45,7 +49,7 @@ class TestPauliExponential:
             PauliExpTestCase(
                 PauliExponentialDAG(
                     pauli_exponentials={
-                        0: PauliExponential(alpha / 2, PauliString({1: Axis.X}, sign=Sign.MINUS)),
+                        0: PauliExponential(alpha / 2, PauliString(dim=1, axes={0: Axis.X}, sign=Sign.MINUS)),
                     },
                     partial_order_layers=[{1}, {0}],
                     output_nodes=[1],
@@ -55,7 +59,7 @@ class TestPauliExponential:
             PauliExpTestCase(
                 PauliExponentialDAG(
                     pauli_exponentials={
-                        0: PauliExponential(alpha / 2, PauliString({1: Axis.Y}, sign=Sign.MINUS)),
+                        0: PauliExponential(alpha / 2, PauliString(dim=1, axes={0: Axis.Y}, sign=Sign.MINUS)),
                     },
                     partial_order_layers=[{1}, {0}],
                     output_nodes=[1],
@@ -65,9 +69,9 @@ class TestPauliExponential:
             PauliExpTestCase(
                 PauliExponentialDAG(
                     pauli_exponentials={
-                        0: PauliExponential(ANGLE_PI / 4, PauliString({3: Axis.Z}, sign=Sign.MINUS)),
-                        1: PauliExponential(ANGLE_PI / 4, PauliString({3: Axis.X}, sign=Sign.MINUS)),
-                        2: PauliExponential(ANGLE_PI / 4, PauliString({3: Axis.Z}, sign=Sign.MINUS)),
+                        0: PauliExponential(ANGLE_PI / 4, PauliString(dim=1, axes={0: Axis.Z}, sign=Sign.MINUS)),
+                        1: PauliExponential(ANGLE_PI / 4, PauliString(dim=1, axes={0: Axis.X}, sign=Sign.MINUS)),
+                        2: PauliExponential(ANGLE_PI / 4, PauliString(dim=1, axes={0: Axis.Z}, sign=Sign.MINUS)),
                     },
                     partial_order_layers=[{3}, {2}, {1}, {0}],
                     output_nodes=[3],
@@ -77,9 +81,11 @@ class TestPauliExponential:
             PauliExpTestCase(
                 PauliExponentialDAG(
                     pauli_exponentials={
-                        0: PauliExponential(ANGLE_PI / 4, PauliString({3: Axis.X})),
-                        1: PauliExponential(ANGLE_PI / 4, PauliString({5: Axis.Z})),
-                        2: PauliExponential(ANGLE_PI / 4, PauliString({3: Axis.X, 5: Axis.Z}, Sign.MINUS)),
+                        0: PauliExponential(ANGLE_PI / 4, PauliString(dim=2, axes={1: Axis.X})),
+                        1: PauliExponential(ANGLE_PI / 4, PauliString(dim=2, axes={0: Axis.Z})),
+                        2: PauliExponential(
+                            ANGLE_PI / 4, PauliString(dim=2, axes={1: Axis.X, 0: Axis.Z}, sign=Sign.MINUS)
+                        ),
                     },
                     partial_order_layers=[{5, 3}, {2}, {0, 1}],
                     output_nodes=[5, 3],  # Node 5 -> qubit 0 (control), node 3 -> qubit 1 (target)
@@ -89,7 +95,7 @@ class TestPauliExponential:
             PauliExpTestCase(
                 PauliExponentialDAG(
                     pauli_exponentials={
-                        0: PauliExponential(alpha / 2, PauliString({1: Axis.X, 2: Axis.Z, 4: Axis.Z})),
+                        0: PauliExponential(alpha / 2, PauliString(dim=4, axes={1: Axis.X, 0: Axis.Z, 3: Axis.Z})),
                     },
                     partial_order_layers=[{1, 2, 3, 4}, {0}],
                     output_nodes=[2, 1, 3, 4],
@@ -100,33 +106,36 @@ class TestPauliExponential:
     )
     def test_to_circuit(self, test_case: PauliExpTestCase, fx_rng: Generator) -> None:
         qc = Circuit(len(test_case.pexp_dag.output_nodes))
-        outputs_mapping = NodeIndex()
-        outputs_mapping.extend(test_case.pexp_dag.output_nodes)
-        pexp_ladder_pass(test_case.pexp_dag.remap(outputs_mapping.index), qc)
+        pexp_ladder_pass(test_case.pexp_dag, qc)
         state = qc.simulate_statevector(rng=fx_rng).statevec
         state_ref = test_case.qc.simulate_statevector(rng=fx_rng).statevec
         assert state.isclose(state_ref)
 
     def test_to_circuit_outputs_order(self, fx_rng: Generator) -> None:
-        pexp_map = {2: PauliExponential(0.1, PauliString({1: Axis.X, 0: Axis.Z}))}
         pol = [{0, 1}, {2}]
 
+        # PauliString is defined as a mapping from qubit to axes.
+        # We represent the Pauli string Z(node=0)X(node=1) in both cases.
+
+        pexp_map_1 = {2: PauliExponential(0.1, PauliString(dim=2, axes={0: Axis.Z, 1: Axis.X}))}
         outputs_1 = [0, 1]
+        pexp_dag_1 = PauliExponentialDAG(
+            pauli_exponentials=pexp_map_1, partial_order_layers=pol, output_nodes=outputs_1
+        )
+
+        pexp_map_2 = {2: PauliExponential(0.1, PauliString(dim=2, axes={0: Axis.X, 1: Axis.Z}))}
         outputs_2 = [1, 0]
+        pexp_dag_2 = PauliExponentialDAG(
+            pauli_exponentials=pexp_map_2, partial_order_layers=pol, output_nodes=outputs_2
+        )
 
-        pexp_dag_1 = PauliExponentialDAG(pauli_exponentials=pexp_map, partial_order_layers=pol, output_nodes=outputs_1)
         qc_1 = Circuit(2)
-        outputs_mapping_1 = NodeIndex()
-        outputs_mapping_1.extend(pexp_dag_1.output_nodes)
-        pexp_ladder_pass(pexp_dag_1.remap(outputs_mapping_1.index), qc_1)
+        pexp_ladder_pass(pexp_dag_1, qc_1)
         s_1 = qc_1.simulate_statevector(rng=fx_rng, input_state=[BasicStates.PLUS, BasicStates.MINUS]).statevec
 
-        pexp_dag_2 = PauliExponentialDAG(pauli_exponentials=pexp_map, partial_order_layers=pol, output_nodes=outputs_2)
         qc_2 = Circuit(2)
         qc_2.swap(0, 1)  # We must swap before and after the Pauli exponential!
-        outputs_mapping_2 = NodeIndex()
-        outputs_mapping_2.extend(pexp_dag_2.output_nodes)
-        pexp_ladder_pass(pexp_dag_2.remap(outputs_mapping_2.index), qc_2)
+        pexp_ladder_pass(pexp_dag_2, qc_2)
 
         s_2 = qc_2.simulate_statevector(rng=fx_rng, input_state=[BasicStates.PLUS, BasicStates.MINUS]).statevec
         assert not s_1.isclose(s_2)
@@ -170,12 +179,14 @@ def test_from_focused_flow(self) -> None:
 
         pexp_dag_ref = PauliExponentialDAG(
             pauli_exponentials={
-                0: PauliExponential(ANGLE_PI * 0.1 / 2, PauliString({6: Axis.X})),
-                1: PauliExponential(ANGLE_PI * 0.2 / 2, PauliString({6: Axis.Y, 5: Axis.Z})),
-                2: PauliExponential(ANGLE_PI * 0.3 / 2, PauliString({5: Axis.Y, 6: Axis.Z}, Sign.MINUS)),
-                3: PauliExponential(ANGLE_PI * 0.4 / 2, PauliString({5: Axis.X})),
+                0: PauliExponential(ANGLE_PI * 0.1 / 2, PauliString(dim=2, axes={1: Axis.X})),
+                1: PauliExponential(ANGLE_PI * 0.2 / 2, PauliString(dim=2, axes={1: Axis.Y, 0: Axis.Z})),
+                2: PauliExponential(
+                    ANGLE_PI * 0.3 / 2, PauliString(dim=2, axes={0: Axis.Y, 1: Axis.Z}, sign=Sign.MINUS)
+                ),
+                3: PauliExponential(ANGLE_PI * 0.4 / 2, PauliString(dim=2, axes={0: Axis.X})),
                 4: PauliExponential(
-                    0, PauliString({6: Axis.X})
+                    0, PauliString(dim=2, axes={1: Axis.X})
                 ),  # The angle is 0 (interpreted from the Pauli measurement).
             },
             partial_order_layers=flow.partial_order_layers,
@@ -185,6 +196,413 @@ def test_from_focused_flow(self) -> None:
         assert pexp_dag == pexp_dag_ref
 
 
+class TestCliffordMap:
+    @pytest.mark.parametrize(
+        ("cm", "tab_ref"),
+        [
+            (
+                CliffordMap(
+                    x_map=(PauliString(dim=2, axes={0: Axis.Z}), PauliString(dim=2, axes={1: Axis.Y})),
+                    z_map=(
+                        PauliString(dim=2, axes={0: Axis.X, 1: Axis.Y}, sign=Sign.MINUS),
+                        PauliString(dim=2, axes={0: Axis.Z, 1: Axis.Z}),
+                    ),
+                    input_nodes=[0, 1],
+                    output_nodes=[0, 1],
+                ),
+                MatGF2([[0, 0, 1, 0, 0], [0, 1, 0, 1, 0], [1, 1, 0, 1, 1], [0, 0, 1, 1, 0]]),
+            ),
+            (
+                CliffordMap(
+                    x_map=(
+                        PauliString(dim=3, axes={0: Axis.Z}),
+                        PauliString(dim=3, axes={1: Axis.X}),
+                        PauliString(dim=3, axes={2: Axis.Y}),
+                    ),
+                    z_map=(
+                        PauliString(dim=3, axes={0: Axis.X, 1: Axis.X}),
+                        PauliString(dim=3, axes={0: Axis.Z, 1: Axis.Z}),
+                        PauliString(dim=3, axes={2: Axis.Z}),
+                    ),
+                    input_nodes=[0, 1, 2],
+                    output_nodes=[0, 1, 2],
+                ),
+                MatGF2(
+                    [
+                        [0, 0, 0, 1, 0, 0, 0],
+                        [0, 1, 0, 0, 0, 0, 0],
+                        [0, 0, 1, 0, 0, 1, 0],
+                        [1, 1, 0, 0, 0, 0, 0],
+                        [0, 0, 0, 1, 1, 0, 0],
+                        [0, 0, 0, 0, 0, 1, 0],
+                    ]
+                ),
+            ),
+        ],
+    )
+    def test_to_tableau(self, cm: CliffordMap, tab_ref: MatGF2) -> None:
+        tab = cm.to_tableau()
+        assert np.all(tab == tab_ref)
+
+    # The CliffordMap test cases were generated using conversion tools in the graphix-stim-compiler plugin and stim.Tableau.random.
+    @pytest.mark.parametrize(
+        ("qc_ref", "cm"),
+        [
+            (
+                Circuit(width=1, instr=[H(0)]),
+                CliffordMap(
+                    x_map=(PauliString(dim=1, axes={0: Axis.Z}, sign=Sign.PLUS),),
+                    z_map=(PauliString(dim=1, axes={0: Axis.X}, sign=Sign.PLUS),),
+                    input_nodes=[0],
+                    output_nodes=[0],
+                ),
+            ),
+            (
+                Circuit(width=1, instr=[S(0)]),
+                CliffordMap(
+                    x_map=(PauliString(dim=1, axes={0: Axis.Y}, sign=Sign.PLUS),),
+                    z_map=(PauliString(dim=1, axes={0: Axis.Z}, sign=Sign.PLUS),),
+                    input_nodes=[0],
+                    output_nodes=[0],
+                ),
+            ),
+            (
+                Circuit(width=2, instr=[CNOT(1, 0)]),
+                CliffordMap(
+                    x_map=(
+                        PauliString(dim=2, axes={0: Axis.X, 1: Axis.X}, sign=Sign.PLUS),
+                        PauliString(dim=2, axes={1: Axis.X}, sign=Sign.PLUS),
+                    ),
+                    z_map=(
+                        PauliString(dim=2, axes={0: Axis.Z}, sign=Sign.PLUS),
+                        PauliString(dim=2, axes={0: Axis.Z, 1: Axis.Z}, sign=Sign.PLUS),
+                    ),
+                    input_nodes=[0, 1],
+                    output_nodes=[0, 1],
+                ),
+            ),
+            (
+                Circuit(width=3, instr=[CNOT(1, 0), H(0), H(1), CNOT(2, 1), S(1), CNOT(2, 0), H(2), S(2)]),
+                CliffordMap(
+                    x_map=(
+                        PauliString(dim=3, axes={0: Axis.Z, 1: Axis.Z}, sign=Sign.PLUS),
+                        PauliString(dim=3, axes={1: Axis.Z}, sign=Sign.PLUS),
+                        PauliString(dim=3, axes={2: Axis.Z}, sign=Sign.PLUS),
+                    ),
+                    z_map=(
+                        PauliString(dim=3, axes={0: Axis.X, 2: Axis.Z}, sign=Sign.PLUS),
+                        PauliString(dim=3, axes={0: Axis.X, 1: Axis.Y}, sign=Sign.PLUS),
+                        PauliString(dim=3, axes={0: Axis.Z, 1: Axis.Z, 2: Axis.Y}, sign=Sign.PLUS),
+                    ),
+                    input_nodes=[0, 1, 2],
+                    output_nodes=[0, 1, 2],
+                ),
+            ),
+            (
+                Circuit(width=1, instr=[S(0), S(0), S(0)]),
+                CliffordMap(
+                    x_map=(PauliString(dim=1, axes={0: Axis.Y}, sign=Sign.MINUS),),
+                    z_map=(PauliString(dim=1, axes={0: Axis.Z}, sign=Sign.PLUS),),
+                    input_nodes=[0],
+                    output_nodes=[0],
+                ),
+            ),
+            (
+                Circuit(width=2, instr=[CNOT(1, 0), H(0), S(0), S(0), H(0), S(1), S(1)]),
+                CliffordMap(
+                    x_map=(
+                        PauliString(dim=2, axes={0: Axis.X, 1: Axis.X}, sign=Sign.MINUS),
+                        PauliString(dim=2, axes={1: Axis.X}, sign=Sign.MINUS),
+                    ),
+                    z_map=(
+                        PauliString(dim=2, axes={0: Axis.Z}, sign=Sign.MINUS),
+                        PauliString(dim=2, axes={0: Axis.Z, 1: Axis.Z}, sign=Sign.MINUS),
+                    ),
+                    input_nodes=[0, 1],
+                    output_nodes=[0, 1],
+                ),
+            ),
+            (
+                Circuit(
+                    width=3,
+                    instr=[
+                        S(0),
+                        H(2),
+                        CNOT(2, 0),
+                        H(1),
+                        H(2),
+                        CNOT(0, 1),
+                        CNOT(0, 2),
+                        S(2),
+                        CNOT(2, 1),
+                        H(2),
+                        CNOT(1, 2),
+                        H(2),
+                        H(1),
+                        S(1),
+                        S(1),
+                        H(1),
+                        S(0),
+                        S(0),
+                        S(1),
+                        S(1),
+                    ],
+                ),
+                CliffordMap(
+                    x_map=(
+                        PauliString(dim=3, axes={0: Axis.Y, 1: Axis.Z, 2: Axis.X}, sign=Sign.PLUS),
+                        PauliString(dim=3, axes={1: Axis.Z, 2: Axis.X}, sign=Sign.MINUS),
+                        PauliString(dim=3, axes={0: Axis.Y, 1: Axis.Z}, sign=Sign.PLUS),
+                    ),
+                    z_map=(
+                        PauliString(dim=3, axes={0: Axis.Z, 1: Axis.X, 2: Axis.Z}, sign=Sign.MINUS),
+                        PauliString(dim=3, axes={0: Axis.X, 1: Axis.X, 2: Axis.X}, sign=Sign.PLUS),
+                        PauliString(dim=3, axes={1: Axis.Y, 2: Axis.Y}, sign=Sign.PLUS),
+                    ),
+                    input_nodes=[0, 1, 2],
+                    output_nodes=[0, 1, 2],
+                ),
+            ),
+            (
+                Circuit(
+                    width=4,
+                    instr=[
+                        CNOT(0, 3),
+                        CNOT(3, 0),
+                        CNOT(0, 3),
+                        S(0),
+                        H(0),
+                        S(0),
+                        H(2),
+                        CNOT(2, 0),
+                        H(1),
+                        H(3),
+                        CNOT(0, 1),
+                        CNOT(0, 3),
+                        S(1),
+                        S(2),
+                        S(3),
+                        CNOT(2, 1),
+                        CNOT(3, 1),
+                        S(3),
+                        H(3),
+                        CNOT(1, 2),
+                        CNOT(1, 3),
+                        CNOT(2, 3),
+                        CNOT(3, 2),
+                        CNOT(2, 3),
+                        S(2),
+                        S(3),
+                        H(3),
+                        S(3),
+                        H(0),
+                        S(0),
+                        S(0),
+                        H(0),
+                        S(0),
+                        S(0),
+                        S(1),
+                        S(1),
+                        S(2),
+                        S(2),
+                    ],
+                ),
+                CliffordMap(
+                    x_map=(
+                        PauliString(dim=4, axes={1: Axis.Y, 2: Axis.X, 3: Axis.Y}, sign=Sign.PLUS),
+                        PauliString(dim=4, axes={1: Axis.Z, 2: Axis.Z, 3: Axis.Y}, sign=Sign.PLUS),
+                        PauliString(dim=4, axes={0: Axis.Z, 1: Axis.Y, 2: Axis.X}, sign=Sign.MINUS),
+                        PauliString(dim=4, axes={0: Axis.X, 1: Axis.Y, 2: Axis.Z, 3: Axis.X}, sign=Sign.PLUS),
+                    ),
+                    z_map=(
+                        PauliString(dim=4, axes={0: Axis.X, 1: Axis.Z, 3: Axis.Y}, sign=Sign.PLUS),
+                        PauliString(dim=4, axes={0: Axis.X, 1: Axis.Y, 2: Axis.Y, 3: Axis.Z}, sign=Sign.MINUS),
+                        PauliString(dim=4, axes={1: Axis.Y, 2: Axis.Z, 3: Axis.X}, sign=Sign.MINUS),
+                        PauliString(dim=4, axes={0: Axis.Y, 1: Axis.Z, 2: Axis.X, 3: Axis.X}, sign=Sign.PLUS),
+                    ),
+                    input_nodes=[0, 1, 2, 3],
+                    output_nodes=[0, 1, 2, 3],
+                ),
+            ),
+        ],
+    )
+    def test_cm_berg_pass(self, qc_ref: Circuit, cm: CliffordMap, fx_rng: Generator) -> None:
+
+        qc = Circuit(qc_ref.width)
+        cm_berg_pass(cm, qc)
+
+        s_test = qc.simulate_statevector(rng=fx_rng).statevec
+        s_ref = qc_ref.simulate_statevector(rng=fx_rng).statevec
+
+        assert s_test.isclose(s_ref)
+
+
+class TestExtraction:
+    @pytest.mark.parametrize("jumps", range(1, 11))
+    def test_extract_rnd_circuit(self, fx_bg: PCG64, jumps: int) -> None:
+        rng = Generator(fx_bg.jumped(jumps))
+        nqubits = 2
+        depth = 2
+        circuit_ref = rand_circuit(nqubits, depth, rng, use_ccx=False)
+        pattern = circuit_ref.transpile().pattern
+
+        circuit = pattern.extract_opengraph().extract_circuit()
+
+        s_ref = circuit.simulate_statevector(rng=rng).statevec
+        s_test = circuit_ref.simulate_statevector(rng=rng).statevec
+        assert s_ref.isclose(s_test)
+
+    @pytest.mark.parametrize(
+        "test_case",
+        [
+            OpenGraph(
+                graph=nx.Graph([(1, 3), (2, 4), (3, 4), (3, 5), (4, 6)]),
+                input_nodes=[1, 2],
+                output_nodes=[5, 6],
+                measurements={
+                    1: Measurement.XY(0.1),
+                    2: Measurement.XY(0.2),
+                    3: Measurement.XY(0.3),
+                    4: Measurement.XY(0.4),
+                },
+            ),
+            OpenGraph(
+                graph=nx.Graph([(1, 4), (1, 6), (2, 4), (2, 5), (2, 6), (3, 5), (3, 6)]),
+                input_nodes=[1, 2, 3],
+                output_nodes=[4, 5, 6],
+                measurements={
+                    1: Measurement.XY(0.1),
+                    2: Measurement.XY(0.2),
+                    3: Measurement.XY(0.3),
+                },
+            ),
+            OpenGraph(
+                graph=nx.Graph([(0, 1), (0, 2), (0, 4), (1, 5), (2, 4), (2, 5), (3, 5)]),
+                input_nodes=[0, 1],
+                output_nodes=[4, 5],
+                measurements={
+                    0: Measurement.XY(0.1),
+                    1: Measurement.XY(0.1),
+                    2: Measurement.XZ(0.2),
+                    3: Measurement.YZ(0.3),
+                },
+            ),
+            OpenGraph(
+                graph=nx.Graph([(0, 1), (1, 2), (1, 4), (2, 3)]),
+                input_nodes=[0],
+                output_nodes=[4],
+                measurements={
+                    0: Measurement.XY(0.1),  # XY
+                    1: Measurement.X,  # X
+                    2: Measurement.XY(0.1),  # XY
+                    3: Measurement.X,  # X
+                },
+            ),
+            OpenGraph(
+                graph=nx.Graph([(0, 1), (1, 2)]),
+                input_nodes=[0],
+                output_nodes=[2],
+                measurements={
+                    0: Measurement.XY(0.1),  # XY
+                    1: Measurement.Y,  # Y
+                },
+            ),
+            OpenGraph(
+                graph=nx.Graph([(0, 2), (1, 3), (2, 3), (2, 4), (3, 5), (4, 5), (4, 6), (5, 7)]),
+                input_nodes=[0, 1],
+                output_nodes=[7, 6],
+                measurements={
+                    0: Measurement.XY(0.1),  # XY
+                    1: Measurement.XY(0.1),  # XY
+                    2: Measurement.X,  # X
+                    3: Measurement.XY(0.1),  # XY
+                    4: Measurement.X,  # X
+                    5: Measurement.Y,  # Y
+                },
+            ),
+        ],
+    )
+    def test_extract_og(self, test_case: OpenGraph[Measurement], fx_rng: Generator) -> None:
+        pattern = test_case.to_pattern()
+        circuit = test_case.extract_circuit()
+
+        state = circuit.simulate_statevector(rng=fx_rng).statevec
+        state_ref = pattern.simulate_pattern(rng=fx_rng)
+        assert state.isclose(state_ref)
+
+    @pytest.mark.parametrize("infer_pauli", [True, False])
+    def test_extract_og_infer_pauli(self, infer_pauli: bool, fx_rng: Generator) -> None:
+        og: OpenGraph[Measurement] = OpenGraph(
+            graph=nx.Graph([(0, 2), (1, 3), (2, 3), (2, 4), (3, 5), (4, 5), (4, 6), (5, 7)]),
+            input_nodes=[0, 1],
+            output_nodes=[6, 7],
+            measurements={
+                0: Measurement.XY(0),
+                1: Measurement.XY(0),
+                2: Measurement.XY(0),
+                3: Measurement.XY(0),
+                4: Measurement.XY(0),
+                5: Measurement.XY(0),
+            },
+        )
+        pattern = og.to_pattern()
+        if infer_pauli:
+            og = og.infer_pauli_measurements()
+
+        circuit = og.extract_circuit()
+
+        state = circuit.simulate_statevector(rng=fx_rng).statevec
+        state_ref = pattern.simulate_pattern(rng=fx_rng)
+        assert state.isclose(state_ref)
+
+    def test_extract_og_gflow(self, fx_rng: Generator) -> None:
+        og = OpenGraph(
+            graph=nx.Graph([(1, 3), (2, 4), (3, 4), (3, 5), (4, 6)]),
+            input_nodes=[1, 2],
+            output_nodes=[5, 6],
+            measurements={
+                1: Measurement.XY(0.1),
+                2: Measurement.XY(0.2),
+                3: Measurement.XY(0.3),
+                4: Measurement.XY(0.4),
+            },
+        )
+        pattern = og.to_pattern()
+        circuit = og.extract_gflow().extract_circuit().to_circuit()
+
+        state = circuit.simulate_statevector(rng=fx_rng).statevec
+        state_ref = pattern.simulate_pattern(rng=fx_rng)
+        assert state.isclose(state_ref)
+
+    @pytest.mark.parametrize("test_case", [0.2, 0.5, 1.0])
+    def test_parametric_angles(self, test_case: float, fx_rng: Generator) -> None:
+        alpha = Placeholder("alpha")
+        alpha_val = test_case
+        og = OpenGraph(
+            graph=nx.Graph([(1, 3), (2, 4), (3, 4), (3, 5), (4, 6)]),
+            input_nodes=[1, 2],
+            output_nodes=[5, 6],
+            measurements={
+                1: Measurement.XY(0.1),
+                2: Measurement.XY(alpha),
+                3: Measurement.XY(0.3),
+                4: Measurement.XY(alpha),
+            },
+        )
+
+        # Substitute parameter at the level of the extracted circuit
+        qc1 = og.extract_circuit()
+        s1 = qc1.subs(alpha, alpha_val).simulate_statevector(rng=fx_rng).statevec
+
+        # Substitute parameter at the level of the open graph object
+        # Calling `infer_pauli_measurements` is not necessary for the test to pass
+        # (and it should not be), but it suppresses the warnings.
+        qc2 = og.subs(alpha, alpha_val).infer_pauli_measurements().extract_circuit()
+        s2 = qc2.simulate_statevector(rng=fx_rng).statevec
+
+        assert s1.isclose(s2)
+
+
 def test_extend_input() -> None:
     og = OpenGraph(
         graph=nx.Graph([(1, 3), (2, 4), (3, 4), (3, 5), (4, 6)]),