From d6d8f96595760804fac7bf971ec96c59468168d2 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tiberiu=20Sab=C4=83u?= <T.Sabau@student.tudelft.nl>
Date: Thu, 5 Feb 2026 17:42:55 +0100
Subject: [PATCH 01/20] Make tests use the BranchSelector feature

---
 tests/test_noisy_density_matrix.py | 196 ++++++++++++++++++-----------
 1 file changed, 125 insertions(+), 71 deletions(-)

diff --git a/tests/test_noisy_density_matrix.py b/tests/test_noisy_density_matrix.py
index 20c5cf91e..87e41828d 100644
--- a/tests/test_noisy_density_matrix.py
+++ b/tests/test_noisy_density_matrix.py
@@ -6,6 +6,7 @@
 import numpy.typing as npt
 import pytest
 
+from graphix.branch_selector import ConstBranchSelector, FixedBranchSelector
 from graphix.fundamentals import angle_to_rad
 from graphix.noise_models import DepolarisingNoiseModel
 from graphix.noise_models.noise_model import NoiselessNoiseModel
@@ -69,14 +70,22 @@ def test_noisy_measure_confuse_hadamard(self, fx_rng: Generator) -> None:
         assert isinstance(res, DensityMatrix)
         assert np.allclose(res.rho, np.array([[0.0, 0.0], [0.0, 1.0]]))
 
-        # arbitrary probability
+    @pytest.mark.parametrize("outcome", [0, 1])
+    def test_noisy_measure_confuse_hadamard_arbitrary(self, fx_rng: Generator, outcome: int) -> None:
+        # arbitrary probability with fixed branch
+        hadamardpattern = hpat()
         measure_error_pr = fx_rng.random()
-        print(f"measure_error_pr = {measure_error_pr}")
+        print(f"measure_error_pr = {measure_error_pr}, outcome = {outcome}")
         res = hadamardpattern.simulate_pattern(
-            backend="densitymatrix", noise_model=DepolarisingNoiseModel(measure_error_prob=measure_error_pr), rng=fx_rng
+            backend="densitymatrix",
+            noise_model=DepolarisingNoiseModel(measure_error_prob=measure_error_pr),
+            branch_selector=ConstBranchSelector(outcome),
+            rng=fx_rng,
         )
-        # result should be |1>
         assert isinstance(res, DensityMatrix)
+        # With measure_error_prob, the outcome might be flipped, resulting in different X corrections
+        # However, we cannot predict the exact result without knowing if the error occurred
+        # So we check both possibilities
         assert np.allclose(res.rho, np.array([[1.0, 0.0], [0.0, 0.0]])) or np.allclose(
             res.rho,
             np.array([[0.0, 0.0], [0.0, 1.0]]),
@@ -100,21 +109,21 @@ def test_noisy_measure_channel_hadamard(self, fx_rng: Generator) -> None:
         )
 
     # test Pauli X error
-    def test_noisy_x_hadamard(self, fx_rng: Generator) -> None:
+    @pytest.mark.parametrize("outcome", [0, 1])
+    def test_noisy_x_hadamard(self, fx_rng: Generator, outcome: int) -> None:
         hadamardpattern = hpat()
         # x error only
         x_error_pr = fx_rng.random()
-        print(f"x_error_pr = {x_error_pr}")
+        print(f"x_error_pr = {x_error_pr}, outcome = {outcome}")
         res = hadamardpattern.simulate_pattern(
             backend="densitymatrix",
             noise_model=DepolarisingNoiseModel(x_error_prob=x_error_pr),
+            branch_selector=ConstBranchSelector(outcome),
             rng=fx_rng,
         )
-        # analytical result since deterministic pattern output is |0>.
-        # if no X applied, no noise. If X applied X noise on |0><0|
-
+        # Both outcomes lead to |0> after correction. ApplyNoise is unconditional, so noise is always applied.
         assert isinstance(res, DensityMatrix)
-        assert np.allclose(res.rho, np.array([[1.0, 0.0], [0.0, 0.0]])) or np.allclose(
+        assert np.allclose(
             res.rho,
             np.array([[1 - 2 * x_error_pr / 3.0, 0.0], [0.0, 2 * x_error_pr / 3.0]]),
         )
@@ -310,26 +319,36 @@ def test_noisy_measure_channel_rz(self, fx_rng: Generator) -> None:
             ),
         )
 
-    def test_noisy_x_rz(self, fx_rng: Generator) -> None:
+    @pytest.mark.parametrize("outcome_z,outcome_x", [(0, 0), (0, 1), (1, 0), (1, 1)])
+    def test_noisy_x_rz(self, fx_rng: Generator, outcome_z: int, outcome_x: int) -> None:
         alpha = fx_rng.random()
         rzpattern = rzpat(alpha)
         # x error only
         x_error_pr = fx_rng.random()
-        print(f"x_error_pr = {x_error_pr}")
+        print(f"x_error_pr = {x_error_pr}, outcome_z = {outcome_z}, outcome_x = {outcome_x}")
+
+        # M(0) determines Z, M(1) determines X
+        results = {}
+        cmd_count = 0
+        for cmd in rzpattern:
+            if cmd.kind.name == "M":
+                if cmd_count == 0:
+                    results[cmd.node] = outcome_z
+                elif cmd_count == 1:
+                    results[cmd.node] = outcome_x
+                cmd_count += 1
+
         res = rzpattern.simulate_pattern(
             backend="densitymatrix",
             noise_model=DepolarisingNoiseModel(x_error_prob=x_error_pr),
+            branch_selector=FixedBranchSelector(results),
             rng=fx_rng,
         )
 
-        # only two cases: if no X correction, Z or no Z correction but exact result.
-        # If X correction the noise result is the same with or without the PERFECT Z correction.
+        # All outcomes lead to the same state after corrections. ApplyNoise is unconditional, so noise is always applied.
         assert isinstance(res, DensityMatrix)
         rad = angle_to_rad(alpha)
         assert np.allclose(
-            res.rho,
-            0.5 * np.array([[1.0, np.exp(-1j * rad)], [np.exp(1j * rad), 1.0]]),
-        ) or np.allclose(
             res.rho,
             0.5
             * np.array(
@@ -340,26 +359,36 @@ def test_noisy_x_rz(self, fx_rng: Generator) -> None:
             ),
         )
 
-    def test_noisy_z_rz(self, fx_rng: Generator) -> None:
+    @pytest.mark.parametrize("outcome_z,outcome_x", [(0, 0), (0, 1), (1, 0), (1, 1)])
+    def test_noisy_z_rz(self, fx_rng: Generator, outcome_z: int, outcome_x: int) -> None:
         alpha = fx_rng.random()
         rzpattern = rzpat(alpha)
         # z error only
         z_error_pr = fx_rng.random()
-        print(f"z_error_pr = {z_error_pr}")
+        print(f"z_error_pr = {z_error_pr}, outcome_z = {outcome_z}, outcome_x = {outcome_x}")
+
+        # M(0) determines Z, M(1) determines X
+        results = {}
+        cmd_count = 0
+        for cmd in rzpattern:
+            if cmd.kind.name == "M":
+                if cmd_count == 0:
+                    results[cmd.node] = outcome_z
+                elif cmd_count == 1:
+                    results[cmd.node] = outcome_x
+                cmd_count += 1
+
         res = rzpattern.simulate_pattern(
             backend="densitymatrix",
             noise_model=DepolarisingNoiseModel(z_error_prob=z_error_pr),
+            branch_selector=FixedBranchSelector(results),
             rng=fx_rng,
         )
 
-        # only two cases: if no Z correction, X or no X correction but exact result.
-        # If Z correction the noise result is the same with or without the PERFECT X correction.
+        # All outcomes lead to the same state after corrections. ApplyNoise is unconditional, so noise is always applied.
         assert isinstance(res, DensityMatrix)
         rad = angle_to_rad(alpha)
         assert np.allclose(
-            res.rho,
-            0.5 * np.array([[1.0, np.exp(-1j * rad)], [np.exp(1j * rad), 1.0]]),
-        ) or np.allclose(
             res.rho,
             0.5
             * np.array(
@@ -370,7 +399,8 @@ def test_noisy_z_rz(self, fx_rng: Generator) -> None:
             ),
         )
 
-    def test_noisy_xz_rz(self, fx_rng: Generator) -> None:
+    @pytest.mark.parametrize("z_outcome,x_outcome", [(0, 0), (0, 1), (1, 0), (1, 1)])
+    def test_noisy_xz_rz(self, fx_rng: Generator, z_outcome: int, x_outcome: int) -> None:
         alpha = fx_rng.random()
         rzpattern = rzpat(alpha)
         # x and z errors
@@ -378,78 +408,102 @@ def test_noisy_xz_rz(self, fx_rng: Generator) -> None:
         print(f"x_error_pr = {x_error_pr}")
         z_error_pr = fx_rng.random()
         print(f"z_error_pr = {z_error_pr}")
+        print(f"z_outcome = {z_outcome}, x_outcome = {x_outcome}")
+
+        # M(0) determines Z correction, M(1) determines X correction
+        results = {}
+        cmd_count = 0
+        for cmd in rzpattern:
+            if cmd.kind.name == "M":
+                if cmd_count == 0:
+                    results[cmd.node] = z_outcome
+                elif cmd_count == 1:
+                    results[cmd.node] = x_outcome
+                cmd_count += 1
+
         res = rzpattern.simulate_pattern(
             backend="densitymatrix",
             noise_model=DepolarisingNoiseModel(x_error_prob=x_error_pr, z_error_prob=z_error_pr),
+            branch_selector=FixedBranchSelector(results),
             rng=fx_rng,
         )
 
-        # 4 cases : no corr, noisy X, noisy Z, noisy XZ.
+        # All outcomes lead to same state after corrections. Both X and Z noise applied unconditionally.
         assert isinstance(res, DensityMatrix)
         rad = angle_to_rad(alpha)
-        assert (
-            np.allclose(res.rho, 0.5 * np.array([[1.0, np.exp(-1j * rad)], [np.exp(1j * rad), 1.0]]))
-            or np.allclose(
-                res.rho,
-                0.5
-                * np.array(
-                    [
-                        [1.0, np.exp(-1j * rad) * (3 - 4 * x_error_pr) * (3 - 4 * z_error_pr) / 9],
-                        [np.exp(1j * rad) * (3 - 4 * x_error_pr) * (3 - 4 * z_error_pr) / 9, 1.0],
-                    ],
-                ),
-            )
-            or np.allclose(
-                res.rho,
-                0.5
-                * np.array(
-                    [
-                        [1.0, np.exp(-1j * rad) * (3 - 4 * z_error_pr) / 3],
-                        [np.exp(1j * rad) * (3 - 4 * z_error_pr) / 3, 1.0],
-                    ],
-                ),
-            )
-            or np.allclose(
-                res.rho,
-                0.5
-                * np.array(
-                    [
-                        [1.0, np.exp(-1j * rad) * (3 - 4 * x_error_pr) / 3],
-                        [np.exp(1j * rad) * (3 - 4 * x_error_pr) / 3, 1.0],
-                    ],
-                ),
-            )
+        assert np.allclose(
+            res.rho,
+            0.5
+            * np.array(
+                [
+                    [1.0, np.exp(-1j * rad) * (3 - 4 * x_error_pr) * (3 - 4 * z_error_pr) / 9],
+                    [np.exp(1j * rad) * (3 - 4 * x_error_pr) * (3 - 4 * z_error_pr) / 9, 1.0],
+                ],
+            ),
         )
 
     # test measurement confuse outcome
-    def test_noisy_measure_confuse_rz(self, fx_rng: Generator) -> None:
+    @pytest.mark.parametrize("z_outcome,x_outcome", [(0, 0), (0, 1), (1, 0), (1, 1)])
+    def test_noisy_measure_confuse_rz(self, fx_rng: Generator, z_outcome: int, x_outcome: int) -> None:
         alpha = fx_rng.random()
         rzpattern = rzpat(alpha)
-        # probability 1 to shift both outcome
+
+        # M(0) determines Z, M(1) determines X
+        results = {}
+        cmd_count = 0
+        for cmd in rzpattern:
+            if cmd.kind.name == "M":
+                if cmd_count == 0:
+                    results[cmd.node] = z_outcome
+                elif cmd_count == 1:
+                    results[cmd.node] = x_outcome
+                cmd_count += 1
+
+        # Test with probability 1 to flip both outcomes
         res = rzpattern.simulate_pattern(
-            backend="densitymatrix", noise_model=DepolarisingNoiseModel(measure_error_prob=1.0), rng=fx_rng
+            backend="densitymatrix",
+            noise_model=DepolarisingNoiseModel(measure_error_prob=1.0),
+            branch_selector=FixedBranchSelector(results),
+            rng=fx_rng,
         )
-        # result X, XZ or Z
 
         exact = rz_exact_res(alpha)
-
         assert isinstance(res, DensityMatrix)
-        assert (
-            np.allclose(res.rho, Ops.X @ exact @ Ops.X)
-            or np.allclose(res.rho, Ops.Z @ exact @ Ops.Z)
-            or np.allclose(res.rho, Ops.Z @ Ops.X @ exact @ Ops.X @ Ops.Z)
-        )
+        # All outcomes lead to same result: both corrections applied due to flipping
+        assert np.allclose(res.rho, Ops.Z @ Ops.X @ exact @ Ops.X @ Ops.Z)
 
-        # arbitrary probability
+    @pytest.mark.parametrize("z_outcome,x_outcome", [(0, 0), (0, 1), (1, 0), (1, 1)])
+    def test_noisy_measure_confuse_rz_arbitrary(self, fx_rng: Generator, z_outcome: int, x_outcome: int) -> None:
+        alpha = fx_rng.random()
+        rzpattern = rzpat(alpha)
+
+        # M(0) determines Z, M(1) determines X
+        results = {}
+        cmd_count = 0
+        for cmd in rzpattern:
+            if cmd.kind.name == "M":
+                if cmd_count == 0:
+                    results[cmd.node] = z_outcome
+                elif cmd_count == 1:
+                    results[cmd.node] = x_outcome
+                cmd_count += 1
+
+        # Test with arbitrary probability
         measure_error_pr = fx_rng.random()
-        print(f"measure_error_pr = {measure_error_pr}")
+        print(f"measure_error_pr = {measure_error_pr}, z_outcome = {z_outcome}, x_outcome = {x_outcome}")
         res = rzpattern.simulate_pattern(
             backend="densitymatrix",
             noise_model=DepolarisingNoiseModel(measure_error_prob=measure_error_pr),
+            branch_selector=FixedBranchSelector(results),
             rng=fx_rng,
         )
-        # just add the case without readout errors
+
+        exact = rz_exact_res(alpha)
         assert isinstance(res, DensityMatrix)
+
+        # With arbitrary measure_error_pr, outcomes may or may not be flipped
+        # The physical result depends on whether the error occurs
+        # We check all possible cases
         assert (
             np.allclose(res.rho, exact)
             or np.allclose(res.rho, Ops.X @ exact @ Ops.X)

From 6d2b0f687f8c68a11d22f98dae8c20eb1c655e58 Mon Sep 17 00:00:00 2001
From: thierry-martinez <thierry.martinez@inria.fr>
Date: Fri, 6 Feb 2026 09:10:41 +0100
Subject: [PATCH 02/20] Prepare release v0.3.4 (#424)

* Update CHANGELOG

* Remove trailing spaces

* Update graphix-qasm-parser constraint
---
 CHANGELOG.md         | 26 ++++++++++++++++----------
 requirements-dev.txt |  2 +-
 2 files changed, 17 insertions(+), 11 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index c739378b2..761e8bce8 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,7 +5,7 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
-## Unreleased
+## [0.3.4] - 2026-02-05
 
 ### Added
 
@@ -20,7 +20,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - #385
   - Introduced `graphix.flow.core.XZCorrections.check_well_formed` which verifies the correctness of an XZ-corrections instance and raises an exception if incorrect.
   - Added XZ-correction exceptions to module `graphix.flow.core.exceptions`.
-  
+
 - #378:
   - Introduced new method `graphix.flow.core.PauliFlow.check_well_formed`, `graphix.flow.core.GFlow.check_well_formed` and `graphix.flow.core.CausalFlow.check_well_formed` which verify the correctness of flow objects and raise exceptions when the flow is incorrect.
   - Introduced new method `graphix.flow.core.PauliFlow.is_well_formed` which verify the correctness of flow objects and returns a boolean when the flow is incorrect.
@@ -37,13 +37,15 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - #402: Support for Python 3.14.
 
-- #407: Introduced new method `graphix.optimization.StandardizedPattern.extract_xzcorrections` and its wrapper `graphix.pattern.Pattern.extract_xzcorrections` which extract an `XZCorrections` instance from a pattern. 
+- #253, #406: Added classes `BaseCommand` and `BaseInstruction`.
+
+- #407: Introduced new method `graphix.optimization.StandardizedPattern.extract_xzcorrections` and its wrapper `graphix.pattern.Pattern.extract_xzcorrections` which extract an `XZCorrections` instance from a pattern.
 
 - #412: Added pretty-print methods (`to_ascii`, `to_latex` and `to_unicode`) for `PauliFlow` and `XZCorrections` classes. Implemented their `__str__` method as a call to `self.to_ascii`.
 
 ### Fixed
 
-- #392: `Pattern.remove_input_nodes` is required before the `Pattern.perform_pauli_measurements` method to ensure input nodes are removed and fixed in the |+> state.
+- #363, #392: `Pattern.remove_input_nodes` is required before the `Pattern.perform_pauli_measurements` method to ensure input nodes are removed and fixed in the |+> state.
 
 - #379: Removed unnecessary `meas_index` from API for rotation instructions `RZ`, `RY` and `RX`.
 
@@ -59,10 +61,14 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   is ensured with normalization passed `incorporate_pauli_results` and
   `single_qubit_domains`.
 
-- #409: Axis labels are shown when visualizing a pattern. Legend is placed outside the plot so that the graph remains visible.
+- #231, #405: `IXYZ` is now defined as `Literal[I] | Axis`.
+
+- #382, #409: Axis labels are shown when visualizing a pattern. Legend is placed outside the plot so that the graph remains visible.
 
 - #407: Fixed an unreported bug in `OpenGraph.is_equal_structurally` which failed to compare open graphs differing on the output nodes only.
 
+- #157, #417: `Pattern.minimize_space` uses `Pattern.extract_causal_flow()` and preserves runnability
+
 ### Changed
 
 - #396: Removed generic `BackendState` from `graphix.sim` modules and methods in `graphix.pattern` and `graphix.simulator` modules.
@@ -82,7 +88,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   - Modified the constructor `XZCorrections.from_measured_nodes_mapping` so that it doesn't need to create an `nx.DiGraph` instance. This fixes an unreported bug in the method.
   - Removed modules `graphix.gflow` and `graphix.find_pflow`.
 
-- #414: Tests are now type-checked.
+- #369, #414: `random_objects.py` and tests are now type-checked.
 
 - #418: `Pattern.extract_measurement_commands` now returns a dictionary. Removed `Pattern.get_meas_plane` and `Pattern.get_angles`.
 
@@ -106,7 +112,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   a pattern.
 
 - #358: Refactor of flow tools - Part I
-  - New module `graphix.flow.core` which introduces classes `PauliFlow`, `GFlow`, `CausalFlow` and `XZCorrections` allowing a finer analysis of MBQC flows. This module subsumes `graphix.generator` which has been removed and part of `graphix.gflow` which will be removed in the future. 
+  - New module `graphix.flow.core` which introduces classes `PauliFlow`, `GFlow`, `CausalFlow` and `XZCorrections` allowing a finer analysis of MBQC flows. This module subsumes `graphix.generator` which has been removed and part of `graphix.gflow` which will be removed in the future.
   - New module `graphix.flow._find_cflow` with the existing causal-flow finding algorithm.
   - New module `graphix.flow._find_gpflow` with the existing g- and Pauli-flow finding algorithm introduced in #337.
   - New abstract types `graphix.fundamentals.AbstractMeasurement` and `graphix.fundamentals.AbstractPlanarMeasurement` which serve as an umbrella of the existing types `graphix.measurements.Measurement`, `graphix.fundamentals.Plane` and `graphix.fundamentals.Axis`.
@@ -215,7 +221,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
   * a function `incorporate_pauli_results` that returns an equivalent
     pattern in which the `results` are incorporated into measurement
-    and correction domains.  
+    and correction domains.
     The resulting pattern is suitable for flow analysis. In
     particular, if a pattern has a flow, it is preserved by
     `perform_pauli_measurements` after applying `standardize` and
@@ -266,11 +272,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - #314, #322: The method `Pattern.standardize()` now places C commands
   after X and Z commands, making the resulting patterns suitable for
-  flow analysis.  
+  flow analysis.
   The `flow_from_pattern` functions now fail if the input pattern is
   not strictly standardized (as checked by
   `Pattern.is_standard(strict=True)`, which requires C commands to be
-  last).  
+  last).
   Note: the method `perform_pauli_measurements` still places C
   commands before X and Z commands.
 
diff --git a/requirements-dev.txt b/requirements-dev.txt
index 5f0fe5002..32d2d6f0c 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -26,4 +26,4 @@ qiskit_qasm3_import
 qiskit-aer; python_version < "3.14"
 
 openqasm-parser>=3.1.0
-graphix-qasm-parser @ git+https://github.com/TeamGraphix/graphix-qasm-parser.git
+graphix-qasm-parser>=0.1.1

From 0115aed95910995f649addb62e8b160790cc3ee0 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tiberiu=20Sab=C4=83u?=
 <96194994+tibisabau@users.noreply.github.com>
Date: Fri, 6 Feb 2026 09:21:19 +0100
Subject: [PATCH 03/20] Update tests/test_noisy_density_matrix.py

Co-authored-by: thierry-martinez <thierry.martinez@inria.fr>
---
 tests/test_noisy_density_matrix.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/tests/test_noisy_density_matrix.py b/tests/test_noisy_density_matrix.py
index 87e41828d..aacb51ea9 100644
--- a/tests/test_noisy_density_matrix.py
+++ b/tests/test_noisy_density_matrix.py
@@ -319,7 +319,8 @@ def test_noisy_measure_channel_rz(self, fx_rng: Generator) -> None:
             ),
         )
 
-    @pytest.mark.parametrize("outcome_z,outcome_x", [(0, 0), (0, 1), (1, 0), (1, 1)])
+    @pytest.mark.parametrize("z_outcome", [0, 1])
+    @pytest.mark.parametrize("x_outcome", [0, 1])
     def test_noisy_x_rz(self, fx_rng: Generator, outcome_z: int, outcome_x: int) -> None:
         alpha = fx_rng.random()
         rzpattern = rzpat(alpha)

From 647adac47c0ba8636a9df98db911b2e39fa53aec Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tiberiu=20Sab=C4=83u?=
 <96194994+tibisabau@users.noreply.github.com>
Date: Fri, 6 Feb 2026 09:21:42 +0100
Subject: [PATCH 04/20] Update tests/test_noisy_density_matrix.py

Co-authored-by: thierry-martinez <thierry.martinez@inria.fr>
---
 tests/test_noisy_density_matrix.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test_noisy_density_matrix.py b/tests/test_noisy_density_matrix.py
index aacb51ea9..be4671319 100644
--- a/tests/test_noisy_density_matrix.py
+++ b/tests/test_noisy_density_matrix.py
@@ -71,7 +71,7 @@ def test_noisy_measure_confuse_hadamard(self, fx_rng: Generator) -> None:
         assert np.allclose(res.rho, np.array([[0.0, 0.0], [0.0, 1.0]]))
 
     @pytest.mark.parametrize("outcome", [0, 1])
-    def test_noisy_measure_confuse_hadamard_arbitrary(self, fx_rng: Generator, outcome: int) -> None:
+    def test_noisy_measure_confuse_hadamard_arbitrary(self, fx_rng: Generator, outcome: Outcome) -> None:
         # arbitrary probability with fixed branch
         hadamardpattern = hpat()
         measure_error_pr = fx_rng.random()

From 44d1798eb434a213bd54e5ed7359385eebdd2fd1 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tiberiu=20Sab=C4=83u?=
 <96194994+tibisabau@users.noreply.github.com>
Date: Fri, 6 Feb 2026 09:22:07 +0100
Subject: [PATCH 05/20] Update tests/test_noisy_density_matrix.py

Co-authored-by: thierry-martinez <thierry.martinez@inria.fr>
---
 tests/test_noisy_density_matrix.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test_noisy_density_matrix.py b/tests/test_noisy_density_matrix.py
index be4671319..3ce9cc627 100644
--- a/tests/test_noisy_density_matrix.py
+++ b/tests/test_noisy_density_matrix.py
@@ -332,7 +332,7 @@ def test_noisy_x_rz(self, fx_rng: Generator, outcome_z: int, outcome_x: int) ->
         results = {}
         cmd_count = 0
         for cmd in rzpattern:
-            if cmd.kind.name == "M":
+            if cmd.kind == CommandKind.M:
                 if cmd_count == 0:
                     results[cmd.node] = outcome_z
                 elif cmd_count == 1:

From b95dfdcdebbadbb8527be36975422aa947b6229b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tiberiu=20Sab=C4=83u?=
 <96194994+tibisabau@users.noreply.github.com>
Date: Fri, 6 Feb 2026 09:22:24 +0100
Subject: [PATCH 06/20] Update tests/test_noisy_density_matrix.py

Co-authored-by: thierry-martinez <thierry.martinez@inria.fr>
---
 tests/test_noisy_density_matrix.py | 11 ++---------
 1 file changed, 2 insertions(+), 9 deletions(-)

diff --git a/tests/test_noisy_density_matrix.py b/tests/test_noisy_density_matrix.py
index 3ce9cc627..9bdd2d438 100644
--- a/tests/test_noisy_density_matrix.py
+++ b/tests/test_noisy_density_matrix.py
@@ -329,15 +329,8 @@ def test_noisy_x_rz(self, fx_rng: Generator, outcome_z: int, outcome_x: int) ->
         print(f"x_error_pr = {x_error_pr}, outcome_z = {outcome_z}, outcome_x = {outcome_x}")
 
         # M(0) determines Z, M(1) determines X
-        results = {}
-        cmd_count = 0
-        for cmd in rzpattern:
-            if cmd.kind == CommandKind.M:
-                if cmd_count == 0:
-                    results[cmd.node] = outcome_z
-                elif cmd_count == 1:
-                    results[cmd.node] = outcome_x
-                cmd_count += 1
+        m_nodes = (cmd.node for cmd in rzpattern if cmd.kind == CommandKind.M)
+        results = {next(m_nodes): outcome_z, next(m_nodes): outcome_x}
 
         res = rzpattern.simulate_pattern(
             backend="densitymatrix",

From f228370d2433fa94c42cd48ec25d0ef39795ebbb Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tiberiu=20Sab=C4=83u?= <T.Sabau@student.tudelft.nl>
Date: Fri, 6 Feb 2026 15:35:38 +0100
Subject: [PATCH 07/20] Fix ruff and mypy

---
 tests/test_noisy_density_matrix.py | 54 ++++++++++++++++--------------
 1 file changed, 29 insertions(+), 25 deletions(-)

diff --git a/tests/test_noisy_density_matrix.py b/tests/test_noisy_density_matrix.py
index 9bdd2d438..13f056af1 100644
--- a/tests/test_noisy_density_matrix.py
+++ b/tests/test_noisy_density_matrix.py
@@ -1,12 +1,13 @@
 from __future__ import annotations
 
-from typing import TYPE_CHECKING
+from typing import TYPE_CHECKING, cast
 
 import numpy as np
 import numpy.typing as npt
 import pytest
 
 from graphix.branch_selector import ConstBranchSelector, FixedBranchSelector
+from graphix.command import CommandKind, M
 from graphix.fundamentals import angle_to_rad
 from graphix.noise_models import DepolarisingNoiseModel
 from graphix.noise_models.noise_model import NoiselessNoiseModel
@@ -18,6 +19,7 @@
     from numpy.random import Generator
 
     from graphix.fundamentals import Angle
+    from graphix.measurements import Outcome
     from graphix.pattern import Pattern
 
 
@@ -110,7 +112,7 @@ def test_noisy_measure_channel_hadamard(self, fx_rng: Generator) -> None:
 
     # test Pauli X error
     @pytest.mark.parametrize("outcome", [0, 1])
-    def test_noisy_x_hadamard(self, fx_rng: Generator, outcome: int) -> None:
+    def test_noisy_x_hadamard(self, fx_rng: Generator, outcome: Outcome) -> None:
         hadamardpattern = hpat()
         # x error only
         x_error_pr = fx_rng.random()
@@ -321,16 +323,16 @@ def test_noisy_measure_channel_rz(self, fx_rng: Generator) -> None:
 
     @pytest.mark.parametrize("z_outcome", [0, 1])
     @pytest.mark.parametrize("x_outcome", [0, 1])
-    def test_noisy_x_rz(self, fx_rng: Generator, outcome_z: int, outcome_x: int) -> None:
+    def test_noisy_x_rz(self, fx_rng: Generator, z_outcome: Outcome, x_outcome: Outcome) -> None:
         alpha = fx_rng.random()
         rzpattern = rzpat(alpha)
         # x error only
         x_error_pr = fx_rng.random()
-        print(f"x_error_pr = {x_error_pr}, outcome_z = {outcome_z}, outcome_x = {outcome_x}")
+        print(f"x_error_pr = {x_error_pr}, outcome_z = {z_outcome}, outcome_x = {x_outcome}")
 
         # M(0) determines Z, M(1) determines X
         m_nodes = (cmd.node for cmd in rzpattern if cmd.kind == CommandKind.M)
-        results = {next(m_nodes): outcome_z, next(m_nodes): outcome_x}
+        results: dict[int, Outcome] = {next(m_nodes): z_outcome, next(m_nodes): x_outcome}
 
         res = rzpattern.simulate_pattern(
             backend="densitymatrix",
@@ -353,8 +355,8 @@ def test_noisy_x_rz(self, fx_rng: Generator, outcome_z: int, outcome_x: int) ->
             ),
         )
 
-    @pytest.mark.parametrize("outcome_z,outcome_x", [(0, 0), (0, 1), (1, 0), (1, 1)])
-    def test_noisy_z_rz(self, fx_rng: Generator, outcome_z: int, outcome_x: int) -> None:
+    @pytest.mark.parametrize(("outcome_z", "outcome_x"), [(0, 0), (0, 1), (1, 0), (1, 1)])
+    def test_noisy_z_rz(self, fx_rng: Generator, outcome_z: Outcome, outcome_x: Outcome) -> None:
         alpha = fx_rng.random()
         rzpattern = rzpat(alpha)
         # z error only
@@ -362,14 +364,14 @@ def test_noisy_z_rz(self, fx_rng: Generator, outcome_z: int, outcome_x: int) ->
         print(f"z_error_pr = {z_error_pr}, outcome_z = {outcome_z}, outcome_x = {outcome_x}")
 
         # M(0) determines Z, M(1) determines X
-        results = {}
+        results: dict[int, Outcome] = {}
         cmd_count = 0
         for cmd in rzpattern:
             if cmd.kind.name == "M":
                 if cmd_count == 0:
-                    results[cmd.node] = outcome_z
+                    results[cast("M", cmd).node] = outcome_z
                 elif cmd_count == 1:
-                    results[cmd.node] = outcome_x
+                    results[cast("M", cmd).node] = outcome_x
                 cmd_count += 1
 
         res = rzpattern.simulate_pattern(
@@ -393,8 +395,8 @@ def test_noisy_z_rz(self, fx_rng: Generator, outcome_z: int, outcome_x: int) ->
             ),
         )
 
-    @pytest.mark.parametrize("z_outcome,x_outcome", [(0, 0), (0, 1), (1, 0), (1, 1)])
-    def test_noisy_xz_rz(self, fx_rng: Generator, z_outcome: int, x_outcome: int) -> None:
+    @pytest.mark.parametrize(("z_outcome", "x_outcome"), [(0, 0), (0, 1), (1, 0), (1, 1)])
+    def test_noisy_xz_rz(self, fx_rng: Generator, z_outcome: Outcome, x_outcome: Outcome) -> None:
         alpha = fx_rng.random()
         rzpattern = rzpat(alpha)
         # x and z errors
@@ -405,14 +407,14 @@ def test_noisy_xz_rz(self, fx_rng: Generator, z_outcome: int, x_outcome: int) ->
         print(f"z_outcome = {z_outcome}, x_outcome = {x_outcome}")
 
         # M(0) determines Z correction, M(1) determines X correction
-        results = {}
+        results: dict[int, Outcome] = {}
         cmd_count = 0
         for cmd in rzpattern:
             if cmd.kind.name == "M":
                 if cmd_count == 0:
-                    results[cmd.node] = z_outcome
+                    results[cast("M", cmd).node] = z_outcome
                 elif cmd_count == 1:
-                    results[cmd.node] = x_outcome
+                    results[cast("M", cmd).node] = x_outcome
                 cmd_count += 1
 
         res = rzpattern.simulate_pattern(
@@ -437,20 +439,20 @@ def test_noisy_xz_rz(self, fx_rng: Generator, z_outcome: int, x_outcome: int) ->
         )
 
     # test measurement confuse outcome
-    @pytest.mark.parametrize("z_outcome,x_outcome", [(0, 0), (0, 1), (1, 0), (1, 1)])
-    def test_noisy_measure_confuse_rz(self, fx_rng: Generator, z_outcome: int, x_outcome: int) -> None:
+    @pytest.mark.parametrize(("z_outcome", "x_outcome"), [(0, 0), (0, 1), (1, 0), (1, 1)])
+    def test_noisy_measure_confuse_rz(self, fx_rng: Generator, z_outcome: Outcome, x_outcome: Outcome) -> None:
         alpha = fx_rng.random()
         rzpattern = rzpat(alpha)
 
         # M(0) determines Z, M(1) determines X
-        results = {}
+        results: dict[int, Outcome] = {}
         cmd_count = 0
         for cmd in rzpattern:
             if cmd.kind.name == "M":
                 if cmd_count == 0:
-                    results[cmd.node] = z_outcome
+                    results[cast("M", cmd).node] = z_outcome
                 elif cmd_count == 1:
-                    results[cmd.node] = x_outcome
+                    results[cast("M", cmd).node] = x_outcome
                 cmd_count += 1
 
         # Test with probability 1 to flip both outcomes
@@ -466,20 +468,22 @@ def test_noisy_measure_confuse_rz(self, fx_rng: Generator, z_outcome: int, x_out
         # All outcomes lead to same result: both corrections applied due to flipping
         assert np.allclose(res.rho, Ops.Z @ Ops.X @ exact @ Ops.X @ Ops.Z)
 
-    @pytest.mark.parametrize("z_outcome,x_outcome", [(0, 0), (0, 1), (1, 0), (1, 1)])
-    def test_noisy_measure_confuse_rz_arbitrary(self, fx_rng: Generator, z_outcome: int, x_outcome: int) -> None:
+    @pytest.mark.parametrize(("z_outcome", "x_outcome"), [(0, 0), (0, 1), (1, 0), (1, 1)])
+    def test_noisy_measure_confuse_rz_arbitrary(
+        self, fx_rng: Generator, z_outcome: Outcome, x_outcome: Outcome
+    ) -> None:
         alpha = fx_rng.random()
         rzpattern = rzpat(alpha)
 
         # M(0) determines Z, M(1) determines X
-        results = {}
+        results: dict[int, Outcome] = {}
         cmd_count = 0
         for cmd in rzpattern:
             if cmd.kind.name == "M":
                 if cmd_count == 0:
-                    results[cmd.node] = z_outcome
+                    results[cast("M", cmd).node] = z_outcome
                 elif cmd_count == 1:
-                    results[cmd.node] = x_outcome
+                    results[cast("M", cmd).node] = x_outcome
                 cmd_count += 1
 
         # Test with arbitrary probability

From e2f54e32a083ff7612094b816190ea99e0266e2e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tiberiu=20Sab=C4=83u?= <T.Sabau@student.tudelft.nl>
Date: Fri, 6 Feb 2026 17:24:23 +0100
Subject: [PATCH 08/20] fix: Remove loops and cast

---
 tests/test_noisy_density_matrix.py | 56 ++++++++----------------------
 1 file changed, 14 insertions(+), 42 deletions(-)

diff --git a/tests/test_noisy_density_matrix.py b/tests/test_noisy_density_matrix.py
index 13f056af1..e4d867b07 100644
--- a/tests/test_noisy_density_matrix.py
+++ b/tests/test_noisy_density_matrix.py
@@ -1,13 +1,13 @@
 from __future__ import annotations
 
-from typing import TYPE_CHECKING, cast
+from typing import TYPE_CHECKING
 
 import numpy as np
 import numpy.typing as npt
 import pytest
 
 from graphix.branch_selector import ConstBranchSelector, FixedBranchSelector
-from graphix.command import CommandKind, M
+from graphix.command import CommandKind
 from graphix.fundamentals import angle_to_rad
 from graphix.noise_models import DepolarisingNoiseModel
 from graphix.noise_models.noise_model import NoiselessNoiseModel
@@ -355,7 +355,8 @@ def test_noisy_x_rz(self, fx_rng: Generator, z_outcome: Outcome, x_outcome: Outc
             ),
         )
 
-    @pytest.mark.parametrize(("outcome_z", "outcome_x"), [(0, 0), (0, 1), (1, 0), (1, 1)])
+    @pytest.mark.parametrize("outcome_z", [0, 1])
+    @pytest.mark.parametrize("outcome_x", [0, 1])
     def test_noisy_z_rz(self, fx_rng: Generator, outcome_z: Outcome, outcome_x: Outcome) -> None:
         alpha = fx_rng.random()
         rzpattern = rzpat(alpha)
@@ -364,15 +365,7 @@ def test_noisy_z_rz(self, fx_rng: Generator, outcome_z: Outcome, outcome_x: Outc
         print(f"z_error_pr = {z_error_pr}, outcome_z = {outcome_z}, outcome_x = {outcome_x}")
 
         # M(0) determines Z, M(1) determines X
-        results: dict[int, Outcome] = {}
-        cmd_count = 0
-        for cmd in rzpattern:
-            if cmd.kind.name == "M":
-                if cmd_count == 0:
-                    results[cast("M", cmd).node] = outcome_z
-                elif cmd_count == 1:
-                    results[cast("M", cmd).node] = outcome_x
-                cmd_count += 1
+        results: dict[int, Outcome] = {0: outcome_z, 1: outcome_x}
 
         res = rzpattern.simulate_pattern(
             backend="densitymatrix",
@@ -395,7 +388,8 @@ def test_noisy_z_rz(self, fx_rng: Generator, outcome_z: Outcome, outcome_x: Outc
             ),
         )
 
-    @pytest.mark.parametrize(("z_outcome", "x_outcome"), [(0, 0), (0, 1), (1, 0), (1, 1)])
+    @pytest.mark.parametrize("z_outcome", [0, 1])
+    @pytest.mark.parametrize("x_outcome", [0, 1])
     def test_noisy_xz_rz(self, fx_rng: Generator, z_outcome: Outcome, x_outcome: Outcome) -> None:
         alpha = fx_rng.random()
         rzpattern = rzpat(alpha)
@@ -407,15 +401,7 @@ def test_noisy_xz_rz(self, fx_rng: Generator, z_outcome: Outcome, x_outcome: Out
         print(f"z_outcome = {z_outcome}, x_outcome = {x_outcome}")
 
         # M(0) determines Z correction, M(1) determines X correction
-        results: dict[int, Outcome] = {}
-        cmd_count = 0
-        for cmd in rzpattern:
-            if cmd.kind.name == "M":
-                if cmd_count == 0:
-                    results[cast("M", cmd).node] = z_outcome
-                elif cmd_count == 1:
-                    results[cast("M", cmd).node] = x_outcome
-                cmd_count += 1
+        results: dict[int, Outcome] = {0: z_outcome, 1: x_outcome}
 
         res = rzpattern.simulate_pattern(
             backend="densitymatrix",
@@ -439,21 +425,14 @@ def test_noisy_xz_rz(self, fx_rng: Generator, z_outcome: Outcome, x_outcome: Out
         )
 
     # test measurement confuse outcome
-    @pytest.mark.parametrize(("z_outcome", "x_outcome"), [(0, 0), (0, 1), (1, 0), (1, 1)])
+    @pytest.mark.parametrize("z_outcome", [0, 1])
+    @pytest.mark.parametrize("x_outcome", [0, 1])
     def test_noisy_measure_confuse_rz(self, fx_rng: Generator, z_outcome: Outcome, x_outcome: Outcome) -> None:
         alpha = fx_rng.random()
         rzpattern = rzpat(alpha)
 
         # M(0) determines Z, M(1) determines X
-        results: dict[int, Outcome] = {}
-        cmd_count = 0
-        for cmd in rzpattern:
-            if cmd.kind.name == "M":
-                if cmd_count == 0:
-                    results[cast("M", cmd).node] = z_outcome
-                elif cmd_count == 1:
-                    results[cast("M", cmd).node] = x_outcome
-                cmd_count += 1
+        results: dict[int, Outcome] = {0: z_outcome, 1: x_outcome}
 
         # Test with probability 1 to flip both outcomes
         res = rzpattern.simulate_pattern(
@@ -468,7 +447,8 @@ def test_noisy_measure_confuse_rz(self, fx_rng: Generator, z_outcome: Outcome, x
         # All outcomes lead to same result: both corrections applied due to flipping
         assert np.allclose(res.rho, Ops.Z @ Ops.X @ exact @ Ops.X @ Ops.Z)
 
-    @pytest.mark.parametrize(("z_outcome", "x_outcome"), [(0, 0), (0, 1), (1, 0), (1, 1)])
+    @pytest.mark.parametrize("z_outcome", [0, 1])
+    @pytest.mark.parametrize("x_outcome", [0, 1])
     def test_noisy_measure_confuse_rz_arbitrary(
         self, fx_rng: Generator, z_outcome: Outcome, x_outcome: Outcome
     ) -> None:
@@ -476,15 +456,7 @@ def test_noisy_measure_confuse_rz_arbitrary(
         rzpattern = rzpat(alpha)
 
         # M(0) determines Z, M(1) determines X
-        results: dict[int, Outcome] = {}
-        cmd_count = 0
-        for cmd in rzpattern:
-            if cmd.kind.name == "M":
-                if cmd_count == 0:
-                    results[cast("M", cmd).node] = z_outcome
-                elif cmd_count == 1:
-                    results[cast("M", cmd).node] = x_outcome
-                cmd_count += 1
+        results: dict[int, Outcome] = {0: z_outcome, 1: x_outcome}
 
         # Test with arbitrary probability
         measure_error_pr = fx_rng.random()

From 35871a0a7196155ca67ff7eaf1c124de549e3e0b Mon Sep 17 00:00:00 2001
From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com>
Date: Mon, 9 Feb 2026 12:39:19 +0100
Subject: [PATCH 09/20] Bump ruff from 0.14.14 to 0.15.0 in the python-packages
 group (#427)

Bumps the python-packages group with 1 update: [ruff](https://github.com/astral-sh/ruff).


Updates `ruff` from 0.14.14 to 0.15.0
- [Release notes](https://github.com/astral-sh/ruff/releases)
- [Changelog](https://github.com/astral-sh/ruff/blob/main/CHANGELOG.md)
- [Commits](https://github.com/astral-sh/ruff/compare/0.14.14...0.15.0)

---
updated-dependencies:
- dependency-name: ruff
  dependency-version: 0.15.0
  dependency-type: direct:development
  update-type: version-update:semver-minor
  dependency-group: python-packages
...

Signed-off-by: dependabot[bot] <support@github.com>
Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com>
---
 requirements-dev.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/requirements-dev.txt b/requirements-dev.txt
index 32d2d6f0c..12c6dcc20 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -2,7 +2,7 @@
 mypy==1.19.1
 pre-commit # for language-agnostic hooks
 pyright
-ruff==0.14.14
+ruff==0.15.0
 
 # Stubs
 types-networkx==3.6.1.20251220

From 643feb8e4cc7a237dbd45e903731840653a1cf2f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tiberiu=20Sab=C4=83u?=
 <96194994+tibisabau@users.noreply.github.com>
Date: Mon, 9 Feb 2026 13:15:30 +0100
Subject: [PATCH 10/20] Replace ValueError with PatternError (#426)

---
 graphix/pattern.py    | 26 +++++++++++++++-----------
 tests/test_pattern.py | 18 ++++++++++--------
 2 files changed, 25 insertions(+), 19 deletions(-)

diff --git a/graphix/pattern.py b/graphix/pattern.py
index 0f4ec4f25..f52313a6d 100644
--- a/graphix/pattern.py
+++ b/graphix/pattern.py
@@ -211,7 +211,7 @@ def compose(
         nodes_p2 = other.extract_nodes() | other.results.keys()
 
         if not mapping.keys() <= nodes_p2:
-            raise ValueError("Keys of `mapping` must correspond to the nodes of `other`.")
+            raise PatternError("Keys of `mapping` must correspond to the nodes of `other`.")
 
         # Cast to set for improved performance in membership test
         mapping_values_set = set(mapping.values())
@@ -219,14 +219,14 @@ def compose(
         i2_set = set(other.input_nodes)
 
         if len(mapping) != len(mapping_values_set):
-            raise ValueError("Values of `mapping` contain duplicates.")
+            raise PatternError("Values of `mapping` contain duplicates.")
 
         if mapping_values_set & nodes_p1 - o1_set:
-            raise ValueError("Values of `mapping` must not contain measured nodes of pattern `self`.")
+            raise PatternError("Values of `mapping` must not contain measured nodes of pattern `self`.")
 
         for k, v in mapping.items():
             if v in o1_set and k not in i2_set:
-                raise ValueError(
+                raise PatternError(
                     f"Mapping {k} -> {v} is not valid. {v} is an output of pattern `self` but {k} is not an input of pattern `other`."
                 )
 
@@ -506,7 +506,7 @@ def shift_signals(self, method: str = "direct") -> dict[int, set[int]]:
                     self._commute_with_following(target)
                 target += 1
             return signal_dict
-        raise ValueError("Invalid method")
+        raise PatternError("Invalid method")
 
     def shift_signals_direct(self) -> dict[int, set[int]]:
         """Perform signal shifting procedure."""
@@ -1153,7 +1153,7 @@ def extract_opengraph(self) -> OpenGraph[Measurement]:
         for cmd in self.__seq:
             if cmd.kind == CommandKind.N:
                 if cmd.state != BasicStates.PLUS:
-                    raise ValueError(
+                    raise PatternError(
                         f"Open graph extraction requires N commands to represent a |+⟩ state. Error found in {cmd}."
                     )
                 nodes.add(cmd.node)
@@ -1411,7 +1411,7 @@ def perform_pauli_measurements(self, ignore_pauli_with_deps: bool = False) -> No
 
         """
         if self.input_nodes:
-            raise ValueError("Remove inputs with `self.remove_input_nodes()` before performing Pauli presimulation.")
+            raise PatternError("Remove inputs with `self.remove_input_nodes()` before performing Pauli presimulation.")
         self.__dict__.update(measure_pauli(self, ignore_pauli_with_deps=ignore_pauli_with_deps).__dict__)
 
     def draw_graph(
@@ -1596,6 +1596,10 @@ def check_measured(cmd: Command, node: int) -> None:
                 check_active(cmd, cmd.node)
 
 
+class PatternError(Exception):
+    """Exception subclass to handle pattern errors."""
+
+
 class RunnabilityErrorReason(Enum):
     """Describe the reason for a pattern not being runnable."""
 
@@ -1616,7 +1620,7 @@ class RunnabilityErrorReason(Enum):
 
 
 @dataclass
-class RunnabilityError(Exception):
+class RunnabilityError(PatternError):
     """Error raised by :method:`Pattern.check_runnability`."""
 
     cmd: Command
@@ -1778,7 +1782,7 @@ def pauli_nodes(pattern: optimization.StandardizedPattern) -> tuple[list[tuple[c
                 else:
                     pauli_node.append((cmd, pm))
             else:
-                raise ValueError("Unknown Pauli measurement basis")
+                raise PatternError("Unknown Pauli measurement basis")
         else:
             non_pauli_node.add(cmd.node)
     return pauli_node, non_pauli_node
@@ -1788,12 +1792,12 @@ def assert_permutation(original: list[int], user: list[int]) -> None:
     """Check that the provided `user` node list is a permutation from `original`."""
     node_set = set(user)
     if node_set != set(original):
-        raise ValueError(f"{node_set} != {set(original)}")
+        raise PatternError(f"{node_set} != {set(original)}")
     for node in user:
         if node in node_set:
             node_set.remove(node)
         else:
-            raise ValueError(f"{node} appears twice")
+            raise PatternError(f"{node} appears twice")
 
 
 @dataclass
diff --git a/tests/test_pattern.py b/tests/test_pattern.py
index ff6306c9b..b0bd65d35 100644
--- a/tests/test_pattern.py
+++ b/tests/test_pattern.py
@@ -19,7 +19,7 @@
 from graphix.fundamentals import ANGLE_PI, Angle, Plane
 from graphix.measurements import Measurement, Outcome, PauliMeasurement
 from graphix.opengraph import OpenGraph
-from graphix.pattern import Pattern, RunnabilityError, RunnabilityErrorReason, shift_outcomes
+from graphix.pattern import Pattern, PatternError, RunnabilityError, RunnabilityErrorReason, shift_outcomes
 from graphix.random_objects import rand_circuit, rand_gate
 from graphix.sim.density_matrix import DensityMatrix
 from graphix.sim.statevec import Statevec
@@ -53,7 +53,7 @@ def test_init(self) -> None:
         pattern = Pattern(input_nodes=[1, 0], cmds=[N(node=2), M(node=1)], output_nodes=[2, 0])
         assert pattern.input_nodes == [1, 0]
         assert pattern.output_nodes == [2, 0]
-        with pytest.raises(ValueError):
+        with pytest.raises(PatternError):
             Pattern(input_nodes=[1, 0], cmds=[N(node=2), M(node=1)], output_nodes=[0, 1, 2])
 
     def test_eq(self) -> None:
@@ -300,7 +300,7 @@ def test_pauli_measurement_error(self, fx_rng: Generator) -> None:
         circuit = rand_circuit(nqubits, depth, fx_rng)
         pattern = circuit.transpile().pattern
         pattern.standardize()
-        with pytest.raises(ValueError):
+        with pytest.raises(PatternError):
             pattern.perform_pauli_measurements()
 
     def test_pauli_measurement_leave_input(self) -> None:
@@ -321,7 +321,7 @@ def test_pauli_measurement_leave_input(self) -> None:
         swap(circuit, 0, 2)
         pattern = circuit.transpile().pattern
         pattern.standardize()
-        with pytest.raises(ValueError):
+        with pytest.raises(PatternError):
             pattern.perform_pauli_measurements()
 
     @pytest.mark.parametrize("jumps", range(1, 6))
@@ -527,17 +527,19 @@ def test_compose_1(self) -> None:
         assert pc == p
         assert mapping_c == {0: 1, 2: 5}
 
-        with pytest.raises(ValueError, match=r"Keys of `mapping` must correspond to the nodes of `other`."):
+        with pytest.raises(PatternError, match=r"Keys of `mapping` must correspond to the nodes of `other`."):
             p1.compose(p2, mapping={0: 1, 2: 5, 1: 2})
 
-        with pytest.raises(ValueError, match=r"Values of `mapping` contain duplicates."):
+        with pytest.raises(PatternError, match=r"Values of `mapping` contain duplicates."):
             p1.compose(p2, mapping={0: 1, 2: 1})
 
-        with pytest.raises(ValueError, match=r"Values of `mapping` must not contain measured nodes of pattern `self`."):
+        with pytest.raises(
+            PatternError, match=r"Values of `mapping` must not contain measured nodes of pattern `self`."
+        ):
             p1.compose(p2, mapping={0: 1, 2: 0})
 
         with pytest.raises(
-            ValueError,
+            PatternError,
             match=r"Mapping 2 -> 1 is not valid. 1 is an output of pattern `self` but 2 is not an input of pattern `other`.",
         ):
             p1.compose(p2, mapping={2: 1})

From f05e6cd62a8eafea3be00e77b3b1e603ccfa3d96 Mon Sep 17 00:00:00 2001
From: thierry-martinez <thierry.martinez@inria.fr>
Date: Thu, 12 Feb 2026 10:51:04 +0100
Subject: [PATCH 11/20] Force test reverse dependencies with the current
 graphix codebase (#431)

This commit installs graphix after the reverse dependency to ensure
that the test uses the current graphix codebase even if another
version is pinned in the reverse dependency.
---
 noxfile.py | 33 +++++++++++++++++++++------------
 1 file changed, 21 insertions(+), 12 deletions(-)

diff --git a/noxfile.py b/noxfile.py
index 07dc6885c..c0de78d54 100644
--- a/noxfile.py
+++ b/noxfile.py
@@ -142,21 +142,30 @@ def tests_reverse_dependencies(session: Session, package: ReverseDependency) ->
             f"{dirname} only supports Python versions {package.version_constraint}; current Python version: {session.python}"
         )
 
-    session.install(".")
     install_pytest(session)
     if package.doctest_modules:
         session.install("nox")
-    # Use `session.cd` as a context manager to ensure that the
-    # working directory is restored afterward. This is important
-    # because Windows cannot delete a temporary directory while it
-    # is the working directory.
-    with TemporaryDirectory() as tmpdir, session.cd(tmpdir):
-        if package.branch is None:
-            session.run("git", "clone", package.repository)
-        else:
-            session.run("git", "clone", "-b", package.branch, package.repository)
-        with session.cd(dirname):
-            session.install(".")
+    with TemporaryDirectory() as tmpdir:
+        with session.cd(tmpdir):
+            if package.branch is None:
+                session.run("git", "clone", package.repository)
+            else:
+                session.run("git", "clone", "-b", package.branch, package.repository)
+            with session.cd(dirname):
+                session.install(".")
+        # Note that `session.cd` is used as a context manager above,
+        # so that the working directory is restored at this point.  We
+        # install now the graphix package from the working directory.
+        # This is done after having installed the reverse dependency,
+        # so that we run the test with the current graphix codebase,
+        # even if another graphix version has been pinned in the
+        # reverse dependendy.
+        session.install(".")
+        # Use `session.cd` as a context manager again to ensure that the
+        # working directory is restored afterward. This is important
+        # because Windows cannot delete a temporary directory while it
+        # is the working directory.
+        with session.cd(tmpdir), session.cd(dirname):
             if package.initialization is not None:
                 package.initialization(session)
             run_pytest(session, doctest_modules=package.doctest_modules)

From be598e2c8cf46261991f71c25accea0fa238e242 Mon Sep 17 00:00:00 2001
From: mgarnier59 <50111289+mgarnier59@users.noreply.github.com>
Date: Thu, 12 Feb 2026 16:13:46 +0100
Subject: [PATCH 12/20] Fix #428 : fix conditionality of `ApplyNoise` commands
 (#429)

* fix noise conditionality and extract it

* finalizing refactor and removing test

* lint

* implement Thierry's comments

* move default `check_domain` to `MeasureMethod`

* Fix `ApplyNoise` docstring

* Update noise_model.py

* Update noise_model.py

* double backticks for doc

* Update CHANGELOG.md
---
 CHANGELOG.md                         | 19 +++++++++++++--
 graphix/noise_models/depolarising.py |  4 ++--
 graphix/noise_models/noise_model.py  | 19 ++++++++++++++-
 graphix/sim/base_backend.py          | 35 ++++++++++++++--------------
 graphix/sim/tensornet.py             | 18 +++-----------
 graphix/simulator.py                 | 19 +++++++++++++--
 tests/test_density_matrix.py         | 32 +------------------------
 7 files changed, 76 insertions(+), 70 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 761e8bce8..561f0a839 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,6 +5,19 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [Unreleased]
+
+### Added
+
+### Fixed
+
+- #429
+  - Modify `graphix.noise_models.noise_model.ApplyNoise` to handle conditionality based on a `domain` attribute (like `command.X` and `command.Z`).
+  - Moved the conditional logic to `graphix.simulator` to remove code duplication in the backends.
+  - Solves [#428](https://github.com/TeamGraphix/graphix/issues/428).
+
+### Changed
+
 ## [0.3.4] - 2026-02-05
 
 ### Added
@@ -45,6 +58,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Fixed
 
+-
+
 - #363, #392: `Pattern.remove_input_nodes` is required before the `Pattern.perform_pauli_measurements` method to ensure input nodes are removed and fixed in the |+> state.
 
 - #379: Removed unnecessary `meas_index` from API for rotation instructions `RZ`, `RY` and `RX`.
@@ -216,10 +231,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - #322: Added a new `optimization` module containing:
 
-  * a functional version of `standardize` that returns a standardized
+  - a functional version of `standardize` that returns a standardized
     pattern as a new object;
 
-  * a function `incorporate_pauli_results` that returns an equivalent
+  - a function `incorporate_pauli_results` that returns an equivalent
     pattern in which the `results` are incorporated into measurement
     and correction domains.
     The resulting pattern is suitable for flow analysis. In
diff --git a/graphix/noise_models/depolarising.py b/graphix/noise_models/depolarising.py
index 84af01815..ff72a416a 100644
--- a/graphix/noise_models/depolarising.py
+++ b/graphix/noise_models/depolarising.py
@@ -118,9 +118,9 @@ def command(self, cmd: CommandOrNoise, rng: Generator | None = None) -> list[Com
         if cmd.kind == CommandKind.M:
             return [ApplyNoise(noise=DepolarisingNoise(self.measure_channel_prob), nodes=[cmd.node]), cmd]
         if cmd.kind == CommandKind.X:
-            return [cmd, ApplyNoise(noise=DepolarisingNoise(self.x_error_prob), nodes=[cmd.node])]
+            return [cmd, ApplyNoise(noise=DepolarisingNoise(self.x_error_prob), nodes=[cmd.node], domain=cmd.domain)]
         if cmd.kind == CommandKind.Z:
-            return [cmd, ApplyNoise(noise=DepolarisingNoise(self.z_error_prob), nodes=[cmd.node])]
+            return [cmd, ApplyNoise(noise=DepolarisingNoise(self.z_error_prob), nodes=[cmd.node], domain=cmd.domain)]
         # Use of `==` here for mypy
         if cmd.kind == CommandKind.C or cmd.kind == CommandKind.T or cmd.kind == CommandKind.ApplyNoise:  # noqa: PLR1714
             return [cmd]
diff --git a/graphix/noise_models/noise_model.py b/graphix/noise_models/noise_model.py
index 9508328db..4c17d8c59 100644
--- a/graphix/noise_models/noise_model.py
+++ b/graphix/noise_models/noise_model.py
@@ -42,11 +42,28 @@ def to_kraus_channel(self) -> KrausChannel:
 
 @dataclass
 class ApplyNoise(_KindChecker):
-    """Apply noise command."""
+    """Apply noise command.
+
+    Parameters
+    ----------
+    noise : Noise
+        noise to be applied
+
+    nodes : list[Node]
+        list of node indices on which to apply noise
+
+    domain: set[Node] | None = None
+        Optional domain for conditional noise.
+        If ``None``, the noise is applied unconditionally.
+        Otherwise, the noise is applied if there is an odd number of nodes among ``domain`` that have been measured with outcome 1 (as for ``X`` and ``Z`` commands).
+        Note that the noise is never applied if ``domain`` is the empty set.
+
+    """
 
     kind: ClassVar[Literal[CommandKind.ApplyNoise]] = dataclasses.field(default=CommandKind.ApplyNoise, init=False)
     noise: Noise
     nodes: list[Node]
+    domain: set[Node] | None = None
 
 
 CommandOrNoise = Command | ApplyNoise
diff --git a/graphix/sim/base_backend.py b/graphix/sim/base_backend.py
index 62934d007..7c77750a1 100644
--- a/graphix/sim/base_backend.py
+++ b/graphix/sim/base_backend.py
@@ -29,10 +29,9 @@
 
     from graphix import command
     from graphix.measurements import Measurement, Outcome
-    from graphix.noise_models.noise_model import Noise
+    from graphix.noise_models.noise_model import ApplyNoise, Noise
     from graphix.parameter import ExpressionOrComplex, ExpressionOrFloat
     from graphix.sim.data import Data
-    from graphix.simulator import MeasureMethod
 
 
 Matrix: TypeAlias = npt.NDArray[np.object_ | np.complex128]
@@ -619,7 +618,7 @@ def add_nodes(self, nodes: Sequence[int], data: Data = BasicStates.PLUS) -> None
         Previously existing nodes remain unchanged.
         """
 
-    def apply_noise(self, nodes: Sequence[int], noise: Noise) -> None:  # noqa: ARG002,PLR6301
+    def apply_noise(self, cmd: ApplyNoise) -> None:  # noqa: ARG002,PLR6301
         """Apply noise.
 
         The default implementation of this method raises
@@ -628,6 +627,8 @@ def apply_noise(self, nodes: Sequence[int], noise: Noise) -> None:  # noqa: ARG0
         `DensityMatrixBackend`) override this method to implement
         the effect of noise.
 
+        Note: the simulator is responsible for checking that the measurement outcomes match the domain condition before calling this method.
+
         Parameters
         ----------
         nodes : sequence of ints.
@@ -642,8 +643,11 @@ def apply_clifford(self, node: int, clifford: Clifford) -> None:
         """Apply single-qubit Clifford gate, specified by vop index specified in graphix.clifford.CLIFFORD."""
 
     @abstractmethod
-    def correct_byproduct(self, cmd: command.X | command.Z, measure_method: MeasureMethod) -> None:
-        """Byproduct correction correct for the X or Z byproduct operators, by applying the X or Z gate."""
+    def correct_byproduct(self, cmd: command.X | command.Z) -> None:
+        """Byproduct correction correct for the X or Z byproduct operators, by applying the X or Z gate.
+
+        Note: the simulator is responsible for checking that the measurement outcomes match the domain condition before calling this method.
+        """
 
     @abstractmethod
     def entangle_nodes(self, edge: tuple[int, int]) -> None:
@@ -782,25 +786,22 @@ def f_expectation0() -> float:
         return outcome
 
     @override
-    def correct_byproduct(self, cmd: command.X | command.Z, measure_method: MeasureMethod) -> None:
+    def correct_byproduct(self, cmd: command.X | command.Z) -> None:
         """Byproduct correction correct for the X or Z byproduct operators, by applying the X or Z gate."""
-        if np.mod(sum(measure_method.measurement_outcome(j) for j in cmd.domain), 2) == 1:
-            op = Ops.X if cmd.kind == CommandKind.X else Ops.Z
-            self.apply_single(node=cmd.node, op=op)
+        op = Ops.X if cmd.kind == CommandKind.X else Ops.Z
+        self.apply_single(node=cmd.node, op=op)
 
     @override
-    def apply_noise(self, nodes: Sequence[int], noise: Noise) -> None:
-        """Apply noise.
+    def apply_noise(self, cmd: ApplyNoise) -> None:
+        """Apply noise for the command `:class: graphix.noise_model.ApplyNoise`.
 
         Parameters
         ----------
-        nodes : sequence of ints.
-            Target qubits
-        noise : Noise
-            Noise to apply
+        cmd : ApplyNoise
+            command ApplyNoise
         """
-        indices = [self.node_index.index(i) for i in nodes]
-        self.state.apply_noise(indices, noise)
+        indices = [self.node_index.index(i) for i in cmd.nodes]
+        self.state.apply_noise(indices, cmd.noise)
 
     def apply_single(self, node: int, op: Matrix) -> None:
         """Apply a single gate to the state."""
diff --git a/graphix/sim/tensornet.py b/graphix/sim/tensornet.py
index 08e4720ae..e86f5913c 100644
--- a/graphix/sim/tensornet.py
+++ b/graphix/sim/tensornet.py
@@ -34,7 +34,6 @@
     from graphix.clifford import Clifford
     from graphix.measurements import Measurement, Outcome
     from graphix.sim import Data
-    from graphix.simulator import MeasureMethod
 
 PrepareState: TypeAlias = str | npt.NDArray[np.complex128]
 
@@ -768,20 +767,9 @@ def measure(self, node: int, measurement: Measurement, rng: Generator | None = N
         return result
 
     @override
-    def correct_byproduct(self, cmd: command.X | command.Z, measure_method: MeasureMethod) -> None:
-        """Perform byproduct correction.
-
-        Parameters
-        ----------
-        cmd : list
-            Byproduct command
-            i.e. ['X' or 'Z', node, signal_domain]
-        measure_method : MeasureMethod
-            The measure method to use
-        """
-        if sum(measure_method.measurement_outcome(j) for j in cmd.domain) % 2 == 1:
-            op = Ops.X if isinstance(cmd, command.X) else Ops.Z
-            self.state.evolve_single(cmd.node, op, str(cmd.kind))
+    def correct_byproduct(self, cmd: command.X | command.Z) -> None:
+        op = Ops.X if isinstance(cmd, command.X) else Ops.Z
+        self.state.evolve_single(cmd.node, op, str(cmd.kind))
 
     @override
     def apply_clifford(self, node: int, clifford: Clifford) -> None:
diff --git a/graphix/simulator.py b/graphix/simulator.py
index f9a2c7906..380b12b5a 100644
--- a/graphix/simulator.py
+++ b/graphix/simulator.py
@@ -136,6 +136,16 @@ def store_measurement_outcome(self, node: int, result: Outcome) -> None:
         """
         ...
 
+    def check_domain(self, domain: Iterable[int]) -> bool:
+        """Check that the measurement outcomes match the domain condition.
+
+        Parameters
+        ----------
+        domain : Iterable[int]
+            domain on which to compute the condition for applying conditional commands.
+        """
+        return sum(self.measurement_outcome(j) for j in domain) % 2 == 1
+
 
 class DefaultMeasureMethod(MeasureMethod):
     """Default measurement method implementing standard measurement plane/angle update for MBQC."""
@@ -160,6 +170,7 @@ def __init__(self, results: Mapping[int, Outcome] | None = None):
         # results is coerced into dict, since `store_measurement_outcome` mutates it.
         self.results = {} if results is None else dict(results)
 
+    @override
     def describe_measurement(self, cmd: BaseM) -> Measurement:
         """Return the description of the measurement performed by ``cmd``.
 
@@ -181,6 +192,7 @@ def describe_measurement(self, cmd: BaseM) -> Measurement:
         angle = cmd.angle * measure_update.coeff + measure_update.add_term
         return Measurement(angle, measure_update.new_plane)
 
+    @override
     def measurement_outcome(self, node: int) -> Outcome:
         """Return the result of a previous measurement.
 
@@ -196,6 +208,7 @@ def measurement_outcome(self, node: int) -> Outcome:
         """
         return self.results[node]
 
+    @override
     def store_measurement_outcome(self, node: int, result: Outcome) -> None:
         """Store the result of a previous measurement.
 
@@ -311,7 +324,8 @@ def run(self, input_state: Data = BasicStates.PLUS, rng: Generator | None = None
                 self.__measure_method.measure(self.backend, cmd, noise_model=self.noise_model, rng=rng)
             # Use of `==` here for mypy
             elif cmd.kind == CommandKind.X or cmd.kind == CommandKind.Z:  # noqa: PLR1714
-                self.backend.correct_byproduct(cmd, self.__measure_method)
+                if self.__measure_method.check_domain(cmd.domain):
+                    self.backend.correct_byproduct(cmd)
             elif cmd.kind == CommandKind.C:
                 self.backend.apply_clifford(cmd.node, cmd.clifford)
             elif cmd.kind == CommandKind.T:
@@ -321,7 +335,8 @@ def run(self, input_state: Data = BasicStates.PLUS, rng: Generator | None = None
                 # handling of ticks during noise transpilation.
                 pass
             elif cmd.kind == CommandKind.ApplyNoise:
-                self.backend.apply_noise(cmd.nodes, cmd.noise)
+                if cmd.domain is None or self.__measure_method.check_domain(cmd.domain):
+                    self.backend.apply_noise(cmd)
             elif cmd.kind == CommandKind.S:
                 raise ValueError("S commands unexpected in simulated patterns.")
             else:
diff --git a/tests/test_density_matrix.py b/tests/test_density_matrix.py
index 3b8ab337d..7d7492584 100644
--- a/tests/test_density_matrix.py
+++ b/tests/test_density_matrix.py
@@ -16,7 +16,7 @@
 from graphix.fundamentals import ANGLE_PI, Plane
 from graphix.ops import Ops
 from graphix.sim.density_matrix import DensityMatrix, DensityMatrixBackend
-from graphix.sim.statevec import CNOT_TENSOR, CZ_TENSOR, SWAP_TENSOR, Statevec, StatevectorBackend
+from graphix.sim.statevec import CNOT_TENSOR, CZ_TENSOR, SWAP_TENSOR, Statevec
 from graphix.simulator import DefaultMeasureMethod
 from graphix.states import BasicStates, PlanarState
 from graphix.transpiler import Circuit
@@ -926,33 +926,3 @@ def test_measure(self, pr_calc: bool) -> None:
         expected_matrix_1 = np.kron(np.array([[1, 0], [0, 0]]), np.ones((2, 2)) / 2)
         expected_matrix_2 = np.kron(np.array([[0, 0], [0, 1]]), np.array([[0.5, -0.5], [-0.5, 0.5]]))
         assert np.allclose(backend.state.rho, expected_matrix_1) or np.allclose(backend.state.rho, expected_matrix_2)
-
-    def test_correct_byproduct(self) -> None:
-        measure_method = DefaultMeasureMethod()
-        dm_backend = DensityMatrixBackend()
-        dm_backend.add_nodes([0])
-        # node 0 initialized in Backend
-        dm_backend.add_nodes([1, 2])
-        dm_backend.entangle_nodes((0, 1))
-        dm_backend.entangle_nodes((1, 2))
-        measure_method.measure(dm_backend, command.M(0))
-        measure_method.measure(dm_backend, command.M(1, angle=-ANGLE_PI / 2, s_domain={0}))
-        dm_backend.correct_byproduct(command.X(2, {1}), measure_method)
-        dm_backend.correct_byproduct(command.Z(2, {0}), measure_method)
-        rho = dm_backend.state.rho
-
-        sv_backend = StatevectorBackend()
-        sv_backend.add_nodes([0])
-        # node 0 initialized in Backend
-        sv_backend.add_nodes([1, 2])
-        sv_backend.entangle_nodes((0, 1))
-        sv_backend.entangle_nodes((1, 2))
-        measure_method.measure(sv_backend, command.M(0))
-        measure_method.measure(sv_backend, command.M(1, angle=-ANGLE_PI / 2, s_domain={0}))
-        sv_backend.correct_byproduct(command.X(2, {1}), measure_method)
-        sv_backend.correct_byproduct(command.Z(2, {0}), measure_method)
-        psi = sv_backend.state.psi
-
-        assert np.allclose(
-            rho, np.outer(psi.astype(np.complex128, copy=False), psi.conj().astype(np.complex128, copy=False))
-        )

From 2ca61818eab8c8ddf7f4d9f92e3bee3c967f7cdd Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tiberiu=20Sab=C4=83u?= <T.Sabau@student.tudelft.nl>
Date: Thu, 5 Feb 2026 17:42:55 +0100
Subject: [PATCH 13/20] Make tests use the BranchSelector feature

---
 tests/test_noisy_density_matrix.py | 196 ++++++++++++++++++-----------
 1 file changed, 125 insertions(+), 71 deletions(-)

diff --git a/tests/test_noisy_density_matrix.py b/tests/test_noisy_density_matrix.py
index 20c5cf91e..87e41828d 100644
--- a/tests/test_noisy_density_matrix.py
+++ b/tests/test_noisy_density_matrix.py
@@ -6,6 +6,7 @@
 import numpy.typing as npt
 import pytest
 
+from graphix.branch_selector import ConstBranchSelector, FixedBranchSelector
 from graphix.fundamentals import angle_to_rad
 from graphix.noise_models import DepolarisingNoiseModel
 from graphix.noise_models.noise_model import NoiselessNoiseModel
@@ -69,14 +70,22 @@ def test_noisy_measure_confuse_hadamard(self, fx_rng: Generator) -> None:
         assert isinstance(res, DensityMatrix)
         assert np.allclose(res.rho, np.array([[0.0, 0.0], [0.0, 1.0]]))
 
-        # arbitrary probability
+    @pytest.mark.parametrize("outcome", [0, 1])
+    def test_noisy_measure_confuse_hadamard_arbitrary(self, fx_rng: Generator, outcome: int) -> None:
+        # arbitrary probability with fixed branch
+        hadamardpattern = hpat()
         measure_error_pr = fx_rng.random()
-        print(f"measure_error_pr = {measure_error_pr}")
+        print(f"measure_error_pr = {measure_error_pr}, outcome = {outcome}")
         res = hadamardpattern.simulate_pattern(
-            backend="densitymatrix", noise_model=DepolarisingNoiseModel(measure_error_prob=measure_error_pr), rng=fx_rng
+            backend="densitymatrix",
+            noise_model=DepolarisingNoiseModel(measure_error_prob=measure_error_pr),
+            branch_selector=ConstBranchSelector(outcome),
+            rng=fx_rng,
         )
-        # result should be |1>
         assert isinstance(res, DensityMatrix)
+        # With measure_error_prob, the outcome might be flipped, resulting in different X corrections
+        # However, we cannot predict the exact result without knowing if the error occurred
+        # So we check both possibilities
         assert np.allclose(res.rho, np.array([[1.0, 0.0], [0.0, 0.0]])) or np.allclose(
             res.rho,
             np.array([[0.0, 0.0], [0.0, 1.0]]),
@@ -100,21 +109,21 @@ def test_noisy_measure_channel_hadamard(self, fx_rng: Generator) -> None:
         )
 
     # test Pauli X error
-    def test_noisy_x_hadamard(self, fx_rng: Generator) -> None:
+    @pytest.mark.parametrize("outcome", [0, 1])
+    def test_noisy_x_hadamard(self, fx_rng: Generator, outcome: int) -> None:
         hadamardpattern = hpat()
         # x error only
         x_error_pr = fx_rng.random()
-        print(f"x_error_pr = {x_error_pr}")
+        print(f"x_error_pr = {x_error_pr}, outcome = {outcome}")
         res = hadamardpattern.simulate_pattern(
             backend="densitymatrix",
             noise_model=DepolarisingNoiseModel(x_error_prob=x_error_pr),
+            branch_selector=ConstBranchSelector(outcome),
             rng=fx_rng,
         )
-        # analytical result since deterministic pattern output is |0>.
-        # if no X applied, no noise. If X applied X noise on |0><0|
-
+        # Both outcomes lead to |0> after correction. ApplyNoise is unconditional, so noise is always applied.
         assert isinstance(res, DensityMatrix)
-        assert np.allclose(res.rho, np.array([[1.0, 0.0], [0.0, 0.0]])) or np.allclose(
+        assert np.allclose(
             res.rho,
             np.array([[1 - 2 * x_error_pr / 3.0, 0.0], [0.0, 2 * x_error_pr / 3.0]]),
         )
@@ -310,26 +319,36 @@ def test_noisy_measure_channel_rz(self, fx_rng: Generator) -> None:
             ),
         )
 
-    def test_noisy_x_rz(self, fx_rng: Generator) -> None:
+    @pytest.mark.parametrize("outcome_z,outcome_x", [(0, 0), (0, 1), (1, 0), (1, 1)])
+    def test_noisy_x_rz(self, fx_rng: Generator, outcome_z: int, outcome_x: int) -> None:
         alpha = fx_rng.random()
         rzpattern = rzpat(alpha)
         # x error only
         x_error_pr = fx_rng.random()
-        print(f"x_error_pr = {x_error_pr}")
+        print(f"x_error_pr = {x_error_pr}, outcome_z = {outcome_z}, outcome_x = {outcome_x}")
+
+        # M(0) determines Z, M(1) determines X
+        results = {}
+        cmd_count = 0
+        for cmd in rzpattern:
+            if cmd.kind.name == "M":
+                if cmd_count == 0:
+                    results[cmd.node] = outcome_z
+                elif cmd_count == 1:
+                    results[cmd.node] = outcome_x
+                cmd_count += 1
+
         res = rzpattern.simulate_pattern(
             backend="densitymatrix",
             noise_model=DepolarisingNoiseModel(x_error_prob=x_error_pr),
+            branch_selector=FixedBranchSelector(results),
             rng=fx_rng,
         )
 
-        # only two cases: if no X correction, Z or no Z correction but exact result.
-        # If X correction the noise result is the same with or without the PERFECT Z correction.
+        # All outcomes lead to the same state after corrections. ApplyNoise is unconditional, so noise is always applied.
         assert isinstance(res, DensityMatrix)
         rad = angle_to_rad(alpha)
         assert np.allclose(
-            res.rho,
-            0.5 * np.array([[1.0, np.exp(-1j * rad)], [np.exp(1j * rad), 1.0]]),
-        ) or np.allclose(
             res.rho,
             0.5
             * np.array(
@@ -340,26 +359,36 @@ def test_noisy_x_rz(self, fx_rng: Generator) -> None:
             ),
         )
 
-    def test_noisy_z_rz(self, fx_rng: Generator) -> None:
+    @pytest.mark.parametrize("outcome_z,outcome_x", [(0, 0), (0, 1), (1, 0), (1, 1)])
+    def test_noisy_z_rz(self, fx_rng: Generator, outcome_z: int, outcome_x: int) -> None:
         alpha = fx_rng.random()
         rzpattern = rzpat(alpha)
         # z error only
         z_error_pr = fx_rng.random()
-        print(f"z_error_pr = {z_error_pr}")
+        print(f"z_error_pr = {z_error_pr}, outcome_z = {outcome_z}, outcome_x = {outcome_x}")
+
+        # M(0) determines Z, M(1) determines X
+        results = {}
+        cmd_count = 0
+        for cmd in rzpattern:
+            if cmd.kind.name == "M":
+                if cmd_count == 0:
+                    results[cmd.node] = outcome_z
+                elif cmd_count == 1:
+                    results[cmd.node] = outcome_x
+                cmd_count += 1
+
         res = rzpattern.simulate_pattern(
             backend="densitymatrix",
             noise_model=DepolarisingNoiseModel(z_error_prob=z_error_pr),
+            branch_selector=FixedBranchSelector(results),
             rng=fx_rng,
         )
 
-        # only two cases: if no Z correction, X or no X correction but exact result.
-        # If Z correction the noise result is the same with or without the PERFECT X correction.
+        # All outcomes lead to the same state after corrections. ApplyNoise is unconditional, so noise is always applied.
         assert isinstance(res, DensityMatrix)
         rad = angle_to_rad(alpha)
         assert np.allclose(
-            res.rho,
-            0.5 * np.array([[1.0, np.exp(-1j * rad)], [np.exp(1j * rad), 1.0]]),
-        ) or np.allclose(
             res.rho,
             0.5
             * np.array(
@@ -370,7 +399,8 @@ def test_noisy_z_rz(self, fx_rng: Generator) -> None:
             ),
         )
 
-    def test_noisy_xz_rz(self, fx_rng: Generator) -> None:
+    @pytest.mark.parametrize("z_outcome,x_outcome", [(0, 0), (0, 1), (1, 0), (1, 1)])
+    def test_noisy_xz_rz(self, fx_rng: Generator, z_outcome: int, x_outcome: int) -> None:
         alpha = fx_rng.random()
         rzpattern = rzpat(alpha)
         # x and z errors
@@ -378,78 +408,102 @@ def test_noisy_xz_rz(self, fx_rng: Generator) -> None:
         print(f"x_error_pr = {x_error_pr}")
         z_error_pr = fx_rng.random()
         print(f"z_error_pr = {z_error_pr}")
+        print(f"z_outcome = {z_outcome}, x_outcome = {x_outcome}")
+
+        # M(0) determines Z correction, M(1) determines X correction
+        results = {}
+        cmd_count = 0
+        for cmd in rzpattern:
+            if cmd.kind.name == "M":
+                if cmd_count == 0:
+                    results[cmd.node] = z_outcome
+                elif cmd_count == 1:
+                    results[cmd.node] = x_outcome
+                cmd_count += 1
+
         res = rzpattern.simulate_pattern(
             backend="densitymatrix",
             noise_model=DepolarisingNoiseModel(x_error_prob=x_error_pr, z_error_prob=z_error_pr),
+            branch_selector=FixedBranchSelector(results),
             rng=fx_rng,
         )
 
-        # 4 cases : no corr, noisy X, noisy Z, noisy XZ.
+        # All outcomes lead to same state after corrections. Both X and Z noise applied unconditionally.
         assert isinstance(res, DensityMatrix)
         rad = angle_to_rad(alpha)
-        assert (
-            np.allclose(res.rho, 0.5 * np.array([[1.0, np.exp(-1j * rad)], [np.exp(1j * rad), 1.0]]))
-            or np.allclose(
-                res.rho,
-                0.5
-                * np.array(
-                    [
-                        [1.0, np.exp(-1j * rad) * (3 - 4 * x_error_pr) * (3 - 4 * z_error_pr) / 9],
-                        [np.exp(1j * rad) * (3 - 4 * x_error_pr) * (3 - 4 * z_error_pr) / 9, 1.0],
-                    ],
-                ),
-            )
-            or np.allclose(
-                res.rho,
-                0.5
-                * np.array(
-                    [
-                        [1.0, np.exp(-1j * rad) * (3 - 4 * z_error_pr) / 3],
-                        [np.exp(1j * rad) * (3 - 4 * z_error_pr) / 3, 1.0],
-                    ],
-                ),
-            )
-            or np.allclose(
-                res.rho,
-                0.5
-                * np.array(
-                    [
-                        [1.0, np.exp(-1j * rad) * (3 - 4 * x_error_pr) / 3],
-                        [np.exp(1j * rad) * (3 - 4 * x_error_pr) / 3, 1.0],
-                    ],
-                ),
-            )
+        assert np.allclose(
+            res.rho,
+            0.5
+            * np.array(
+                [
+                    [1.0, np.exp(-1j * rad) * (3 - 4 * x_error_pr) * (3 - 4 * z_error_pr) / 9],
+                    [np.exp(1j * rad) * (3 - 4 * x_error_pr) * (3 - 4 * z_error_pr) / 9, 1.0],
+                ],
+            ),
         )
 
     # test measurement confuse outcome
-    def test_noisy_measure_confuse_rz(self, fx_rng: Generator) -> None:
+    @pytest.mark.parametrize("z_outcome,x_outcome", [(0, 0), (0, 1), (1, 0), (1, 1)])
+    def test_noisy_measure_confuse_rz(self, fx_rng: Generator, z_outcome: int, x_outcome: int) -> None:
         alpha = fx_rng.random()
         rzpattern = rzpat(alpha)
-        # probability 1 to shift both outcome
+
+        # M(0) determines Z, M(1) determines X
+        results = {}
+        cmd_count = 0
+        for cmd in rzpattern:
+            if cmd.kind.name == "M":
+                if cmd_count == 0:
+                    results[cmd.node] = z_outcome
+                elif cmd_count == 1:
+                    results[cmd.node] = x_outcome
+                cmd_count += 1
+
+        # Test with probability 1 to flip both outcomes
         res = rzpattern.simulate_pattern(
-            backend="densitymatrix", noise_model=DepolarisingNoiseModel(measure_error_prob=1.0), rng=fx_rng
+            backend="densitymatrix",
+            noise_model=DepolarisingNoiseModel(measure_error_prob=1.0),
+            branch_selector=FixedBranchSelector(results),
+            rng=fx_rng,
         )
-        # result X, XZ or Z
 
         exact = rz_exact_res(alpha)
-
         assert isinstance(res, DensityMatrix)
-        assert (
-            np.allclose(res.rho, Ops.X @ exact @ Ops.X)
-            or np.allclose(res.rho, Ops.Z @ exact @ Ops.Z)
-            or np.allclose(res.rho, Ops.Z @ Ops.X @ exact @ Ops.X @ Ops.Z)
-        )
+        # All outcomes lead to same result: both corrections applied due to flipping
+        assert np.allclose(res.rho, Ops.Z @ Ops.X @ exact @ Ops.X @ Ops.Z)
 
-        # arbitrary probability
+    @pytest.mark.parametrize("z_outcome,x_outcome", [(0, 0), (0, 1), (1, 0), (1, 1)])
+    def test_noisy_measure_confuse_rz_arbitrary(self, fx_rng: Generator, z_outcome: int, x_outcome: int) -> None:
+        alpha = fx_rng.random()
+        rzpattern = rzpat(alpha)
+
+        # M(0) determines Z, M(1) determines X
+        results = {}
+        cmd_count = 0
+        for cmd in rzpattern:
+            if cmd.kind.name == "M":
+                if cmd_count == 0:
+                    results[cmd.node] = z_outcome
+                elif cmd_count == 1:
+                    results[cmd.node] = x_outcome
+                cmd_count += 1
+
+        # Test with arbitrary probability
         measure_error_pr = fx_rng.random()
-        print(f"measure_error_pr = {measure_error_pr}")
+        print(f"measure_error_pr = {measure_error_pr}, z_outcome = {z_outcome}, x_outcome = {x_outcome}")
         res = rzpattern.simulate_pattern(
             backend="densitymatrix",
             noise_model=DepolarisingNoiseModel(measure_error_prob=measure_error_pr),
+            branch_selector=FixedBranchSelector(results),
             rng=fx_rng,
         )
-        # just add the case without readout errors
+
+        exact = rz_exact_res(alpha)
         assert isinstance(res, DensityMatrix)
+
+        # With arbitrary measure_error_pr, outcomes may or may not be flipped
+        # The physical result depends on whether the error occurs
+        # We check all possible cases
         assert (
             np.allclose(res.rho, exact)
             or np.allclose(res.rho, Ops.X @ exact @ Ops.X)

From 8803c11ab754645b4794778ceef6e11ec1062475 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tiberiu=20Sab=C4=83u?=
 <96194994+tibisabau@users.noreply.github.com>
Date: Fri, 6 Feb 2026 09:21:19 +0100
Subject: [PATCH 14/20] Update tests/test_noisy_density_matrix.py

Co-authored-by: thierry-martinez <thierry.martinez@inria.fr>
---
 tests/test_noisy_density_matrix.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/tests/test_noisy_density_matrix.py b/tests/test_noisy_density_matrix.py
index 87e41828d..aacb51ea9 100644
--- a/tests/test_noisy_density_matrix.py
+++ b/tests/test_noisy_density_matrix.py
@@ -319,7 +319,8 @@ def test_noisy_measure_channel_rz(self, fx_rng: Generator) -> None:
             ),
         )
 
-    @pytest.mark.parametrize("outcome_z,outcome_x", [(0, 0), (0, 1), (1, 0), (1, 1)])
+    @pytest.mark.parametrize("z_outcome", [0, 1])
+    @pytest.mark.parametrize("x_outcome", [0, 1])
     def test_noisy_x_rz(self, fx_rng: Generator, outcome_z: int, outcome_x: int) -> None:
         alpha = fx_rng.random()
         rzpattern = rzpat(alpha)

From efc76d51a574b580b677fec4011ef85a1908b1cb Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tiberiu=20Sab=C4=83u?=
 <96194994+tibisabau@users.noreply.github.com>
Date: Fri, 6 Feb 2026 09:21:42 +0100
Subject: [PATCH 15/20] Update tests/test_noisy_density_matrix.py

Co-authored-by: thierry-martinez <thierry.martinez@inria.fr>
---
 tests/test_noisy_density_matrix.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test_noisy_density_matrix.py b/tests/test_noisy_density_matrix.py
index aacb51ea9..be4671319 100644
--- a/tests/test_noisy_density_matrix.py
+++ b/tests/test_noisy_density_matrix.py
@@ -71,7 +71,7 @@ def test_noisy_measure_confuse_hadamard(self, fx_rng: Generator) -> None:
         assert np.allclose(res.rho, np.array([[0.0, 0.0], [0.0, 1.0]]))
 
     @pytest.mark.parametrize("outcome", [0, 1])
-    def test_noisy_measure_confuse_hadamard_arbitrary(self, fx_rng: Generator, outcome: int) -> None:
+    def test_noisy_measure_confuse_hadamard_arbitrary(self, fx_rng: Generator, outcome: Outcome) -> None:
         # arbitrary probability with fixed branch
         hadamardpattern = hpat()
         measure_error_pr = fx_rng.random()

From 931d3ec09bab2b3a6d188d24c3c148a64a08b49f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tiberiu=20Sab=C4=83u?=
 <96194994+tibisabau@users.noreply.github.com>
Date: Fri, 6 Feb 2026 09:22:07 +0100
Subject: [PATCH 16/20] Update tests/test_noisy_density_matrix.py

Co-authored-by: thierry-martinez <thierry.martinez@inria.fr>
---
 tests/test_noisy_density_matrix.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test_noisy_density_matrix.py b/tests/test_noisy_density_matrix.py
index be4671319..3ce9cc627 100644
--- a/tests/test_noisy_density_matrix.py
+++ b/tests/test_noisy_density_matrix.py
@@ -332,7 +332,7 @@ def test_noisy_x_rz(self, fx_rng: Generator, outcome_z: int, outcome_x: int) ->
         results = {}
         cmd_count = 0
         for cmd in rzpattern:
-            if cmd.kind.name == "M":
+            if cmd.kind == CommandKind.M:
                 if cmd_count == 0:
                     results[cmd.node] = outcome_z
                 elif cmd_count == 1:

From bfc66a0f56d17f7cfdd335bdc2b41583e77fc0ac Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tiberiu=20Sab=C4=83u?=
 <96194994+tibisabau@users.noreply.github.com>
Date: Fri, 6 Feb 2026 09:22:24 +0100
Subject: [PATCH 17/20] Update tests/test_noisy_density_matrix.py

Co-authored-by: thierry-martinez <thierry.martinez@inria.fr>
---
 tests/test_noisy_density_matrix.py | 11 ++---------
 1 file changed, 2 insertions(+), 9 deletions(-)

diff --git a/tests/test_noisy_density_matrix.py b/tests/test_noisy_density_matrix.py
index 3ce9cc627..9bdd2d438 100644
--- a/tests/test_noisy_density_matrix.py
+++ b/tests/test_noisy_density_matrix.py
@@ -329,15 +329,8 @@ def test_noisy_x_rz(self, fx_rng: Generator, outcome_z: int, outcome_x: int) ->
         print(f"x_error_pr = {x_error_pr}, outcome_z = {outcome_z}, outcome_x = {outcome_x}")
 
         # M(0) determines Z, M(1) determines X
-        results = {}
-        cmd_count = 0
-        for cmd in rzpattern:
-            if cmd.kind == CommandKind.M:
-                if cmd_count == 0:
-                    results[cmd.node] = outcome_z
-                elif cmd_count == 1:
-                    results[cmd.node] = outcome_x
-                cmd_count += 1
+        m_nodes = (cmd.node for cmd in rzpattern if cmd.kind == CommandKind.M)
+        results = {next(m_nodes): outcome_z, next(m_nodes): outcome_x}
 
         res = rzpattern.simulate_pattern(
             backend="densitymatrix",

From ce5e9a27a18f129df7ca872ce820e4bbf2e238a5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tiberiu=20Sab=C4=83u?= <T.Sabau@student.tudelft.nl>
Date: Fri, 6 Feb 2026 15:35:38 +0100
Subject: [PATCH 18/20] Fix ruff and mypy

---
 tests/test_noisy_density_matrix.py | 54 ++++++++++++++++--------------
 1 file changed, 29 insertions(+), 25 deletions(-)

diff --git a/tests/test_noisy_density_matrix.py b/tests/test_noisy_density_matrix.py
index 9bdd2d438..13f056af1 100644
--- a/tests/test_noisy_density_matrix.py
+++ b/tests/test_noisy_density_matrix.py
@@ -1,12 +1,13 @@
 from __future__ import annotations
 
-from typing import TYPE_CHECKING
+from typing import TYPE_CHECKING, cast
 
 import numpy as np
 import numpy.typing as npt
 import pytest
 
 from graphix.branch_selector import ConstBranchSelector, FixedBranchSelector
+from graphix.command import CommandKind, M
 from graphix.fundamentals import angle_to_rad
 from graphix.noise_models import DepolarisingNoiseModel
 from graphix.noise_models.noise_model import NoiselessNoiseModel
@@ -18,6 +19,7 @@
     from numpy.random import Generator
 
     from graphix.fundamentals import Angle
+    from graphix.measurements import Outcome
     from graphix.pattern import Pattern
 
 
@@ -110,7 +112,7 @@ def test_noisy_measure_channel_hadamard(self, fx_rng: Generator) -> None:
 
     # test Pauli X error
     @pytest.mark.parametrize("outcome", [0, 1])
-    def test_noisy_x_hadamard(self, fx_rng: Generator, outcome: int) -> None:
+    def test_noisy_x_hadamard(self, fx_rng: Generator, outcome: Outcome) -> None:
         hadamardpattern = hpat()
         # x error only
         x_error_pr = fx_rng.random()
@@ -321,16 +323,16 @@ def test_noisy_measure_channel_rz(self, fx_rng: Generator) -> None:
 
     @pytest.mark.parametrize("z_outcome", [0, 1])
     @pytest.mark.parametrize("x_outcome", [0, 1])
-    def test_noisy_x_rz(self, fx_rng: Generator, outcome_z: int, outcome_x: int) -> None:
+    def test_noisy_x_rz(self, fx_rng: Generator, z_outcome: Outcome, x_outcome: Outcome) -> None:
         alpha = fx_rng.random()
         rzpattern = rzpat(alpha)
         # x error only
         x_error_pr = fx_rng.random()
-        print(f"x_error_pr = {x_error_pr}, outcome_z = {outcome_z}, outcome_x = {outcome_x}")
+        print(f"x_error_pr = {x_error_pr}, outcome_z = {z_outcome}, outcome_x = {x_outcome}")
 
         # M(0) determines Z, M(1) determines X
         m_nodes = (cmd.node for cmd in rzpattern if cmd.kind == CommandKind.M)
-        results = {next(m_nodes): outcome_z, next(m_nodes): outcome_x}
+        results: dict[int, Outcome] = {next(m_nodes): z_outcome, next(m_nodes): x_outcome}
 
         res = rzpattern.simulate_pattern(
             backend="densitymatrix",
@@ -353,8 +355,8 @@ def test_noisy_x_rz(self, fx_rng: Generator, outcome_z: int, outcome_x: int) ->
             ),
         )
 
-    @pytest.mark.parametrize("outcome_z,outcome_x", [(0, 0), (0, 1), (1, 0), (1, 1)])
-    def test_noisy_z_rz(self, fx_rng: Generator, outcome_z: int, outcome_x: int) -> None:
+    @pytest.mark.parametrize(("outcome_z", "outcome_x"), [(0, 0), (0, 1), (1, 0), (1, 1)])
+    def test_noisy_z_rz(self, fx_rng: Generator, outcome_z: Outcome, outcome_x: Outcome) -> None:
         alpha = fx_rng.random()
         rzpattern = rzpat(alpha)
         # z error only
@@ -362,14 +364,14 @@ def test_noisy_z_rz(self, fx_rng: Generator, outcome_z: int, outcome_x: int) ->
         print(f"z_error_pr = {z_error_pr}, outcome_z = {outcome_z}, outcome_x = {outcome_x}")
 
         # M(0) determines Z, M(1) determines X
-        results = {}
+        results: dict[int, Outcome] = {}
         cmd_count = 0
         for cmd in rzpattern:
             if cmd.kind.name == "M":
                 if cmd_count == 0:
-                    results[cmd.node] = outcome_z
+                    results[cast("M", cmd).node] = outcome_z
                 elif cmd_count == 1:
-                    results[cmd.node] = outcome_x
+                    results[cast("M", cmd).node] = outcome_x
                 cmd_count += 1
 
         res = rzpattern.simulate_pattern(
@@ -393,8 +395,8 @@ def test_noisy_z_rz(self, fx_rng: Generator, outcome_z: int, outcome_x: int) ->
             ),
         )
 
-    @pytest.mark.parametrize("z_outcome,x_outcome", [(0, 0), (0, 1), (1, 0), (1, 1)])
-    def test_noisy_xz_rz(self, fx_rng: Generator, z_outcome: int, x_outcome: int) -> None:
+    @pytest.mark.parametrize(("z_outcome", "x_outcome"), [(0, 0), (0, 1), (1, 0), (1, 1)])
+    def test_noisy_xz_rz(self, fx_rng: Generator, z_outcome: Outcome, x_outcome: Outcome) -> None:
         alpha = fx_rng.random()
         rzpattern = rzpat(alpha)
         # x and z errors
@@ -405,14 +407,14 @@ def test_noisy_xz_rz(self, fx_rng: Generator, z_outcome: int, x_outcome: int) ->
         print(f"z_outcome = {z_outcome}, x_outcome = {x_outcome}")
 
         # M(0) determines Z correction, M(1) determines X correction
-        results = {}
+        results: dict[int, Outcome] = {}
         cmd_count = 0
         for cmd in rzpattern:
             if cmd.kind.name == "M":
                 if cmd_count == 0:
-                    results[cmd.node] = z_outcome
+                    results[cast("M", cmd).node] = z_outcome
                 elif cmd_count == 1:
-                    results[cmd.node] = x_outcome
+                    results[cast("M", cmd).node] = x_outcome
                 cmd_count += 1
 
         res = rzpattern.simulate_pattern(
@@ -437,20 +439,20 @@ def test_noisy_xz_rz(self, fx_rng: Generator, z_outcome: int, x_outcome: int) ->
         )
 
     # test measurement confuse outcome
-    @pytest.mark.parametrize("z_outcome,x_outcome", [(0, 0), (0, 1), (1, 0), (1, 1)])
-    def test_noisy_measure_confuse_rz(self, fx_rng: Generator, z_outcome: int, x_outcome: int) -> None:
+    @pytest.mark.parametrize(("z_outcome", "x_outcome"), [(0, 0), (0, 1), (1, 0), (1, 1)])
+    def test_noisy_measure_confuse_rz(self, fx_rng: Generator, z_outcome: Outcome, x_outcome: Outcome) -> None:
         alpha = fx_rng.random()
         rzpattern = rzpat(alpha)
 
         # M(0) determines Z, M(1) determines X
-        results = {}
+        results: dict[int, Outcome] = {}
         cmd_count = 0
         for cmd in rzpattern:
             if cmd.kind.name == "M":
                 if cmd_count == 0:
-                    results[cmd.node] = z_outcome
+                    results[cast("M", cmd).node] = z_outcome
                 elif cmd_count == 1:
-                    results[cmd.node] = x_outcome
+                    results[cast("M", cmd).node] = x_outcome
                 cmd_count += 1
 
         # Test with probability 1 to flip both outcomes
@@ -466,20 +468,22 @@ def test_noisy_measure_confuse_rz(self, fx_rng: Generator, z_outcome: int, x_out
         # All outcomes lead to same result: both corrections applied due to flipping
         assert np.allclose(res.rho, Ops.Z @ Ops.X @ exact @ Ops.X @ Ops.Z)
 
-    @pytest.mark.parametrize("z_outcome,x_outcome", [(0, 0), (0, 1), (1, 0), (1, 1)])
-    def test_noisy_measure_confuse_rz_arbitrary(self, fx_rng: Generator, z_outcome: int, x_outcome: int) -> None:
+    @pytest.mark.parametrize(("z_outcome", "x_outcome"), [(0, 0), (0, 1), (1, 0), (1, 1)])
+    def test_noisy_measure_confuse_rz_arbitrary(
+        self, fx_rng: Generator, z_outcome: Outcome, x_outcome: Outcome
+    ) -> None:
         alpha = fx_rng.random()
         rzpattern = rzpat(alpha)
 
         # M(0) determines Z, M(1) determines X
-        results = {}
+        results: dict[int, Outcome] = {}
         cmd_count = 0
         for cmd in rzpattern:
             if cmd.kind.name == "M":
                 if cmd_count == 0:
-                    results[cmd.node] = z_outcome
+                    results[cast("M", cmd).node] = z_outcome
                 elif cmd_count == 1:
-                    results[cmd.node] = x_outcome
+                    results[cast("M", cmd).node] = x_outcome
                 cmd_count += 1
 
         # Test with arbitrary probability

From 4acdf57e585a9038950b617a29aa5b5f684296b6 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tiberiu=20Sab=C4=83u?= <T.Sabau@student.tudelft.nl>
Date: Fri, 6 Feb 2026 17:24:23 +0100
Subject: [PATCH 19/20] fix: Remove loops and cast

---
 tests/test_noisy_density_matrix.py | 56 ++++++++----------------------
 1 file changed, 14 insertions(+), 42 deletions(-)

diff --git a/tests/test_noisy_density_matrix.py b/tests/test_noisy_density_matrix.py
index 13f056af1..e4d867b07 100644
--- a/tests/test_noisy_density_matrix.py
+++ b/tests/test_noisy_density_matrix.py
@@ -1,13 +1,13 @@
 from __future__ import annotations
 
-from typing import TYPE_CHECKING, cast
+from typing import TYPE_CHECKING
 
 import numpy as np
 import numpy.typing as npt
 import pytest
 
 from graphix.branch_selector import ConstBranchSelector, FixedBranchSelector
-from graphix.command import CommandKind, M
+from graphix.command import CommandKind
 from graphix.fundamentals import angle_to_rad
 from graphix.noise_models import DepolarisingNoiseModel
 from graphix.noise_models.noise_model import NoiselessNoiseModel
@@ -355,7 +355,8 @@ def test_noisy_x_rz(self, fx_rng: Generator, z_outcome: Outcome, x_outcome: Outc
             ),
         )
 
-    @pytest.mark.parametrize(("outcome_z", "outcome_x"), [(0, 0), (0, 1), (1, 0), (1, 1)])
+    @pytest.mark.parametrize("outcome_z", [0, 1])
+    @pytest.mark.parametrize("outcome_x", [0, 1])
     def test_noisy_z_rz(self, fx_rng: Generator, outcome_z: Outcome, outcome_x: Outcome) -> None:
         alpha = fx_rng.random()
         rzpattern = rzpat(alpha)
@@ -364,15 +365,7 @@ def test_noisy_z_rz(self, fx_rng: Generator, outcome_z: Outcome, outcome_x: Outc
         print(f"z_error_pr = {z_error_pr}, outcome_z = {outcome_z}, outcome_x = {outcome_x}")
 
         # M(0) determines Z, M(1) determines X
-        results: dict[int, Outcome] = {}
-        cmd_count = 0
-        for cmd in rzpattern:
-            if cmd.kind.name == "M":
-                if cmd_count == 0:
-                    results[cast("M", cmd).node] = outcome_z
-                elif cmd_count == 1:
-                    results[cast("M", cmd).node] = outcome_x
-                cmd_count += 1
+        results: dict[int, Outcome] = {0: outcome_z, 1: outcome_x}
 
         res = rzpattern.simulate_pattern(
             backend="densitymatrix",
@@ -395,7 +388,8 @@ def test_noisy_z_rz(self, fx_rng: Generator, outcome_z: Outcome, outcome_x: Outc
             ),
         )
 
-    @pytest.mark.parametrize(("z_outcome", "x_outcome"), [(0, 0), (0, 1), (1, 0), (1, 1)])
+    @pytest.mark.parametrize("z_outcome", [0, 1])
+    @pytest.mark.parametrize("x_outcome", [0, 1])
     def test_noisy_xz_rz(self, fx_rng: Generator, z_outcome: Outcome, x_outcome: Outcome) -> None:
         alpha = fx_rng.random()
         rzpattern = rzpat(alpha)
@@ -407,15 +401,7 @@ def test_noisy_xz_rz(self, fx_rng: Generator, z_outcome: Outcome, x_outcome: Out
         print(f"z_outcome = {z_outcome}, x_outcome = {x_outcome}")
 
         # M(0) determines Z correction, M(1) determines X correction
-        results: dict[int, Outcome] = {}
-        cmd_count = 0
-        for cmd in rzpattern:
-            if cmd.kind.name == "M":
-                if cmd_count == 0:
-                    results[cast("M", cmd).node] = z_outcome
-                elif cmd_count == 1:
-                    results[cast("M", cmd).node] = x_outcome
-                cmd_count += 1
+        results: dict[int, Outcome] = {0: z_outcome, 1: x_outcome}
 
         res = rzpattern.simulate_pattern(
             backend="densitymatrix",
@@ -439,21 +425,14 @@ def test_noisy_xz_rz(self, fx_rng: Generator, z_outcome: Outcome, x_outcome: Out
         )
 
     # test measurement confuse outcome
-    @pytest.mark.parametrize(("z_outcome", "x_outcome"), [(0, 0), (0, 1), (1, 0), (1, 1)])
+    @pytest.mark.parametrize("z_outcome", [0, 1])
+    @pytest.mark.parametrize("x_outcome", [0, 1])
     def test_noisy_measure_confuse_rz(self, fx_rng: Generator, z_outcome: Outcome, x_outcome: Outcome) -> None:
         alpha = fx_rng.random()
         rzpattern = rzpat(alpha)
 
         # M(0) determines Z, M(1) determines X
-        results: dict[int, Outcome] = {}
-        cmd_count = 0
-        for cmd in rzpattern:
-            if cmd.kind.name == "M":
-                if cmd_count == 0:
-                    results[cast("M", cmd).node] = z_outcome
-                elif cmd_count == 1:
-                    results[cast("M", cmd).node] = x_outcome
-                cmd_count += 1
+        results: dict[int, Outcome] = {0: z_outcome, 1: x_outcome}
 
         # Test with probability 1 to flip both outcomes
         res = rzpattern.simulate_pattern(
@@ -468,7 +447,8 @@ def test_noisy_measure_confuse_rz(self, fx_rng: Generator, z_outcome: Outcome, x
         # All outcomes lead to same result: both corrections applied due to flipping
         assert np.allclose(res.rho, Ops.Z @ Ops.X @ exact @ Ops.X @ Ops.Z)
 
-    @pytest.mark.parametrize(("z_outcome", "x_outcome"), [(0, 0), (0, 1), (1, 0), (1, 1)])
+    @pytest.mark.parametrize("z_outcome", [0, 1])
+    @pytest.mark.parametrize("x_outcome", [0, 1])
     def test_noisy_measure_confuse_rz_arbitrary(
         self, fx_rng: Generator, z_outcome: Outcome, x_outcome: Outcome
     ) -> None:
@@ -476,15 +456,7 @@ def test_noisy_measure_confuse_rz_arbitrary(
         rzpattern = rzpat(alpha)
 
         # M(0) determines Z, M(1) determines X
-        results: dict[int, Outcome] = {}
-        cmd_count = 0
-        for cmd in rzpattern:
-            if cmd.kind.name == "M":
-                if cmd_count == 0:
-                    results[cast("M", cmd).node] = z_outcome
-                elif cmd_count == 1:
-                    results[cast("M", cmd).node] = x_outcome
-                cmd_count += 1
+        results: dict[int, Outcome] = {0: z_outcome, 1: x_outcome}
 
         # Test with arbitrary probability
         measure_error_pr = fx_rng.random()

From b8c62a79ddc9176cb400783871ebdf9628d82292 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tiberiu=20Sab=C4=83u?= <T.Sabau@student.tudelft.nl>
Date: Sat, 14 Feb 2026 23:58:22 +0100
Subject: [PATCH 20/20] Update test_noisy_density_matrix.py based on issue 428
 fix

---
 tests/test_noisy_density_matrix.py | 120 ++++++++++++++++++++---------
 1 file changed, 82 insertions(+), 38 deletions(-)

diff --git a/tests/test_noisy_density_matrix.py b/tests/test_noisy_density_matrix.py
index e4d867b07..dc207931b 100644
--- a/tests/test_noisy_density_matrix.py
+++ b/tests/test_noisy_density_matrix.py
@@ -123,12 +123,17 @@ def test_noisy_x_hadamard(self, fx_rng: Generator, outcome: Outcome) -> None:
             branch_selector=ConstBranchSelector(outcome),
             rng=fx_rng,
         )
-        # Both outcomes lead to |0> after correction. ApplyNoise is unconditional, so noise is always applied.
+        # Pattern has X(1, {0}), so X error noise only applied when outcome=1
         assert isinstance(res, DensityMatrix)
-        assert np.allclose(
-            res.rho,
-            np.array([[1 - 2 * x_error_pr / 3.0, 0.0], [0.0, 2 * x_error_pr / 3.0]]),
-        )
+        if outcome == 0:
+            # No X correction → no X error noise
+            assert np.allclose(res.rho, np.array([[1.0, 0.0], [0.0, 0.0]]))
+        else:
+            # X correction applied → X error noise applied
+            assert np.allclose(
+                res.rho,
+                np.array([[1 - 2 * x_error_pr / 3.0, 0.0], [0.0, 2 * x_error_pr / 3.0]]),
+            )
 
     # test entanglement error
     def test_noisy_entanglement_hadamard(self, fx_rng: Generator) -> None:
@@ -341,19 +346,24 @@ def test_noisy_x_rz(self, fx_rng: Generator, z_outcome: Outcome, x_outcome: Outc
             rng=fx_rng,
         )
 
-        # All outcomes lead to the same state after corrections. ApplyNoise is unconditional, so noise is always applied.
+        # Pattern has X(2, {1}), so X error noise only applied when x_outcome=1
         assert isinstance(res, DensityMatrix)
         rad = angle_to_rad(alpha)
-        assert np.allclose(
-            res.rho,
-            0.5
-            * np.array(
-                [
-                    [1.0, np.exp(-1j * rad) * (3 - 4 * x_error_pr) / 3],
-                    [np.exp(1j * rad) * (3 - 4 * x_error_pr) / 3, 1.0],
-                ],
-            ),
-        )
+        if x_outcome == 0:
+            # No X correction → no X error noise
+            assert np.allclose(res.rho, rz_exact_res(alpha))
+        else:
+            # X correction applied → X error noise applied
+            assert np.allclose(
+                res.rho,
+                0.5
+                * np.array(
+                    [
+                        [1.0, np.exp(-1j * rad) * (3 - 4 * x_error_pr) / 3],
+                        [np.exp(1j * rad) * (3 - 4 * x_error_pr) / 3, 1.0],
+                    ],
+                ),
+            )
 
     @pytest.mark.parametrize("outcome_z", [0, 1])
     @pytest.mark.parametrize("outcome_x", [0, 1])
@@ -374,19 +384,24 @@ def test_noisy_z_rz(self, fx_rng: Generator, outcome_z: Outcome, outcome_x: Outc
             rng=fx_rng,
         )
 
-        # All outcomes lead to the same state after corrections. ApplyNoise is unconditional, so noise is always applied.
+        # Pattern has Z(2, {0}), so Z error noise only applied when outcome_z=1
         assert isinstance(res, DensityMatrix)
         rad = angle_to_rad(alpha)
-        assert np.allclose(
-            res.rho,
-            0.5
-            * np.array(
-                [
-                    [1.0, np.exp(-1j * rad) * (3 - 4 * z_error_pr) / 3],
-                    [np.exp(1j * rad) * (3 - 4 * z_error_pr) / 3, 1.0],
-                ],
-            ),
-        )
+        if outcome_z == 0:
+            # No Z correction → no Z error noise
+            assert np.allclose(res.rho, rz_exact_res(alpha))
+        else:
+            # Z correction applied → Z error noise applied
+            assert np.allclose(
+                res.rho,
+                0.5
+                * np.array(
+                    [
+                        [1.0, np.exp(-1j * rad) * (3 - 4 * z_error_pr) / 3],
+                        [np.exp(1j * rad) * (3 - 4 * z_error_pr) / 3, 1.0],
+                    ],
+                ),
+            )
 
     @pytest.mark.parametrize("z_outcome", [0, 1])
     @pytest.mark.parametrize("x_outcome", [0, 1])
@@ -410,19 +425,48 @@ def test_noisy_xz_rz(self, fx_rng: Generator, z_outcome: Outcome, x_outcome: Out
             rng=fx_rng,
         )
 
-        # All outcomes lead to same state after corrections. Both X and Z noise applied unconditionally.
+        # Pattern has X(2, {1}) and Z(2, {0}), noise applied conditionally
         assert isinstance(res, DensityMatrix)
         rad = angle_to_rad(alpha)
-        assert np.allclose(
-            res.rho,
-            0.5
-            * np.array(
-                [
-                    [1.0, np.exp(-1j * rad) * (3 - 4 * x_error_pr) * (3 - 4 * z_error_pr) / 9],
-                    [np.exp(1j * rad) * (3 - 4 * x_error_pr) * (3 - 4 * z_error_pr) / 9, 1.0],
-                ],
-            ),
-        )
+        if z_outcome == 0 and x_outcome == 0:
+            # No corrections → no noise
+            assert np.allclose(res.rho, rz_exact_res(alpha))
+        elif z_outcome == 0 and x_outcome == 1:
+            # Only X correction → only X noise
+            assert np.allclose(
+                res.rho,
+                0.5
+                * np.array(
+                    [
+                        [1.0, np.exp(-1j * rad) * (3 - 4 * x_error_pr) / 3],
+                        [np.exp(1j * rad) * (3 - 4 * x_error_pr) / 3, 1.0],
+                    ],
+                ),
+            )
+        elif z_outcome == 1 and x_outcome == 0:
+            # Only Z correction → only Z noise
+            assert np.allclose(
+                res.rho,
+                0.5
+                * np.array(
+                    [
+                        [1.0, np.exp(-1j * rad) * (3 - 4 * z_error_pr) / 3],
+                        [np.exp(1j * rad) * (3 - 4 * z_error_pr) / 3, 1.0],
+                    ],
+                ),
+            )
+        else:  # z_outcome == 1 and x_outcome == 1
+            # Both corrections → both noises
+            assert np.allclose(
+                res.rho,
+                0.5
+                * np.array(
+                    [
+                        [1.0, np.exp(-1j * rad) * (3 - 4 * x_error_pr) * (3 - 4 * z_error_pr) / 9],
+                        [np.exp(1j * rad) * (3 - 4 * x_error_pr) * (3 - 4 * z_error_pr) / 9, 1.0],
+                    ],
+                ),
+            )
 
     # test measurement confuse outcome
     @pytest.mark.parametrize("z_outcome", [0, 1])