From e2fc2388af5482746835f0d01666cef5f7917f4a Mon Sep 17 00:00:00 2001
From: Thierry Martinez <thierry.martinez@inria.fr>
Date: Tue, 10 Jun 2025 15:41:58 +0200
Subject: [PATCH 01/19] Branch selection

This commit adds an abstract class `BranchSelector` and the following
branch selection strategies:

- `RandomBranchSelector`, which corresponds to the strategy that was
  already implemented.

- `FixedBranchSelector`, for a selector where a dictionary sets the
  outcome for each measurement.

- `ConstBranchSelector`, for a selector that always returns the same
  outcome, whatever the measurement is.
---
 CHANGELOG.md                 |  10 +++
 docs/source/simulator.rst    |  23 ++++++
 graphix/branch_selector.py   | 137 +++++++++++++++++++++++++++++++++++
 graphix/sim/base_backend.py  |  65 +++++++++++------
 graphix/transpiler.py        |  25 ++++++-
 tests/test_density_matrix.py |   3 +-
 tests/test_pattern.py        |  28 +++----
 7 files changed, 247 insertions(+), 44 deletions(-)
 create mode 100644 graphix/branch_selector.py

diff --git a/CHANGELOG.md b/CHANGELOG.md
index c133b924b..b2dc7e1ea 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -12,6 +12,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - #277: Methods for pretty-printing `Pattern`: `to_ascii`,
   `to_unicode`, `to_latex`.
 
+- #299: Branch selection in simulation: in addition to
+  `RandomBranchSelector` which corresponds to the strategy that was
+  already implemented, the user can use `FixedBranchSelector`,
+  `ConstBranchSelector`, or define a custom branch selection by
+  deriving the abstract class `BranchSelector`.
+
 ### Fixed
 
 - #277: The result of `repr()` for `Pattern`, `Circuit`, `Command`,
@@ -27,6 +33,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - #277: The method `Pattern.print_pattern` is now deprecated.
 
+- #299: `pr_calc` parameter is deprecated in back-end initializers.
+  The user can specify `pr_calc` in the constructor of
+  `RandomBranchSelector` instead.
+
 ## [0.3.1] - 2025-04-21
 
 ### Added
diff --git a/docs/source/simulator.rst b/docs/source/simulator.rst
index da35237a1..fd4d6c97f 100644
--- a/docs/source/simulator.rst
+++ b/docs/source/simulator.rst
@@ -12,6 +12,29 @@ Pattern Simulation
 
     .. automethod:: run
 
+Branch Selection (:mod:`graphix.branch_selector` module)
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+.. currentmodule:: graphix.branch_selector
+
+.. autoclass:: BranchSelector
+
+    .. automethod:: measure
+
+.. autoclass:: RandomBranchSelector
+    :members:
+
+    .. automethod:: measure
+
+.. autoclass:: FixedBranchSelector
+    :members:
+
+    .. automethod:: measure
+
+.. autoclass:: ConstBranchSelector
+    :members:
+
+    .. automethod:: measure
 
 Simulator backends
 ++++++++++++++++++
diff --git a/graphix/branch_selector.py b/graphix/branch_selector.py
new file mode 100644
index 000000000..49b130ece
--- /dev/null
+++ b/graphix/branch_selector.py
@@ -0,0 +1,137 @@
+"""Branch selector.
+
+Branch selectors determine the computation branch that is explored
+during a simulation, meaning the choice of measurement outcomes.  The
+branch selection can be random (see :class:`RandomBranchSelector`) or
+deterministic (see :class:`ConstBranchSelector`).
+
+"""
+
+from __future__ import annotations
+
+from abc import ABC, abstractmethod
+from dataclasses import dataclass
+from typing import TYPE_CHECKING, Callable
+
+from typing_extensions import override
+
+from graphix.rng import ensure_rng
+
+if TYPE_CHECKING:
+    from collections.abc import Mapping
+
+    from numpy.random import Generator
+
+
+class BranchSelector(ABC):
+    """Abstract class for branch selectors.
+
+    A branch selector provides the method `measure`, which returns the
+    measurement outcome (0 or 1) for a given qubit.
+    """
+
+    @abstractmethod
+    def measure(self, qubit: int, compute_expectation_0: Callable[[], float]) -> bool:
+        """
+        Return the measurement outcome of `qubit`.
+
+        This method may use the function `compute_expectation_0` to
+        retrieve the expected probability of outcome 0. The probability is
+        computed only if this function is called (lazy computation),
+        ensuring no unnecessary computational cost.
+        """
+
+
+@dataclass
+class RandomBranchSelector(BranchSelector):
+    """Random branch selector.
+
+    Parameters
+    ----------
+    pr_calc : bool, optional
+        Whether to compute the probability distribution before selecting the measurement result.
+        If False, measurements yield 0/1 with equal probability (50% each).
+        Default is `True`.
+    rng : Generator | None, optional
+        Random-number generator for measurements.
+        If `None`, a default random-number generator is used.
+        Default is `None`.
+    """
+
+    pr_calc: bool = True
+    rng: Generator | None = None
+
+    @override
+    def measure(self, qubit: int, compute_expectation_0: Callable[[], float]) -> bool:
+        """
+        Return the measurement outcome of `qubit`.
+
+        If `pr_calc` is `True`, the measurement outcome is determined based on the
+        computed probability of outcome 0. Otherwise, the result is randomly chosen
+        with a 50% chance for either outcome.
+        """
+        self.rng = ensure_rng(self.rng)
+        if self.pr_calc:
+            prob_0 = compute_expectation_0()
+            return self.rng.random() > prob_0
+        outcome: int = self.rng.choice([0, 1])
+        return outcome == 1
+
+
+@dataclass
+class FixedBranchSelector(BranchSelector):
+    """Branch selector with predefined measurement outcomes.
+
+    The mapping is fixed in `results`. By default, an error is raised if
+    a qubit is measured without a predefined outcome. However, another
+    branch selector can be specified in `default` to handle such cases.
+
+    Parameters
+    ----------
+    results : Mapping[int, bool]
+        A dictionary mapping qubits to their measurement outcomes.
+        If a qubit is not present in this mapping, the `default` branch
+        selector is used.
+    default : BranchSelector | None, optional
+        Branch selector to use for qubits not present in `results`.
+        If `None`, an error is raised when an unmapped qubit is measured.
+        Default is `None`.
+    """
+
+    results: Mapping[int, bool]
+    default: BranchSelector | None = None
+
+    @override
+    def measure(self, qubit: int, compute_expectation_0: Callable[[], float]) -> bool:
+        """
+        Return the predefined measurement outcome of `qubit`, if available.
+
+        If the qubit is not present in `results`, the `default` branch selector
+        is used. If no default is provided, an error is raised.
+        """
+        result = self.results.get(qubit)
+        if result is None:
+            if self.default is None:
+                raise ValueError(f"Unexpected measurement of qubit {qubit}.")
+            return self.default.measure(qubit, compute_expectation_0)
+        return result
+
+
+@dataclass
+class ConstBranchSelector(BranchSelector):
+    """Branch selector with a constant measurement outcome.
+
+    The value `result` is returned for every qubit.
+
+    Parameters
+    ----------
+    result : bool
+        The fixed measurement outcome for all qubits.
+    """
+
+    result: bool
+
+    @override
+    def measure(self, qubit: int, compute_expectation_0: Callable[[], float]) -> bool:
+        """Return the constant measurement outcome `result` for any qubit."""
+        return self.result
diff --git a/graphix/sim/base_backend.py b/graphix/sim/base_backend.py
index 4dfcd0e7b..e0d6a483d 100644
--- a/graphix/sim/base_backend.py
+++ b/graphix/sim/base_backend.py
@@ -2,11 +2,13 @@
 
 from __future__ import annotations
 
+import warnings
 from abc import ABC, abstractmethod
 from typing import TYPE_CHECKING
 
 import numpy as np
 
+from graphix.branch_selector import BranchSelector, RandomBranchSelector
 from graphix.clifford import Clifford
 from graphix.command import CommandKind
 from graphix.ops import Ops
@@ -160,21 +162,32 @@ def _op_mat_from_result(vec: tuple[float, float, float], result: bool, symbolic:
 
 
 def perform_measure(
-    qubit: int, plane: Plane, angle: float, state, rng, pr_calc: bool = True, symbolic: bool = False
+    qubit_node: int,
+    qubit_loc: int,
+    plane: Plane,
+    angle: float,
+    state,
+    branch_selector: BranchSelector,
+    symbolic: bool = False,
 ) -> npt.NDArray:
     """Perform measurement of a qubit."""
     vec = plane.polar(angle)
-    if pr_calc:
-        op_mat = _op_mat_from_result(vec, False, symbolic=symbolic)
-        prob_0 = state.expectation_single(op_mat, qubit)
-        result = rng.random() > prob_0
-        if result:
-            op_mat = _op_mat_from_result(vec, True, symbolic=symbolic)
-    else:
-        # choose the measurement result randomly
-        result = rng.choice([0, 1])
-        op_mat = _op_mat_from_result(vec, result, symbolic=symbolic)
-    state.evolve_single(op_mat, qubit)
+    # op_mat_0 may contain the matrix operator associated with the outcome 0,
+    # but the value is computed lazily, i.e., only if needed.
+    op_mat_0 = None
+
+    def get_op_mat_0() -> np.ndarray:
+        nonlocal op_mat_0
+        if op_mat_0 is None:
+            op_mat_0 = _op_mat_from_result(vec, False, symbolic=symbolic)
+        return op_mat_0
+
+    def compute_expectation_0() -> float:
+        return state.expectation_single(get_op_mat_0(), qubit_loc)
+
+    result = branch_selector.measure(qubit_node, compute_expectation_0)
+    op_mat = _op_mat_from_result(vec, True, symbolic=symbolic) if result else get_op_mat_0()
+    state.evolve_single(op_mat, qubit_loc)
     return result
 
 
@@ -185,7 +198,8 @@ def __init__(
         self,
         state: State,
         node_index: NodeIndex | None = None,
-        pr_calc: bool = True,
+        branch_selector: BranchSelector | None = None,
+        pr_calc: bool | None = None,
         rng: Generator | None = None,
         symbolic: bool = False,
     ):
@@ -211,17 +225,24 @@ def __init__(
             self.__node_index = NodeIndex()
         else:
             self.__node_index = node_index.copy()
-        if not isinstance(pr_calc, bool):
-            raise TypeError("`pr_calc` should be bool")
         # whether to compute the probability
-        self.__pr_calc = pr_calc
+        if branch_selector is None:
+            if pr_calc is None:
+                pr_calc = True
+            else:
+                warnings.warn(
+                    "Setting `pr_calc` in `Backend` is deprecated. Use a `RandomBranchSelector` instead.",
+                    DeprecationWarning,
+                    stacklevel=1,
+                )
+            self.__branch_selector: BranchSelector = RandomBranchSelector(pr_calc=pr_calc, rng=rng)
+        else:
+            if pr_calc is not None or rng is not None:
+                raise ValueError("Cannot specify both branch selector and pr_calc/rng")
+            self.__branch_selector = branch_selector
         self.__rng = ensure_rng(rng)
         self.__symbolic = symbolic
 
-    def copy(self) -> Backend:
-        """Return a copy of the backend."""
-        return Backend(self.__state, self.__node_index, self.__pr_calc, self.__rng)
-
     @property
     def rng(self) -> Generator:
         """Return the associated random-number generator."""
@@ -274,12 +295,12 @@ def measure(self, node: int, measurement: Measurement) -> bool:
         """
         loc = self.node_index.index(node)
         result = perform_measure(
+            node,
             loc,
             measurement.plane,
             measurement.angle,
             self.state,
-            rng=self.__rng,
-            pr_calc=self.__pr_calc,
+            self.__branch_selector,
             symbolic=self.__symbolic,
         )
         self.node_index.remove(node)
diff --git a/graphix/transpiler.py b/graphix/transpiler.py
index 1de8b05c3..63624e97e 100644
--- a/graphix/transpiler.py
+++ b/graphix/transpiler.py
@@ -14,6 +14,7 @@
 from typing_extensions import assert_never
 
 from graphix import command, instruction, parameter
+from graphix.branch_selector import BranchSelector, RandomBranchSelector
 from graphix.command import CommandKind, E, M, N, X, Z
 from graphix.fundamentals import Plane
 from graphix.instruction import Instruction, InstructionKind
@@ -26,6 +27,8 @@
 if TYPE_CHECKING:
     from collections.abc import Iterable, Mapping, Sequence
 
+    from numpy.random import Generator
+
 
 @dataclasses.dataclass
 class TranspileResult:
@@ -885,18 +888,34 @@ def _sort_outputs(cls, pattern: Pattern, output_nodes: Sequence[int]):
             elif cmd.nodes in old_out:
                 cmd.nodes = output_nodes[old_out.index(cmd.nodes)]
 
-    def simulate_statevector(self, input_state: Data | None = None) -> SimulateResult:
+    def simulate_statevector(
+        self,
+        input_state: Data | None = None,
+        branch_selector: BranchSelector | None = None,
+        rng: Generator | None = None,
+    ) -> SimulateResult:
         """Run statevector simulation of the gate sequence.
 
         Parameters
         ----------
         input_state : :class:`graphix.sim.statevec.Statevec`
 
+        branch_selector: :class:`graphix.sim.base_backend.BranchSelector`
+            branch selector for measures (default: `graphix.sim.base_backend.RandomBranchSelector()`)
+
+        rng: :class:`np.random.Generator`
+            random number generator for `RandomBranchSelector` (should only be used with default branch selector)
+
         Returns
         -------
         result : :class:`SimulateResult`
             output state of the statevector simulation and results of classical measures.
         """
+        if branch_selector is None:
+            branch_selector = RandomBranchSelector(rng=rng)
+        elif rng is not None:
+            raise ValueError("Cannot specify both branch selector and rng")
+
         state = Statevec(nqubit=self.width) if input_state is None else Statevec(nqubit=self.width, data=input_state)
 
         classical_measures = []
@@ -931,7 +950,9 @@ def simulate_statevector(self, input_state: Data | None = None) -> SimulateResul
             elif kind == instruction.InstructionKind.CCX:
                 state.evolve(Ops.CCX, [instr.controls[0], instr.controls[1], instr.target])
             elif kind == instruction.InstructionKind.M:
-                result = base_backend.perform_measure(instr.target, instr.plane, instr.angle * np.pi, state, np.random)
+                result = base_backend.perform_measure(
+                    instr.target, instr.target, instr.plane, instr.angle * np.pi, state, branch_selector
+                )
                 classical_measures.append(result)
             else:
                 raise ValueError(f"Unknown instruction: {instr}")
diff --git a/tests/test_density_matrix.py b/tests/test_density_matrix.py
index 38bc348d5..37354b67a 100644
--- a/tests/test_density_matrix.py
+++ b/tests/test_density_matrix.py
@@ -10,6 +10,7 @@
 import pytest
 
 import graphix.random_objects as randobj
+from graphix.branch_selector import RandomBranchSelector
 from graphix.channels import KrausChannel, dephasing_channel, depolarising_channel
 from graphix.fundamentals import Plane
 from graphix.ops import Ops
@@ -917,7 +918,7 @@ def test_measure(self, pr_calc) -> None:
         pattern = circ.transpile().pattern
 
         measure_method = DefaultMeasureMethod()
-        backend = DensityMatrixBackend(pr_calc=pr_calc)
+        backend = DensityMatrixBackend(branch_selector=RandomBranchSelector(pr_calc=pr_calc))
         backend.add_nodes(pattern.input_nodes)
         backend.add_nodes([1, 2])
         backend.entangle_nodes((0, 1))
diff --git a/tests/test_pattern.py b/tests/test_pattern.py
index 58650824a..550f33736 100644
--- a/tests/test_pattern.py
+++ b/tests/test_pattern.py
@@ -10,6 +10,7 @@
 import pytest
 from numpy.random import PCG64, Generator
 
+from graphix.branch_selector import ConstBranchSelector, FixedBranchSelector
 from graphix.clifford import Clifford
 from graphix.command import C, CommandKind, E, M, N, X, Z
 from graphix.fundamentals import Plane
@@ -22,7 +23,6 @@
 from graphix.transpiler import Circuit
 
 if TYPE_CHECKING:
-    import collections.abc
     from collections.abc import Sequence
 
     from graphix.sim.base_backend import Backend
@@ -37,17 +37,6 @@ def compare_backend_result_with_statevec(backend: str, backend_state, statevec:
     raise NotImplementedError(backend)
 
 
-Outcome = typing.Literal[0, 1]
-
-
-class IterGenerator:
-    def __init__(self, it: collections.abc.Iterable[Outcome]) -> None:
-        self.__it = iter(it)
-
-    def choice(self, _outcomes: list[Outcome]) -> Outcome:
-        return next(self.__it)
-
-
 class TestPattern:
     def test_manual_generation(self) -> None:
         pattern = Pattern()
@@ -222,7 +211,7 @@ def test_pauli_measurement_single(self, plane: Plane, angle: float) -> None:
         pattern_ref = pattern.copy()
         pattern.perform_pauli_measurements()
         state = pattern.simulate_pattern()
-        state_ref = pattern_ref.simulate_pattern(pr_calc=False, rng=IterGenerator([0]))
+        state_ref = pattern_ref.simulate_pattern(branch_selector=ConstBranchSelector(False))
         assert np.abs(np.dot(state.flatten().conjugate(), state_ref.flatten())) == pytest.approx(1)
 
     @pytest.mark.parametrize("jumps", range(1, 11))
@@ -341,12 +330,13 @@ def test_shift_signals_plane(self, plane: Plane, method: str) -> None:
         pattern.standardize(method="mc")
         signal_dict = pattern.shift_signals(method=method)
         # Test for every possible outcome of each measure
-        for outcomes_ref in itertools.product(*([[0, 1]] * 3)):
-            state_ref = pattern_ref.simulate_pattern(pr_calc=False, rng=IterGenerator(iter(outcomes_ref)))
-            outcomes_p = shift_outcomes(dict(enumerate(outcomes_ref)), signal_dict)
-            state_p = pattern.simulate_pattern(
-                pr_calc=False, rng=IterGenerator(outcomes_p[i] for i in range(len(outcomes_p)))
-            )
+        for outcomes_ref_list in itertools.product(*([[0, 1]] * 3)):
+            outcomes_ref = dict(enumerate(outcomes_ref_list))
+            branch_selector = FixedBranchSelector(results=outcomes_ref)
+            state_ref = pattern_ref.simulate_pattern(branch_selector=branch_selector)
+            outcomes_p = shift_outcomes(outcomes_ref, signal_dict)
+            branch_selector = FixedBranchSelector(results=outcomes_p)
+            state_p = pattern.simulate_pattern(branch_selector=branch_selector)
             assert np.abs(np.dot(state_p.flatten().conjugate(), state_ref.flatten())) == pytest.approx(1)
 
     @pytest.mark.parametrize("jumps", range(1, 11))

From 8c0e4957ba8155f9643b44e41258dfd9b6b87ffb Mon Sep 17 00:00:00 2001
From: Thierry Martinez <thierry.martinez@inria.fr>
Date: Tue, 10 Jun 2025 15:49:06 +0200
Subject: [PATCH 02/19] Update PR number

---
 CHANGELOG.md | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index b2dc7e1ea..25f24622b 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -12,7 +12,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - #277: Methods for pretty-printing `Pattern`: `to_ascii`,
   `to_unicode`, `to_latex`.
 
-- #299: Branch selection in simulation: in addition to
+- #300: Branch selection in simulation: in addition to
   `RandomBranchSelector` which corresponds to the strategy that was
   already implemented, the user can use `FixedBranchSelector`,
   `ConstBranchSelector`, or define a custom branch selection by
@@ -33,7 +33,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - #277: The method `Pattern.print_pattern` is now deprecated.
 
-- #299: `pr_calc` parameter is deprecated in back-end initializers.
+- #300: `pr_calc` parameter is deprecated in back-end initializers.
   The user can specify `pr_calc` in the constructor of
   `RandomBranchSelector` instead.
 

From 22aae84d020dd6be86f3b54f8f0bc1a2010cf3d5 Mon Sep 17 00:00:00 2001
From: Thierry Martinez <thierry.martinez@inria.fr>
Date: Tue, 10 Jun 2025 17:59:03 +0200
Subject: [PATCH 03/19] Fix documentation

---
 docs/source/simulator.rst  | 24 ------------------
 graphix/branch_selector.py | 50 +++++++++++++++++++++-----------------
 graphix/transpiler.py      |  4 +--
 3 files changed, 30 insertions(+), 48 deletions(-)

diff --git a/docs/source/simulator.rst b/docs/source/simulator.rst
index fd4d6c97f..78b103dc5 100644
--- a/docs/source/simulator.rst
+++ b/docs/source/simulator.rst
@@ -12,30 +12,6 @@ Pattern Simulation
 
     .. automethod:: run
 
-Branch Selection (:mod:`graphix.branch_selector` module)
-++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-.. currentmodule:: graphix.branch_selector
-
-.. autoclass:: BranchSelector
-
-    .. automethod:: measure
-
-.. autoclass:: RandomBranchSelector
-    :members:
-
-    .. automethod:: measure
-
-.. autoclass:: FixedBranchSelector
-    :members:
-
-    .. automethod:: measure
-
-.. autoclass:: ConstBranchSelector
-    :members:
-
-    .. automethod:: measure
-
 Simulator backends
 ++++++++++++++++++
 
diff --git a/graphix/branch_selector.py b/graphix/branch_selector.py
index 49b130ece..cf301b452 100644
--- a/graphix/branch_selector.py
+++ b/graphix/branch_selector.py
@@ -32,13 +32,19 @@ class BranchSelector(ABC):
 
     @abstractmethod
     def measure(self, qubit: int, compute_expectation_0: Callable[[], float]) -> bool:
-        """
-        Return the measurement outcome of `qubit`.
+        """Return the measurement outcome of ``qubit``.
+
+
+        Parameters
+        ----------
+        qubit : int
+            Index of qubit to measure
 
-        This method may use the function `compute_expectation_0` to
-        retrieve the expected probability of outcome 0. The probability is
-        computed only if this function is called (lazy computation),
-        ensuring no unnecessary computational cost.
+        compute_expectation_0 : Callable[[], float]
+            A function that the method can use to retrieve the expected
+            probability of outcome 0. The probability is computed only if
+            this function is called (lazy computation), ensuring no
+            unnecessary computational cost.
         """
 
 
@@ -50,12 +56,12 @@ class RandomBranchSelector(BranchSelector):
     ----------
     pr_calc : bool, optional
         Whether to compute the probability distribution before selecting the measurement result.
-        If False, measurements yield 0/1 with equal probability (50% each).
-        Default is `True`.
+        If ``False``, measurements yield 0/1 with equal probability (50% each).
+        Default is ``True``.
     rng : Generator | None, optional
         Random-number generator for measurements.
-        If `None`, a default random-number generator is used.
-        Default is `None`.
+        If ``None``, a default random-number generator is used.
+        Default is ``None``.
     """
 
     pr_calc: bool = True
@@ -64,9 +70,9 @@ class RandomBranchSelector(BranchSelector):
     @override
     def measure(self, qubit: int, compute_expectation_0: Callable[[], float]) -> bool:
         """
-        Return the measurement outcome of `qubit`.
+        Return the measurement outcome of ``qubit``.
 
-        If `pr_calc` is `True`, the measurement outcome is determined based on the
+        If ``pr_calc`` is ``True``, the measurement outcome is determined based on the
         computed probability of outcome 0. Otherwise, the result is randomly chosen
         with a 50% chance for either outcome.
         """
@@ -82,20 +88,20 @@ def measure(self, qubit: int, compute_expectation_0: Callable[[], float]) -> boo
 class FixedBranchSelector(BranchSelector):
     """Branch selector with predefined measurement outcomes.
 
-    The mapping is fixed in `results`. By default, an error is raised if
+    The mapping is fixed in ``results``. By default, an error is raised if
     a qubit is measured without a predefined outcome. However, another
-    branch selector can be specified in `default` to handle such cases.
+    branch selector can be specified in ``default`` to handle such cases.
 
     Parameters
     ----------
     results : Mapping[int, bool]
         A dictionary mapping qubits to their measurement outcomes.
-        If a qubit is not present in this mapping, the `default` branch
+        If a qubit is not present in this mapping, the ``default`` branch
         selector is used.
     default : BranchSelector | None, optional
-        Branch selector to use for qubits not present in `results`.
-        If `None`, an error is raised when an unmapped qubit is measured.
-        Default is `None`.
+        Branch selector to use for qubits not present in ``results``.
+        If ``None``, an error is raised when an unmapped qubit is measured.
+        Default is ``None``.
     """
 
     results: Mapping[int, bool]
@@ -104,9 +110,9 @@ class FixedBranchSelector(BranchSelector):
     @override
     def measure(self, qubit: int, compute_expectation_0: Callable[[], float]) -> bool:
         """
-        Return the predefined measurement outcome of `qubit`, if available.
+        Return the predefined measurement outcome of ``qubit``, if available.
 
-        If the qubit is not present in `results`, the `default` branch selector
+        If the qubit is not present in ``results``, the ``default`` branch selector
         is used. If no default is provided, an error is raised.
         """
         result = self.results.get(qubit)
@@ -121,7 +127,7 @@ def measure(self, qubit: int, compute_expectation_0: Callable[[], float]) -> boo
 class ConstBranchSelector(BranchSelector):
     """Branch selector with a constant measurement outcome.
 
-    The value `result` is returned for every qubit.
+    The value ``result`` is returned for every qubit.
 
     Parameters
     ----------
@@ -133,5 +139,5 @@ class ConstBranchSelector(BranchSelector):
 
     @override
     def measure(self, qubit: int, compute_expectation_0: Callable[[], float]) -> bool:
-        """Return the constant measurement outcome `result` for any qubit."""
+        """Return the constant measurement outcome ``result`` for any qubit."""
         return self.result
diff --git a/graphix/transpiler.py b/graphix/transpiler.py
index 63624e97e..e0432da77 100644
--- a/graphix/transpiler.py
+++ b/graphix/transpiler.py
@@ -900,8 +900,8 @@ def simulate_statevector(
         ----------
         input_state : :class:`graphix.sim.statevec.Statevec`
 
-        branch_selector: :class:`graphix.sim.base_backend.BranchSelector`
-            branch selector for measures (default: `graphix.sim.base_backend.RandomBranchSelector()`)
+        branch_selector: :class:`graphix.branch_selector.BranchSelector`
+            branch selector for measures (default: `graphix.branch_selector.RandomBranchSelector()`)
 
         rng: :class:`np.random.Generator`
             random number generator for `RandomBranchSelector` (should only be used with default branch selector)

From 43ffe627d4aafc827374aa5ec5980c55f81c2e29 Mon Sep 17 00:00:00 2001
From: Thierry Martinez <thierry.martinez@inria.fr>
Date: Wed, 11 Jun 2025 21:45:32 +0200
Subject: [PATCH 04/19] ruff fix

---
 graphix/branch_selector.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/graphix/branch_selector.py b/graphix/branch_selector.py
index cf301b452..b5a79e92d 100644
--- a/graphix/branch_selector.py
+++ b/graphix/branch_selector.py
@@ -34,7 +34,6 @@ class BranchSelector(ABC):
     def measure(self, qubit: int, compute_expectation_0: Callable[[], float]) -> bool:
         """Return the measurement outcome of ``qubit``.
 
-
         Parameters
         ----------
         qubit : int

From a7e4609520674b0e13391340769e9f8418ab99c7 Mon Sep 17 00:00:00 2001
From: Thierry Martinez <thierry.martinez@inria.fr>
Date: Wed, 11 Jun 2025 21:50:24 +0200
Subject: [PATCH 05/19] Fix documentation

---
 graphix/transpiler.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/graphix/transpiler.py b/graphix/transpiler.py
index e0432da77..e4ade9868 100644
--- a/graphix/transpiler.py
+++ b/graphix/transpiler.py
@@ -901,10 +901,10 @@ def simulate_statevector(
         input_state : :class:`graphix.sim.statevec.Statevec`
 
         branch_selector: :class:`graphix.branch_selector.BranchSelector`
-            branch selector for measures (default: `graphix.branch_selector.RandomBranchSelector()`)
+            branch selector for measures (default: :class:`graphix.branch_selector.RandomBranchSelector`)
 
         rng: :class:`np.random.Generator`
-            random number generator for `RandomBranchSelector` (should only be used with default branch selector)
+            random number generator for :class:`graphix.branch_selector.RandomBranchSelector` (should only be used with default branch selector)
 
         Returns
         -------

From 3a4fb46da3e5a02e28fff87371c1e02fe2e467ee Mon Sep 17 00:00:00 2001
From: Thierry Martinez <thierry.martinez@inria.fr>
Date: Fri, 25 Jul 2025 13:37:46 +0200
Subject: [PATCH 06/19] Add branch selector support for tensor network backend

---
 graphix/sim/base_backend.py |  7 ++---
 graphix/sim/tensornet.py    | 37 ++++++++++++++-------------
 graphix/simulator.py        | 51 +++++++++++++++++++------------------
 tests/test_tnsim.py         |  5 ++--
 4 files changed, 51 insertions(+), 49 deletions(-)

diff --git a/graphix/sim/base_backend.py b/graphix/sim/base_backend.py
index af3246bd5..d23d927e7 100644
--- a/graphix/sim/base_backend.py
+++ b/graphix/sim/base_backend.py
@@ -742,11 +742,8 @@ class DenseStateBackend(Backend[_DenseStateT_co], Generic[_DenseStateT_co]):
     ----------
     node_index : NodeIndex, optional
         Mapping between node numbers and qubit indices in the internal state of the backend.
-    pr_calc : bool, optional
-        Whether or not to compute the probability distribution before choosing the measurement outcome.
-        If False, measurements yield results 0/1 with 50% probabilities each.
-    rng : Generator, optional
-        Random number generator used to sample measurement outcomes.
+    branch_selector: :class:`graphix.branch_selector.BranchSelector`, optional
+        Branch selector used for measurements.  Default is :class:`RandomBranchSelector`.
     symbolic : bool, optional
         If True, support arbitrary objects (typically, symbolic expressions) in matrices.
 
diff --git a/graphix/sim/tensornet.py b/graphix/sim/tensornet.py
index fb6284dfe..77a89cd38 100644
--- a/graphix/sim/tensornet.py
+++ b/graphix/sim/tensornet.py
@@ -20,9 +20,9 @@
 from typing_extensions import TypeAlias, override
 
 from graphix import command
+from graphix.branch_selector import BranchSelector, RandomBranchSelector
 from graphix.ops import Ops
 from graphix.parameter import Expression
-from graphix.rng import ensure_rng
 from graphix.sim.base_backend import Backend, BackendState
 from graphix.states import BasicStates, PlanarState
 
@@ -30,7 +30,6 @@
     from collections.abc import Iterable, Sequence
 
     from cotengra.oe import PathOptimizer
-    from numpy.random import Generator
 
     from graphix import Pattern
     from graphix.clifford import Clifford
@@ -53,7 +52,7 @@ class MBQCTensorNet(BackendState, TensorNetwork):
 
     def __init__(
         self,
-        rng: Generator | None = None,
+        branch_selector: BranchSelector,
         graph_nodes: Iterable[int] | None = None,
         graph_edges: Iterable[tuple[int, int]] | None = None,
         default_output_nodes: Iterable[int] | None = None,
@@ -82,7 +81,7 @@ def __init__(
         # prepare the graph state if graph_nodes and graph_edges are given
         if graph_nodes is not None and graph_edges is not None:
             self.set_graph_state(graph_nodes, graph_edges)
-        self.__rng = ensure_rng(rng)
+        self.__branch_selector = branch_selector
 
     def get_open_tensor_from_index(self, index: int | str) -> npt.NDArray[np.complex128]:
         """Get tensor specified by node index. The tensor has a dangling edge.
@@ -231,7 +230,7 @@ def measure_single(
             measurement result.
         """
         if bypass_probability_calculation:
-            result = outcome if outcome is not None else self.__rng.choice([0, 1])
+            result = outcome if outcome is not None else self.__branch_selector.measure(index, lambda: 0.5)
             # Basis state to be projected
             if isinstance(basis, np.ndarray):
                 if outcome is not None:
@@ -537,7 +536,7 @@ def copy(self, virtual: bool = False, deep: bool = False) -> MBQCTensorNet:
         """
         if deep:
             return deepcopy(self)
-        return self.__class__(rng=self.__rng, ts=self)
+        return self.__class__(branch_selector=self.__branch_selector, ts=self)
 
 
 def _get_decomposed_cz() -> list[npt.NDArray[np.complex128]]:
@@ -581,7 +580,7 @@ class _AbstractTensorNetworkBackend(Backend[MBQCTensorNet], ABC):
     pattern: Pattern
     graph_prep: str
     input_state: Data
-    rng: Generator
+    branch_selector: BranchSelector
     output_nodes: list[int]
     results: dict[int, Outcome]
     _decomposed_cz: list[npt.NDArray[np.complex128]]
@@ -609,12 +608,16 @@ class TensorNetworkBackend(_AbstractTensorNetworkBackend):
         'auto'(default) :
             Automatically select a preparation strategy based on the max degree of a graph
     input_state : preparation for input states (only BasicStates.PLUS is supported for tensor networks yet),
-    rng: :class:`np.random.Generator` (default: *None*)
-        random number generator to use for measurements
+    branch_selector: :class:`graphix.branch_selector.BranchSelector`, optional
+        Branch selector to be used for measurements.
     """
 
     def __init__(
-        self, pattern: Pattern, graph_prep: str = "auto", input_state: Data | None = None, rng: Generator | None = None
+        self,
+        pattern: Pattern,
+        graph_prep: str = "auto",
+        input_state: Data | None = None,
+        branch_selector: BranchSelector | None = None,
     ) -> None:
         """Construct a tensor network backend."""
         if input_state is None:
@@ -622,7 +625,8 @@ def __init__(
         elif input_state != BasicStates.PLUS:
             msg = "TensorNetworkBackend currently only supports BasicStates.PLUS as input state."
             raise NotImplementedError(msg)
-        rng = ensure_rng(rng)
+        if branch_selector is None:
+            branch_selector = RandomBranchSelector()
         if graph_prep in {"parallel", "sequential"}:
             pass
         elif graph_prep == "opt":
@@ -646,11 +650,11 @@ def __init__(
                 graph_nodes=nodes,
                 graph_edges=edges,
                 default_output_nodes=pattern.output_nodes,
-                rng=rng,
+                branch_selector=branch_selector,
             )
             decomposed_cz = []
         else:  # graph_prep == "sequential":
-            state = MBQCTensorNet(default_output_nodes=pattern.output_nodes, rng=rng)
+            state = MBQCTensorNet(default_output_nodes=pattern.output_nodes, branch_selector=branch_selector)
             decomposed_cz = _get_decomposed_cz()
         isolated_nodes = pattern.get_isolated_nodes()
         super().__init__(
@@ -658,7 +662,7 @@ def __init__(
             pattern,
             graph_prep,
             input_state,
-            rng,
+            branch_selector,
             pattern.output_nodes,
             results,
             decomposed_cz,
@@ -752,12 +756,11 @@ def measure(self, node: int, measurement: Measurement) -> Outcome:
             vector: npt.NDArray[np.complex128] = self.state.get_open_tensor_from_index(node)
             probs = (np.abs(vector) ** 2).astype(np.float64)
             probs /= np.sum(probs)
-            result: Outcome = self.rng.choice([0, 1], p=probs)
+            result: Outcome = self.branch_selector.measure(node, lambda: probs[0])
             self.results[node] = result
             buffer = 1 / probs[result] ** 0.5
         else:
-            # choose the measurement result randomly
-            result = self.rng.choice([0, 1])
+            result = self.branch_selector.measure(node, lambda: 0.5)
             self.results[node] = result
             buffer = 2**0.5
         if isinstance(measurement.angle, Expression):
diff --git a/graphix/simulator.py b/graphix/simulator.py
index 6096dd40d..aa39b747c 100644
--- a/graphix/simulator.py
+++ b/graphix/simulator.py
@@ -202,7 +202,7 @@ def __init__(
         noise_model: :class:`graphix.noise_models.noise_model.NoiseModel`, optional
             [Density matrix backend only] Noise model used by the simulator.
         branch_selector: :class:`graphix.branch_selector.BranchSelector`, optional
-            Branch selector used by the backend. Can only be specified if `backend` is not an already instantiated :class:`graphix.sim.backend.Backend` object.  Default is :class:`RandomBranchSelector`.
+            Branch selector used for measurements. Can only be specified if `backend` is not an already instantiated :class:`graphix.sim.backend.Backend` object.  Default is :class:`RandomBranchSelector`.
         rng: :class:`numpy.random.Generator`, optional
             Random number generator to be used by the default `RandomBranchSelector`. Can only be specified if `backend` is not an already instantiated :class:`graphix.sim.backend.Backend` object and if `branch_selector` is not specified.
         graph_prep: str, optional
@@ -211,44 +211,45 @@ def __init__(
             :class:`graphix.sim.tensornet.TensorNetworkBackend`\
             :class:`graphix.sim.density_matrix.DensityMatrixBackend`\
         """
-        if isinstance(backend, Backend):
-            if noise_model is not None:
-                raise ValueError("`noise_model` cannot be specified if `backend` is already instantiated.")
-            if branch_selector is not None:
-                raise ValueError("`branch_selector` cannot be specified if `backend` is already instantiated.")
-            if rng is not None:
-                raise ValueError("`rng` cannot be specified if `backend` is already instantiated.")
-            if graph_prep is not None:
-                raise ValueError("`graph_prep` cannot be specified if `backend` is already instantiated.")
-            self.backend = backend
-        elif backend in {"tensornetwork", "mps"}:
-            if noise_model is not None:
-                raise ValueError("`noise_model` cannot be specified for tensor network backend.")
-            if branch_selector is not None:
-                raise ValueError("`branch_selector` cannot be specified for tensor network backend.")
-            if graph_prep is None:
-                graph_prep = "auto"
-            self.backend = TensorNetworkBackend(pattern, rng=rng, graph_prep=graph_prep)
-        else:
+
+        def initialize_backend() -> Backend[BackendState]:
+            nonlocal backend, branch_selector, rng, graph_prep, noise_model
+            if isinstance(backend, Backend):
+                if noise_model is not None:
+                    raise ValueError("`noise_model` cannot be specified if `backend` is already instantiated.")
+                if branch_selector is not None:
+                    raise ValueError("`branch_selector` cannot be specified if `backend` is already instantiated.")
+                if rng is not None:
+                    raise ValueError("`rng` cannot be specified if `backend` is already instantiated.")
+                if graph_prep is not None:
+                    raise ValueError("`graph_prep` cannot be specified if `backend` is already instantiated.")
+                return backend
             if branch_selector is None:
                 branch_selector = RandomBranchSelector(rng=rng)
             elif rng is not None:
                 raise ValueError("`rng` and `branch_selector` cannot be specified simultaneously.")
+            if backend in {"tensornetwork", "mps"}:
+                if noise_model is not None:
+                    raise ValueError("`noise_model` cannot be specified for tensor network backend.")
+                if graph_prep is None:
+                    graph_prep = "auto"
+                return TensorNetworkBackend(pattern, branch_selector=branch_selector, graph_prep=graph_prep)
             if graph_prep is not None:
                 raise ValueError("`graph_prep` can only be specified for tensor network backend.")
             if backend == "statevector":
                 if noise_model is not None:
                     raise ValueError("`noise_model` cannot be specified for state vector backend.")
-                self.backend = StatevectorBackend(branch_selector=branch_selector)
-            elif backend == "densitymatrix":
+                return StatevectorBackend(branch_selector=branch_selector)
+            if backend == "densitymatrix":
                 if noise_model is None:
                     warnings.warn(
                         "Simulating using densitymatrix backend with no noise. To add noise to the simulation, give an object of `graphix.noise_models.Noisemodel` to `noise_model` keyword argument.",
                         stacklevel=1,
                     )
-                self.backend = DensityMatrixBackend(branch_selector=branch_selector)
-            else:
-                raise ValueError(f"Unknown backend {backend}.")
+                return DensityMatrixBackend(branch_selector=branch_selector)
+            raise ValueError(f"Unknown backend {backend}.")
+
+        self.backend = initialize_backend()
         self.noise_model = noise_model
         self.__pattern = pattern
         if measure_method is None:
diff --git a/tests/test_tnsim.py b/tests/test_tnsim.py
index 0410f7f71..14f6daa1e 100644
--- a/tests/test_tnsim.py
+++ b/tests/test_tnsim.py
@@ -8,6 +8,7 @@
 from numpy.random import PCG64, Generator
 from quimb.tensor import Tensor
 
+from graphix.branch_selector import RandomBranchSelector
 from graphix.clifford import Clifford
 from graphix.command import C, E, X, Z
 from graphix.ops import Ops
@@ -33,7 +34,7 @@ def random_op(sites: int, dtype: type, rng: Generator) -> npt.NDArray:
 class TestTN:
     def test_add_node(self, fx_rng: Generator) -> None:
         node_index = fx_rng.integers(0, 1000)
-        tn = MBQCTensorNet(rng=fx_rng)
+        tn = MBQCTensorNet(branch_selector=RandomBranchSelector(rng=fx_rng))
 
         tn.add_qubit(node_index)
 
@@ -42,7 +43,7 @@ def test_add_node(self, fx_rng: Generator) -> None:
 
     def test_add_nodes(self, fx_rng: Generator) -> None:
         node_index = set(fx_rng.integers(0, 1000, 20))
-        tn = MBQCTensorNet(rng=fx_rng)
+        tn = MBQCTensorNet(branch_selector=RandomBranchSelector(rng=fx_rng))
 
         tn.graph_prep = "sequential"
         tn.add_qubits(node_index)

From 285b0778b34e02b1e8265f13d2106606f924f9b2 Mon Sep 17 00:00:00 2001
From: Thierry Martinez <thierry.martinez@inria.fr>
Date: Fri, 25 Jul 2025 16:05:52 +0200
Subject: [PATCH 07/19] Add specific tests for branch selectors

---
 tests/test_branch_selector.py | 162 ++++++++++++++++++++++++++++++++++
 1 file changed, 162 insertions(+)
 create mode 100644 tests/test_branch_selector.py

diff --git a/tests/test_branch_selector.py b/tests/test_branch_selector.py
new file mode 100644
index 000000000..e91ecc89a
--- /dev/null
+++ b/tests/test_branch_selector.py
@@ -0,0 +1,162 @@
+from __future__ import annotations
+
+import dataclasses
+import itertools
+import math
+from dataclasses import dataclass
+from typing import TYPE_CHECKING, Callable
+
+import pytest
+from typing_extensions import override
+
+from graphix import Pattern
+from graphix.branch_selector import ConstBranchSelector, FixedBranchSelector, RandomBranchSelector
+from graphix.command import M, N
+from graphix.simulator import DefaultMeasureMethod
+
+if TYPE_CHECKING:
+    from collections.abc import Mapping
+
+    from numpy.random import Generator
+
+    from graphix.measurements import Outcome
+
+NB_ROUNDS = 100
+
+
+@dataclass
+class CheckedBranchSelector(RandomBranchSelector):
+    """Random branch selector that verifies that expectation values match the expected ones."""
+
+    rng: Generator | None = None
+    expected: Mapping[int, float] = dataclasses.field(default_factory=dict)
+
+    @override
+    def measure(self, qubit: int, compute_expectation_0: Callable[[], float]) -> Outcome:
+        """Return the measurement outcome of ``qubit``."""
+        expectation_0 = compute_expectation_0()
+        assert math.isclose(expectation_0, self.expected[qubit])
+        return super().measure(qubit, lambda: expectation_0)
+
+
+@pytest.mark.filterwarnings("ignore:Simulating using densitymatrix backend with no noise.")
+@pytest.mark.parametrize(
+    "backend",
+    [
+        "statevector",
+        "densitymatrix",
+        pytest.param(
+            "tensornetwork",
+            marks=pytest.mark.xfail(
+                reason="[Bug]: TensorNetworkBackend computes incorrect measurement probabilities #325"
+            ),
+        ),
+    ],
+)
+def test_expectation_value(fx_rng: Generator, backend: str) -> None:
+    # Pattern that measures 0 on qubit 0 with probability 1.
+    pattern = Pattern(cmds=[N(0), M(0)])
+    branch_selector = CheckedBranchSelector(rng=fx_rng, expected={0: 1.0})
+    pattern.simulate_pattern(backend, branch_selector=branch_selector)
+
+
+@pytest.mark.filterwarnings("ignore:Simulating using densitymatrix backend with no noise.")
+@pytest.mark.parametrize(
+    "backend",
+    [
+        "statevector",
+        "densitymatrix",
+        pytest.param(
+            "tensornetwork",
+            marks=pytest.mark.xfail(
+                reason="[Bug]: TensorNetworkBackend computes incorrect measurement probabilities #325"
+            ),
+        ),
+    ],
+)
+def test_random_branch_selector(fx_rng: Generator, backend: str) -> None:
+    branch_selector = RandomBranchSelector(rng=fx_rng)
+    pattern = Pattern(cmds=[N(0), M(0)])
+    for _ in range(NB_ROUNDS):
+        measure_method = DefaultMeasureMethod()
+        pattern.simulate_pattern(backend, branch_selector=branch_selector, measure_method=measure_method)
+        assert measure_method.results[0] == 0
+
+
+@pytest.mark.filterwarnings("ignore:Simulating using densitymatrix backend with no noise.")
+@pytest.mark.parametrize(
+    "backend",
+    [
+        "statevector",
+        "densitymatrix",
+        "tensornetwork",
+    ],
+)
+def test_random_branch_selector_without_pr_calc(backend: str) -> None:
+    branch_selector = RandomBranchSelector(pr_calc=False)
+    # Pattern that measures 0 on qubit 0 with probability > 0.999999999, to avoid numerical errors when exploring impossible branches.
+    pattern = Pattern(cmds=[N(0), M(0, angle=1e-5)])
+    nb_outcome_1 = 0
+    for _ in range(NB_ROUNDS):
+        measure_method = DefaultMeasureMethod()
+        pattern.simulate_pattern(backend, branch_selector=branch_selector, measure_method=measure_method)
+        if measure_method.results[0]:
+            nb_outcome_1 += 1
+    assert abs(nb_outcome_1 - NB_ROUNDS / 2) < NB_ROUNDS / 10
+
+
+@pytest.mark.filterwarnings("ignore:Simulating using densitymatrix backend with no noise.")
+@pytest.mark.parametrize(
+    "backend",
+    [
+        "statevector",
+        "densitymatrix",
+        "tensornetwork",
+    ],
+)
+@pytest.mark.parametrize("outcome", itertools.product([0, 1], repeat=3))
+def test_fixed_branch_selector(backend: str, outcome: list[Outcome]) -> None:
+    branch_selector = FixedBranchSelector(
+        results=dict(enumerate(outcome[:-1])), default=FixedBranchSelector({2: outcome[2]})
+    )
+    pattern = Pattern(cmds=[cmd for qubit in range(3) for cmd in (N(qubit), M(qubit, angle=0.1))])
+    measure_method = DefaultMeasureMethod()
+    pattern.simulate_pattern(backend, branch_selector=branch_selector, measure_method=measure_method)
+    for qubit, value in enumerate(outcome):
+        assert measure_method.results[qubit] == value
+
+
+@pytest.mark.filterwarnings("ignore:Simulating using densitymatrix backend with no noise.")
+@pytest.mark.parametrize(
+    "backend",
+    [
+        "statevector",
+        "densitymatrix",
+        "tensornetwork",
+    ],
+)
+def test_fixed_branch_selector_no_default(backend: str) -> None:
+    branch_selector = FixedBranchSelector(results={})
+    pattern = Pattern(cmds=[N(0), M(0, angle=1e-5)])
+    measure_method = DefaultMeasureMethod()
+    with pytest.raises(ValueError):
+        pattern.simulate_pattern(backend, branch_selector=branch_selector, measure_method=measure_method)
+
+
+@pytest.mark.filterwarnings("ignore:Simulating using densitymatrix backend with no noise.")
+@pytest.mark.parametrize(
+    "backend",
+    [
+        "statevector",
+        "densitymatrix",
+        "tensornetwork",
+    ],
+)
+@pytest.mark.parametrize("outcome", [0, 1])
+def test_const_branch_selector(backend: str, outcome: Outcome) -> None:
+    branch_selector = ConstBranchSelector(outcome)
+    pattern = Pattern(cmds=[N(0), M(0, angle=1e-5)])
+    for _ in range(NB_ROUNDS):
+        measure_method = DefaultMeasureMethod()
+        pattern.simulate_pattern(backend, branch_selector=branch_selector, measure_method=measure_method)
+        assert measure_method.results[0] == outcome

From 5bafa13f832e37cbf26f4e16a22e3dd267896a0f Mon Sep 17 00:00:00 2001
From: Thierry Martinez <thierry.martinez@inria.fr>
Date: Fri, 25 Jul 2025 16:29:57 +0200
Subject: [PATCH 08/19] Fix symbolic

---
 CHANGELOG.md         |  2 +-
 graphix/simulator.py | 12 ++++++++++--
 noxfile.py           |  3 ++-
 3 files changed, 13 insertions(+), 4 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 8a4ee4917..542796872 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -67,7 +67,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - #277: The method `Pattern.print_pattern` is now deprecated.
 
-- #300: `pr_calc` parameter is deprecated in back-end initializers.
+- #300: `pr_calc` parameter is removed in back-end initializers.
   The user can specify `pr_calc` in the constructor of
   `RandomBranchSelector` instead.
 
diff --git a/graphix/simulator.py b/graphix/simulator.py
index aa39b747c..c31359e44 100644
--- a/graphix/simulator.py
+++ b/graphix/simulator.py
@@ -186,6 +186,7 @@ def __init__(
         branch_selector: BranchSelector | None = None,
         rng: Generator | None = None,
         graph_prep: str | None = None,
+        symbolic: bool = False,
     ) -> None:
         """
         Construct a pattern simulator.
@@ -207,6 +208,9 @@ def __init__(
             Random number generator to be used by the default `RandomBranchSelector`. Can only be specified if `backend` is not an already instantiated :class:`graphix.sim.backend.Backend` object and if `branch_selector` is not specified.
         graph_prep: str, optional
             [Tensor network backend only] Strategy for preparing the graph state.  See :class:`graphix.sim.tensornet.TensorNetworkBackend`.
+        symbolic : bool, optional
+            [State vector and density matrix backends only] If True, support arbitrary objects (typically, symbolic expressions) in measurement angles.
+
         .. seealso:: :class:`graphix.sim.statevec.StatevectorBackend`\
             :class:`graphix.sim.tensornet.TensorNetworkBackend`\
             :class:`graphix.sim.density_matrix.DensityMatrixBackend`\
@@ -223,6 +227,8 @@ def initialize_backend() -> Backend[BackendState]:
                     raise ValueError("`rng` cannot be specified if `backend` is already instantiated.")
                 if graph_prep is not None:
                     raise ValueError("`graph_prep` cannot be specified if `backend` is already instantiated.")
+                if symbolic:
+                    raise ValueError("`symbolic` cannot be specified if `backend` is already instantiated.")
                 return backend
             if branch_selector is None:
                 branch_selector = RandomBranchSelector(rng=rng)
@@ -231,6 +237,8 @@ def initialize_backend() -> Backend[BackendState]:
             if backend in {"tensornetwork", "mps"}:
                 if noise_model is not None:
                     raise ValueError("`noise_model` cannot be specified for tensor network backend.")
+                if symbolic:
+                    raise ValueError("`symbolic` cannot be specified for tensor network backend.")
                 if graph_prep is None:
                     graph_prep = "auto"
                 return TensorNetworkBackend(pattern, branch_selector=branch_selector, graph_prep=graph_prep)
@@ -239,14 +247,14 @@ def initialize_backend() -> Backend[BackendState]:
             if backend == "statevector":
                 if noise_model is not None:
                     raise ValueError("`noise_model` cannot be specified for state vector backend.")
-                return StatevectorBackend(branch_selector=branch_selector)
+                return StatevectorBackend(branch_selector=branch_selector, symbolic=symbolic)
             if backend == "densitymatrix":
                 if noise_model is None:
                     warnings.warn(
                         "Simulating using densitymatrix backend with no noise. To add noise to the simulation, give an object of `graphix.noise_models.Noisemodel` to `noise_model` keyword argument.",
                         stacklevel=1,
                     )
-                return DensityMatrixBackend(branch_selector=branch_selector)
+                return DensityMatrixBackend(branch_selector=branch_selector, symbolic=symbolic)
             raise ValueError(f"Unknown backend {backend}.")
 
         self.backend = initialize_backend()
diff --git a/noxfile.py b/noxfile.py
index b64c9a64a..0d3c994ae 100644
--- a/noxfile.py
+++ b/noxfile.py
@@ -62,6 +62,7 @@ def tests_symbolic(session: Session) -> None:
     with TemporaryDirectory() as tmpdir, session.cd(tmpdir):
         # If you need a specific branch:
         # session.run("git", "clone", "-b", "branch-name", "https://github.com/TeamGraphix/graphix-symbolic")
-        session.run("git", "clone", "https://github.com/TeamGraphix/graphix-symbolic")
+        session.run("git", "clone", "-b", "branch_selector", "https://github.com/thierry-martinez/graphix-symbolic")
+        # session.run("git", "clone", "https://github.com/TeamGraphix/graphix-symbolic")
         with session.cd("graphix-symbolic"):
             session.run("pytest", "--doctest-modules")

From bddda3c81f19ed1bfb9ff49628510785620ab3d2 Mon Sep 17 00:00:00 2001
From: Thierry Martinez <thierry.martinez@inria.fr>
Date: Fri, 25 Jul 2025 16:31:20 +0200
Subject: [PATCH 09/19] Add comment for graphix-symbolic

---
 noxfile.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/noxfile.py b/noxfile.py
index 0d3c994ae..4f19886e6 100644
--- a/noxfile.py
+++ b/noxfile.py
@@ -62,6 +62,7 @@ def tests_symbolic(session: Session) -> None:
     with TemporaryDirectory() as tmpdir, session.cd(tmpdir):
         # If you need a specific branch:
         # session.run("git", "clone", "-b", "branch-name", "https://github.com/TeamGraphix/graphix-symbolic")
+        # See https://github.com/TeamGraphix/graphix-symbolic/pull/4
         session.run("git", "clone", "-b", "branch_selector", "https://github.com/thierry-martinez/graphix-symbolic")
         # session.run("git", "clone", "https://github.com/TeamGraphix/graphix-symbolic")
         with session.cd("graphix-symbolic"):

From 776edc637e957bab36a57f04b800c7a95ce1b6cd Mon Sep 17 00:00:00 2001
From: Thierry Martinez <thierry.martinez@inria.fr>
Date: Fri, 25 Jul 2025 16:40:27 +0200
Subject: [PATCH 10/19] Fix documentation

---
 graphix/simulator.py | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/graphix/simulator.py b/graphix/simulator.py
index c31359e44..4ed3153e3 100644
--- a/graphix/simulator.py
+++ b/graphix/simulator.py
@@ -193,21 +193,21 @@ def __init__(
 
         Parameters
         ----------
-        pattern: :class:`graphix.pattern.Pattern` object
+        pattern: :class:`Pattern` object
             MBQC pattern to be simulated.
-        backend: :class:`graphix.sim.backend.Backend` object,
+        backend: :class:`Backend` object,
             or 'statevector', or 'densitymatrix', or 'tensornetwork'
             simulation backend (optional), default is 'statevector'.
         measure_method: :class:`MeasureMethod`, optional
             Measure method used by the simulator. Default is :class:`DefaultMeasureMethod`.
-        noise_model: :class:`graphix.noise_models.noise_model.NoiseModel`, optional
+        noise_model: :class:`NoiseModel`, optional
             [Density matrix backend only] Noise model used by the simulator.
-        branch_selector: :class:`graphix.branch_selector.BranchSelector`, optional
-            Branch selector used for measurements. Can only be specified if `backend` is not an already instantiated :class:`graphix.sim.backend.Backend` object.  Default is :class:`RandomBranchSelector`.
+        branch_selector: :class:`BranchSelector`, optional
+            Branch selector used for measurements. Can only be specified if ``backend`` is not an already instantiated :class:`Backend` object.  Default is :class:`RandomBranchSelector`.
         rng: :class:`numpy.random.Generator`, optional
-            Random number generator to be used by the default `RandomBranchSelector`. Can only be specified if `backend` is not an already instantiated :class:`graphix.sim.backend.Backend` object and if `branch_selector` is not specified.
+            Random number generator to be used by the default :class:`RandomBranchSelector`. Can only be specified if ``backend`` is not an already instantiated :class:`Backend` object and if ``branch_selector`` is not specified.
         graph_prep: str, optional
-            [Tensor network backend only] Strategy for preparing the graph state.  See :class:`graphix.sim.tensornet.TensorNetworkBackend`.
+            [Tensor network backend only] Strategy for preparing the graph state.  See :class:`TensorNetworkBackend`.
         symbolic : bool, optional
             [State vector and density matrix backends only] If True, support arbitrary objects (typically, symbolic expressions) in measurement angles.
 

From 379b8dfd6df6791e7cadd0ff22ced353a3779ef5 Mon Sep 17 00:00:00 2001
From: Thierry Martinez <thierry.martinez@inria.fr>
Date: Fri, 25 Jul 2025 16:42:18 +0200
Subject: [PATCH 11/19] Relax test condition

---
 tests/test_branch_selector.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test_branch_selector.py b/tests/test_branch_selector.py
index e91ecc89a..1ad47bcd9 100644
--- a/tests/test_branch_selector.py
+++ b/tests/test_branch_selector.py
@@ -102,7 +102,7 @@ def test_random_branch_selector_without_pr_calc(backend: str) -> None:
         pattern.simulate_pattern(backend, branch_selector=branch_selector, measure_method=measure_method)
         if measure_method.results[0]:
             nb_outcome_1 += 1
-    assert abs(nb_outcome_1 - NB_ROUNDS / 2) < NB_ROUNDS / 10
+    assert abs(nb_outcome_1 - NB_ROUNDS / 2) < NB_ROUNDS / 5
 
 
 @pytest.mark.filterwarnings("ignore:Simulating using densitymatrix backend with no noise.")

From b651128d09b0dabaad3af20047a93d8c60c2d466 Mon Sep 17 00:00:00 2001
From: Thierry Martinez <thierry.martinez@inria.fr>
Date: Fri, 25 Jul 2025 16:59:02 +0200
Subject: [PATCH 12/19] Remove useless redefinition of the `rng` field

---
 tests/test_branch_selector.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/tests/test_branch_selector.py b/tests/test_branch_selector.py
index 1ad47bcd9..8ea539d8f 100644
--- a/tests/test_branch_selector.py
+++ b/tests/test_branch_selector.py
@@ -28,7 +28,6 @@
 class CheckedBranchSelector(RandomBranchSelector):
     """Random branch selector that verifies that expectation values match the expected ones."""
 
-    rng: Generator | None = None
     expected: Mapping[int, float] = dataclasses.field(default_factory=dict)
 
     @override

From 5471c0b9eacc29c87558014ee144af131c844142 Mon Sep 17 00:00:00 2001
From: Thierry Martinez <thierry.martinez@inria.fr>
Date: Mon, 28 Jul 2025 21:08:58 +0200
Subject: [PATCH 13/19] Rename `compute_expectation_0` into `f_expectation0`

---
 graphix/branch_selector.py    | 14 +++++++-------
 graphix/sim/base_backend.py   | 22 +++++++++++-----------
 tests/test_branch_selector.py |  8 ++++----
 3 files changed, 22 insertions(+), 22 deletions(-)

diff --git a/graphix/branch_selector.py b/graphix/branch_selector.py
index 6ffdb9a32..467636d07 100644
--- a/graphix/branch_selector.py
+++ b/graphix/branch_selector.py
@@ -34,7 +34,7 @@ class BranchSelector(ABC):
     """
 
     @abstractmethod
-    def measure(self, qubit: int, compute_expectation_0: Callable[[], float]) -> Outcome:
+    def measure(self, qubit: int, f_expectation0: Callable[[], float]) -> Outcome:
         """Return the measurement outcome of ``qubit``.
 
         Parameters
@@ -42,7 +42,7 @@ def measure(self, qubit: int, compute_expectation_0: Callable[[], float]) -> Out
         qubit : int
             Index of qubit to measure
 
-        compute_expectation_0 : Callable[[], float]
+        f_expectation0 : Callable[[], float]
             A function that the method can use to retrieve the expected
             probability of outcome 0. The probability is computed only if
             this function is called (lazy computation), ensuring no
@@ -70,7 +70,7 @@ class RandomBranchSelector(BranchSelector):
     rng: Generator | None = None
 
     @override
-    def measure(self, qubit: int, compute_expectation_0: Callable[[], float]) -> Outcome:
+    def measure(self, qubit: int, f_expectation0: Callable[[], float]) -> Outcome:
         """
         Return the measurement outcome of ``qubit``.
 
@@ -80,7 +80,7 @@ def measure(self, qubit: int, compute_expectation_0: Callable[[], float]) -> Out
         """
         self.rng = ensure_rng(self.rng)
         if self.pr_calc:
-            prob_0 = compute_expectation_0()
+            prob_0 = f_expectation0()
             return outcome(self.rng.random() > prob_0)
         result: Outcome = self.rng.choice([0, 1])
         return result
@@ -110,7 +110,7 @@ class FixedBranchSelector(BranchSelector):
     default: BranchSelector | None = None
 
     @override
-    def measure(self, qubit: int, compute_expectation_0: Callable[[], float]) -> Outcome:
+    def measure(self, qubit: int, f_expectation0: Callable[[], float]) -> Outcome:
         """
         Return the predefined measurement outcome of ``qubit``, if available.
 
@@ -121,7 +121,7 @@ def measure(self, qubit: int, compute_expectation_0: Callable[[], float]) -> Out
         if result is None:
             if self.default is None:
                 raise ValueError(f"Unexpected measurement of qubit {qubit}.")
-            return self.default.measure(qubit, compute_expectation_0)
+            return self.default.measure(qubit, f_expectation0)
         return result
 
 
@@ -140,6 +140,6 @@ class ConstBranchSelector(BranchSelector):
     result: Outcome
 
     @override
-    def measure(self, qubit: int, compute_expectation_0: Callable[[], float]) -> Outcome:
+    def measure(self, qubit: int, f_expectation0: Callable[[], float]) -> Outcome:
         """Return the constant measurement outcome ``result`` for any qubit."""
         return self.result
diff --git a/graphix/sim/base_backend.py b/graphix/sim/base_backend.py
index d23d927e7..3ce9ebd29 100644
--- a/graphix/sim/base_backend.py
+++ b/graphix/sim/base_backend.py
@@ -529,23 +529,23 @@ def perform_measure(
 ) -> Outcome:
     """Perform measurement of a qubit."""
     vec = plane.polar(angle)
-    # op_mat_0 may contain the matrix operator associated with the outcome 0,
+    # op_mat0 may contain the matrix operator associated with the outcome 0,
     # but the value is computed lazily, i.e., only if needed.
-    op_mat_0 = None
+    op_mat0 = None
 
-    def get_op_mat_0() -> Matrix:
-        nonlocal op_mat_0
-        if op_mat_0 is None:
-            op_mat_0 = _op_mat_from_result(vec, 0, symbolic=symbolic)
-        return op_mat_0
+    def get_op_mat0() -> Matrix:
+        nonlocal op_mat0
+        if op_mat0 is None:
+            op_mat0 = _op_mat_from_result(vec, 0, symbolic=symbolic)
+        return op_mat0
 
-    def compute_expectation_0() -> float:
-        exp_val = state.expectation_single(get_op_mat_0(), qubit_loc)
+    def f_expectation0() -> float:
+        exp_val = state.expectation_single(get_op_mat0(), qubit_loc)
         assert math.isclose(exp_val.imag, 0, abs_tol=1e-10)
         return exp_val.real
 
-    result = branch_selector.measure(qubit_node, compute_expectation_0)
-    op_mat = _op_mat_from_result(vec, 1, symbolic=symbolic) if result else get_op_mat_0()
+    result = branch_selector.measure(qubit_node, f_expectation0)
+    op_mat = _op_mat_from_result(vec, 1, symbolic=symbolic) if result else get_op_mat0()
     state.evolve_single(op_mat, qubit_loc)
     return result
 
diff --git a/tests/test_branch_selector.py b/tests/test_branch_selector.py
index 8ea539d8f..07e7d36a2 100644
--- a/tests/test_branch_selector.py
+++ b/tests/test_branch_selector.py
@@ -31,11 +31,11 @@ class CheckedBranchSelector(RandomBranchSelector):
     expected: Mapping[int, float] = dataclasses.field(default_factory=dict)
 
     @override
-    def measure(self, qubit: int, compute_expectation_0: Callable[[], float]) -> Outcome:
+    def measure(self, qubit: int, f_expectation0: Callable[[], float]) -> Outcome:
         """Return the measurement outcome of ``qubit``."""
-        expectation_0 = compute_expectation_0()
-        assert math.isclose(expectation_0, self.expected[qubit])
-        return super().measure(qubit, lambda: expectation_0)
+        expectation0 = f_expectation0()
+        assert math.isclose(expectation0, self.expected[qubit])
+        return super().measure(qubit, lambda: expectation0)
 
 
 @pytest.mark.filterwarnings("ignore:Simulating using densitymatrix backend with no noise.")

From b529a18f699cd5c996769238264f208c68e7e4f3 Mon Sep 17 00:00:00 2001
From: Thierry Martinez <thierry.martinez@inria.fr>
Date: Mon, 28 Jul 2025 21:49:55 +0200
Subject: [PATCH 14/19] Use a generic for FixedBranchSelector.results

---
 graphix/branch_selector.py    | 16 ++++++++--------
 tests/test_branch_selector.py |  3 ++-
 2 files changed, 10 insertions(+), 9 deletions(-)

diff --git a/graphix/branch_selector.py b/graphix/branch_selector.py
index 467636d07..90af7f404 100644
--- a/graphix/branch_selector.py
+++ b/graphix/branch_selector.py
@@ -10,21 +10,18 @@
 from __future__ import annotations
 
 from abc import ABC, abstractmethod
+from collections.abc import Mapping
 from dataclasses import dataclass
-from typing import TYPE_CHECKING, Callable
+from typing import TYPE_CHECKING, Callable, Generic, TypeVar
 
 from typing_extensions import override
 
-from graphix.measurements import outcome
+from graphix.measurements import Outcome, outcome
 from graphix.rng import ensure_rng
 
 if TYPE_CHECKING:
-    from collections.abc import Mapping
-
     from numpy.random import Generator
 
-    from graphix.measurements import Outcome
-
 
 class BranchSelector(ABC):
     """Abstract class for branch selectors.
@@ -86,8 +83,11 @@ def measure(self, qubit: int, f_expectation0: Callable[[], float]) -> Outcome:
         return result
 
 
+_T = TypeVar("_T", bound=Mapping[int, Outcome])
+
+
 @dataclass
-class FixedBranchSelector(BranchSelector):
+class FixedBranchSelector(BranchSelector, Generic[_T]):
     """Branch selector with predefined measurement outcomes.
 
     The mapping is fixed in ``results``. By default, an error is raised if
@@ -106,7 +106,7 @@ class FixedBranchSelector(BranchSelector):
         Default is ``None``.
     """
 
-    results: Mapping[int, Outcome]
+    results: _T
     default: BranchSelector | None = None
 
     @override
diff --git a/tests/test_branch_selector.py b/tests/test_branch_selector.py
index 07e7d36a2..83e35f2ba 100644
--- a/tests/test_branch_selector.py
+++ b/tests/test_branch_selector.py
@@ -135,7 +135,8 @@ def test_fixed_branch_selector(backend: str, outcome: list[Outcome]) -> None:
     ],
 )
 def test_fixed_branch_selector_no_default(backend: str) -> None:
-    branch_selector = FixedBranchSelector(results={})
+    results: dict[int, Outcome] = {}
+    branch_selector = FixedBranchSelector(results)
     pattern = Pattern(cmds=[N(0), M(0, angle=1e-5)])
     measure_method = DefaultMeasureMethod()
     with pytest.raises(ValueError):

From 092e63f851daa5a69479373c1be0b4fdb44d67ea Mon Sep 17 00:00:00 2001
From: Thierry Martinez <thierry.martinez@inria.fr>
Date: Mon, 28 Jul 2025 22:01:17 +0200
Subject: [PATCH 15/19] Import `Callable` from `collections.abc`

---
 graphix/branch_selector.py    | 4 +++-
 tests/test_branch_selector.py | 4 ++--
 2 files changed, 5 insertions(+), 3 deletions(-)

diff --git a/graphix/branch_selector.py b/graphix/branch_selector.py
index 90af7f404..d0a927276 100644
--- a/graphix/branch_selector.py
+++ b/graphix/branch_selector.py
@@ -12,7 +12,7 @@
 from abc import ABC, abstractmethod
 from collections.abc import Mapping
 from dataclasses import dataclass
-from typing import TYPE_CHECKING, Callable, Generic, TypeVar
+from typing import TYPE_CHECKING, Generic, TypeVar
 
 from typing_extensions import override
 
@@ -20,6 +20,8 @@
 from graphix.rng import ensure_rng
 
 if TYPE_CHECKING:
+    from collections.abc import Callable
+
     from numpy.random import Generator
 
 
diff --git a/tests/test_branch_selector.py b/tests/test_branch_selector.py
index 83e35f2ba..61e7f386c 100644
--- a/tests/test_branch_selector.py
+++ b/tests/test_branch_selector.py
@@ -4,7 +4,7 @@
 import itertools
 import math
 from dataclasses import dataclass
-from typing import TYPE_CHECKING, Callable
+from typing import TYPE_CHECKING
 
 import pytest
 from typing_extensions import override
@@ -15,7 +15,7 @@
 from graphix.simulator import DefaultMeasureMethod
 
 if TYPE_CHECKING:
-    from collections.abc import Mapping
+    from collections.abc import Callable, Mapping
 
     from numpy.random import Generator
 

From fd98b17895d061b87bd4e245e29bd8abeddb4fff Mon Sep 17 00:00:00 2001
From: Thierry Martinez <thierry.martinez@inria.fr>
Date: Mon, 28 Jul 2025 22:27:04 +0200
Subject: [PATCH 16/19] Pass `rng` as argument instead of keeping it in data
 structure

---
 graphix/branch_selector.py    | 26 ++++++++++++++------------
 graphix/pattern.py            | 14 ++++++++++++--
 graphix/sim/base_backend.py   | 15 ++++++++++++---
 graphix/sim/tensornet.py      | 12 +++++++-----
 graphix/simulator.py          | 34 +++++++++++++++++++---------------
 graphix/transpiler.py         | 25 ++++++++++++++-----------
 tests/test_branch_selector.py | 10 +++++-----
 tests/test_tnsim.py           |  4 ++--
 8 files changed, 85 insertions(+), 55 deletions(-)

diff --git a/graphix/branch_selector.py b/graphix/branch_selector.py
index d0a927276..8e7e54f9f 100644
--- a/graphix/branch_selector.py
+++ b/graphix/branch_selector.py
@@ -33,7 +33,7 @@ class BranchSelector(ABC):
     """
 
     @abstractmethod
-    def measure(self, qubit: int, f_expectation0: Callable[[], float]) -> Outcome:
+    def measure(self, qubit: int, f_expectation0: Callable[[], float], rng: Generator | None = None) -> Outcome:
         """Return the measurement outcome of ``qubit``.
 
         Parameters
@@ -46,6 +46,13 @@ def measure(self, qubit: int, f_expectation0: Callable[[], float]) -> Outcome:
             probability of outcome 0. The probability is computed only if
             this function is called (lazy computation), ensuring no
             unnecessary computational cost.
+
+        rng: Generator, optional
+            Random-number generator for measurements.
+            This generator is used only in case of random branch selection
+            (see :class:`RandomBranchSelector`).
+            If ``None``, a default random-number generator is used.
+            Default is ``None``.
         """
 
 
@@ -59,17 +66,12 @@ class RandomBranchSelector(BranchSelector):
         Whether to compute the probability distribution before selecting the measurement result.
         If ``False``, measurements yield 0/1 with equal probability (50% each).
         Default is ``True``.
-    rng : Generator | None, optional
-        Random-number generator for measurements.
-        If ``None``, a default random-number generator is used.
-        Default is ``None``.
     """
 
     pr_calc: bool = True
-    rng: Generator | None = None
 
     @override
-    def measure(self, qubit: int, f_expectation0: Callable[[], float]) -> Outcome:
+    def measure(self, qubit: int, f_expectation0: Callable[[], float], rng: Generator | None = None) -> Outcome:
         """
         Return the measurement outcome of ``qubit``.
 
@@ -77,11 +79,11 @@ def measure(self, qubit: int, f_expectation0: Callable[[], float]) -> Outcome:
         computed probability of outcome 0. Otherwise, the result is randomly chosen
         with a 50% chance for either outcome.
         """
-        self.rng = ensure_rng(self.rng)
+        rng = ensure_rng(rng)
         if self.pr_calc:
             prob_0 = f_expectation0()
-            return outcome(self.rng.random() > prob_0)
-        result: Outcome = self.rng.choice([0, 1])
+            return outcome(rng.random() > prob_0)
+        result: Outcome = rng.choice([0, 1])
         return result
 
 
@@ -112,7 +114,7 @@ class FixedBranchSelector(BranchSelector, Generic[_T]):
     default: BranchSelector | None = None
 
     @override
-    def measure(self, qubit: int, f_expectation0: Callable[[], float]) -> Outcome:
+    def measure(self, qubit: int, f_expectation0: Callable[[], float], rng: Generator | None = None) -> Outcome:
         """
         Return the predefined measurement outcome of ``qubit``, if available.
 
@@ -142,6 +144,6 @@ class ConstBranchSelector(BranchSelector):
     result: Outcome
 
     @override
-    def measure(self, qubit: int, f_expectation0: Callable[[], float]) -> Outcome:
+    def measure(self, qubit: int, f_expectation0: Callable[[], float], rng: Generator | None = None) -> Outcome:
         """Return the constant measurement outcome ``result`` for any qubit."""
         return self.result
diff --git a/graphix/pattern.py b/graphix/pattern.py
index 0cc7ec6ff..dfaae7e88 100644
--- a/graphix/pattern.py
+++ b/graphix/pattern.py
@@ -34,6 +34,8 @@
     from collections.abc import Set as AbstractSet
     from typing import Any
 
+    from numpy.random import Generator
+
     from graphix.parameter import ExpressionOrFloat, ExpressionOrSupportsFloat, Parameter
     from graphix.sim import Backend, BackendState, Data
 
@@ -1355,7 +1357,11 @@ def space_list(self) -> list[int]:
         return n_list
 
     def simulate_pattern(
-        self, backend: Backend[_StateT_co] | str = "statevector", input_state: Data = BasicStates.PLUS, **kwargs: Any
+        self,
+        backend: Backend[_StateT_co] | str = "statevector",
+        input_state: Data = BasicStates.PLUS,
+        rng: Generator | None = None,
+        **kwargs: Any,
     ) -> BackendState:
         """Simulate the execution of the pattern by using :class:`graphix.simulator.PatternSimulator`.
 
@@ -1365,6 +1371,10 @@ def simulate_pattern(
         ----------
         backend : str
             optional parameter to select simulator backend.
+        rng: Generator, optional
+            Random-number generator for measurements.
+            This generator is used only in case of random branch selection
+            (see :class:`RandomBranchSelector`).
         kwargs: keyword args for specified backend.
 
         Returns
@@ -1375,7 +1385,7 @@ def simulate_pattern(
         .. seealso:: :class:`graphix.simulator.PatternSimulator`
         """
         sim = PatternSimulator(self, backend=backend, **kwargs)
-        sim.run(input_state)
+        sim.run(input_state, rng=rng)
         return sim.backend.state
 
     def perform_pauli_measurements(self, leave_input: bool = False, ignore_pauli_with_deps: bool = False) -> None:
diff --git a/graphix/sim/base_backend.py b/graphix/sim/base_backend.py
index 3ce9ebd29..9ab263137 100644
--- a/graphix/sim/base_backend.py
+++ b/graphix/sim/base_backend.py
@@ -26,6 +26,8 @@
 if TYPE_CHECKING:
     from collections.abc import Collection, Iterable, Iterator, Sequence
 
+    from numpy.random import Generator
+
     from graphix import command
     from graphix.channels import KrausChannel
     from graphix.fundamentals import Plane
@@ -525,6 +527,7 @@ def perform_measure(
     angle: ExpressionOrFloat,
     state: DenseState,
     branch_selector: BranchSelector,
+    rng: Generator | None = None,
     symbolic: bool = False,
 ) -> Outcome:
     """Perform measurement of a qubit."""
@@ -544,7 +547,7 @@ def f_expectation0() -> float:
         assert math.isclose(exp_val.imag, 0, abs_tol=1e-10)
         return exp_val.real
 
-    result = branch_selector.measure(qubit_node, f_expectation0)
+    result = branch_selector.measure(qubit_node, f_expectation0, rng)
     op_mat = _op_mat_from_result(vec, 1, symbolic=symbolic) if result else get_op_mat0()
     state.evolve_single(op_mat, qubit_loc)
     return result
@@ -707,13 +710,17 @@ def finalize(self, output_nodes: Iterable[int]) -> None:
         """To be run at the end of pattern simulation to convey the order of output nodes."""
 
     @abstractmethod
-    def measure(self, node: int, measurement: Measurement) -> Outcome:
+    def measure(self, node: int, measurement: Measurement, rng: Generator | None = None) -> Outcome:
         """Perform measurement of a node and trace out the qubit.
 
         Parameters
         ----------
         node: int
         measurement: Measurement
+        rng: Generator, optional
+            Random-number generator for measurements.
+            This generator is used only in case of random branch selection
+            (see :class:`RandomBranchSelector`).
         """
 
 
@@ -791,13 +798,14 @@ def entangle_nodes(self, edge: tuple[int, int]) -> None:
         self.state.entangle((target, control))
 
     @override
-    def measure(self, node: int, measurement: Measurement) -> Outcome:
+    def measure(self, node: int, measurement: Measurement, rng: Generator | None = None) -> Outcome:
         """Perform measurement of a node and trace out the qubit.
 
         Parameters
         ----------
         node: int
         measurement: Measurement
+        rng: Generator, optional
         """
         loc = self.node_index.index(node)
         result = perform_measure(
@@ -807,6 +815,7 @@ def measure(self, node: int, measurement: Measurement) -> Outcome:
             measurement.angle,
             self.state,
             self.branch_selector,
+            rng=rng,
             symbolic=self.symbolic,
         )
         self.node_index.remove(node)
diff --git a/graphix/sim/tensornet.py b/graphix/sim/tensornet.py
index 77a89cd38..8892f62f6 100644
--- a/graphix/sim/tensornet.py
+++ b/graphix/sim/tensornet.py
@@ -30,6 +30,7 @@
     from collections.abc import Iterable, Sequence
 
     from cotengra.oe import PathOptimizer
+    from numpy.random import Generator
 
     from graphix import Pattern
     from graphix.clifford import Clifford
@@ -205,6 +206,7 @@ def measure_single(
         basis: str | npt.NDArray[np.complex128] = "Z",
         bypass_probability_calculation: bool = True,
         outcome: Outcome | None = None,
+        rng: Generator | None = None,
     ) -> Outcome:
         """Measure a node in specified basis. Note this does not perform the partial trace.
 
@@ -230,7 +232,7 @@ def measure_single(
             measurement result.
         """
         if bypass_probability_calculation:
-            result = outcome if outcome is not None else self.__branch_selector.measure(index, lambda: 0.5)
+            result = outcome if outcome is not None else self.__branch_selector.measure(index, lambda: 0.5, rng=rng)
             # Basis state to be projected
             if isinstance(basis, np.ndarray):
                 if outcome is not None:
@@ -739,7 +741,7 @@ def entangle_nodes(self, edge: tuple[int, int]) -> None:
             pass
 
     @override
-    def measure(self, node: int, measurement: Measurement) -> Outcome:
+    def measure(self, node: int, measurement: Measurement, rng: Generator | None = None) -> Outcome:
         """Perform measurement of the node.
 
         In the context of tensornetwork, performing measurement equals to
@@ -756,11 +758,11 @@ def measure(self, node: int, measurement: Measurement) -> Outcome:
             vector: npt.NDArray[np.complex128] = self.state.get_open_tensor_from_index(node)
             probs = (np.abs(vector) ** 2).astype(np.float64)
             probs /= np.sum(probs)
-            result: Outcome = self.branch_selector.measure(node, lambda: probs[0])
+            result: Outcome = self.branch_selector.measure(node, lambda: probs[0], rng=rng)
             self.results[node] = result
             buffer = 1 / probs[result] ** 0.5
         else:
-            result = self.branch_selector.measure(node, lambda: 0.5)
+            result = self.branch_selector.measure(node, lambda: 0.5, rng=rng)
             self.results[node] = result
             buffer = 2**0.5
         if isinstance(measurement.angle, Expression):
@@ -769,7 +771,7 @@ def measure(self, node: int, measurement: Measurement) -> Outcome:
         if result:
             vec = measurement.plane.orth.matrix @ vec
         proj_vec = vec * buffer
-        self.state.measure_single(node, basis=proj_vec)
+        self.state.measure_single(node, basis=proj_vec, rng=rng)
         return result
 
     @override
diff --git a/graphix/simulator.py b/graphix/simulator.py
index 4ed3153e3..222a71997 100644
--- a/graphix/simulator.py
+++ b/graphix/simulator.py
@@ -44,10 +44,16 @@ class MeasureMethod(abc.ABC):
     Example: class `ClientMeasureMethod` in https://github.com/qat-inria/veriphix
     """
 
-    def measure(self, backend: Backend[_StateT_co], cmd: BaseM, noise_model: NoiseModel | None = None) -> None:
+    def measure(
+        self,
+        backend: Backend[_StateT_co],
+        cmd: BaseM,
+        noise_model: NoiseModel | None = None,
+        rng: Generator | None = None,
+    ) -> None:
         """Perform a measure."""
         description = self.get_measurement_description(cmd)
-        result = backend.measure(cmd.node, description)
+        result = backend.measure(cmd.node, description, rng=rng)
         if noise_model is not None:
             result = noise_model.confuse_result(result)
         self.set_measure_result(cmd.node, result)
@@ -184,7 +190,6 @@ def __init__(
         measure_method: MeasureMethod | None = None,
         noise_model: NoiseModel | None = None,
         branch_selector: BranchSelector | None = None,
-        rng: Generator | None = None,
         graph_prep: str | None = None,
         symbolic: bool = False,
     ) -> None:
@@ -204,8 +209,6 @@ def __init__(
             [Density matrix backend only] Noise model used by the simulator.
         branch_selector: :class:`BranchSelector`, optional
             Branch selector used for measurements. Can only be specified if ``backend`` is not an already instantiated :class:`Backend` object.  Default is :class:`RandomBranchSelector`.
-        rng: :class:`numpy.random.Generator`, optional
-            Random number generator to be used by the default :class:`RandomBranchSelector`. Can only be specified if ``backend`` is not an already instantiated :class:`Backend` object and if ``branch_selector`` is not specified.
         graph_prep: str, optional
             [Tensor network backend only] Strategy for preparing the graph state.  See :class:`TensorNetworkBackend`.
         symbolic : bool, optional
@@ -217,23 +220,19 @@ def __init__(
         """
 
         def initialize_backend() -> Backend[BackendState]:
-            nonlocal backend, branch_selector, rng, graph_prep, noise_model
+            nonlocal backend, branch_selector, graph_prep, noise_model
             if isinstance(backend, Backend):
                 if noise_model is not None:
                     raise ValueError("`noise_model` cannot be specified if `backend` is already instantiated.")
                 if branch_selector is not None:
                     raise ValueError("`branch_selector` cannot be specified if `backend` is already instantiated.")
-                if rng is not None:
-                    raise ValueError("`rng` cannot be specified if `backend` is already instantiated.")
                 if graph_prep is not None:
                     raise ValueError("`graph_prep` cannot be specified if `backend` is already instantiated.")
                 if symbolic:
                     raise ValueError("`symbolic` cannot be specified if `backend` is already instantiated.")
                 return backend
             if branch_selector is None:
-                branch_selector = RandomBranchSelector(rng=rng)
-            elif rng is not None:
-                raise ValueError("`rng` and `branch_selector` cannot be specified simultaneously.")
+                branch_selector = RandomBranchSelector()
             if backend in {"tensornetwork", "mps"}:
                 if noise_model is not None:
                     raise ValueError("`noise_model` cannot be specified for tensor network backend.")
@@ -281,14 +280,19 @@ def set_noise_model(self, model: NoiseModel | None) -> None:
             raise ValueError(f"The backend {self.backend} doesn't support noise but noisemodel was provided.")
         self.noise_model = model
 
-    def run(self, input_state: Data = BasicStates.PLUS) -> None:
+    def run(self, input_state: Data = BasicStates.PLUS, rng: Generator | None = None) -> None:
         """Perform the simulation.
 
         Returns
         -------
-        state :
+        input_state: Data, optional
             the output quantum state,
             in the representation depending on the backend used.
+            Default: |+>.
+        rng: Generator, optional
+            Random-number generator for measurements.
+            This generator is used only in case of random branch selection
+            (see :class:`RandomBranchSelector`).
         """
         if input_state is not None:
             self.backend.add_nodes(self.pattern.input_nodes, input_state)
@@ -299,7 +303,7 @@ def run(self, input_state: Data = BasicStates.PLUS) -> None:
                 elif cmd.kind == CommandKind.E:
                     self.backend.entangle_nodes(edge=cmd.nodes)
                 elif cmd.kind == CommandKind.M:
-                    self.__measure_method.measure(self.backend, cmd)
+                    self.__measure_method.measure(self.backend, cmd, 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)
@@ -321,7 +325,7 @@ def run(self, input_state: Data = BasicStates.PLUS) -> None:
                     self.backend.apply_channel(self.noise_model.entangle(), cmd.nodes)
                 elif cmd.kind == CommandKind.M:
                     self.backend.apply_channel(self.noise_model.measure(), [cmd.node])
-                    self.__measure_method.measure(self.backend, cmd, noise_model=self.noise_model)
+                    self.__measure_method.measure(self.backend, cmd, noise_model=self.noise_model, rng=rng)
                 elif cmd.kind == CommandKind.X:
                     self.backend.correct_byproduct(cmd, self.__measure_method)
                     if np.mod(sum(self.__measure_method.get_measure_result(j) for j in cmd.domain), 2) == 1:
diff --git a/graphix/transpiler.py b/graphix/transpiler.py
index 34d46885a..0c288cb1d 100644
--- a/graphix/transpiler.py
+++ b/graphix/transpiler.py
@@ -887,14 +887,12 @@ def simulate_statevector(
         Parameters
         ----------
         input_state : Data
-        rng : Generator
-            Random number generator used to sample measurement outcomes.
-
         branch_selector: :class:`graphix.branch_selector.BranchSelector`
-            branch selector for measures (default: :class:`graphix.branch_selector.RandomBranchSelector`)
-
-        rng: :class:`np.random.Generator`
-            random number generator for :class:`graphix.branch_selector.RandomBranchSelector` (should only be used with default branch selector)
+            branch selector for measures (default: :class:`RandomBranchSelector`).
+        rng: Generator, optional
+            Random-number generator for measurements.
+            This generator is used only in case of random branch selection
+            (see :class:`RandomBranchSelector`).
 
         Returns
         -------
@@ -903,9 +901,7 @@ def simulate_statevector(
         """
         symbolic = self.is_parameterized()
         if branch_selector is None:
-            branch_selector = RandomBranchSelector(rng=rng)
-        elif rng is not None:
-            raise ValueError("Cannot specify both branch selector and rng")
+            branch_selector = RandomBranchSelector()
 
         state = Statevec(nqubit=self.width) if input_state is None else Statevec(nqubit=self.width, data=input_state)
 
@@ -941,7 +937,14 @@ def simulate_statevector(
                 state.evolve(Ops.CCX, [instr.controls[0], instr.controls[1], instr.target])
             elif instr.kind == instruction.InstructionKind.M:
                 result = base_backend.perform_measure(
-                    instr.target, instr.target, instr.plane, instr.angle * np.pi, state, branch_selector, symbolic
+                    instr.target,
+                    instr.target,
+                    instr.plane,
+                    instr.angle * np.pi,
+                    state,
+                    branch_selector,
+                    rng=rng,
+                    symbolic=symbolic,
                 )
                 classical_measures.append(result)
             else:
diff --git a/tests/test_branch_selector.py b/tests/test_branch_selector.py
index 61e7f386c..e801b9e69 100644
--- a/tests/test_branch_selector.py
+++ b/tests/test_branch_selector.py
@@ -31,7 +31,7 @@ class CheckedBranchSelector(RandomBranchSelector):
     expected: Mapping[int, float] = dataclasses.field(default_factory=dict)
 
     @override
-    def measure(self, qubit: int, f_expectation0: Callable[[], float]) -> Outcome:
+    def measure(self, qubit: int, f_expectation0: Callable[[], float], rng: Generator | None = None) -> Outcome:
         """Return the measurement outcome of ``qubit``."""
         expectation0 = f_expectation0()
         assert math.isclose(expectation0, self.expected[qubit])
@@ -55,8 +55,8 @@ def measure(self, qubit: int, f_expectation0: Callable[[], float]) -> Outcome:
 def test_expectation_value(fx_rng: Generator, backend: str) -> None:
     # Pattern that measures 0 on qubit 0 with probability 1.
     pattern = Pattern(cmds=[N(0), M(0)])
-    branch_selector = CheckedBranchSelector(rng=fx_rng, expected={0: 1.0})
-    pattern.simulate_pattern(backend, branch_selector=branch_selector)
+    branch_selector = CheckedBranchSelector(expected={0: 1.0})
+    pattern.simulate_pattern(backend, branch_selector=branch_selector, rng=fx_rng)
 
 
 @pytest.mark.filterwarnings("ignore:Simulating using densitymatrix backend with no noise.")
@@ -74,11 +74,11 @@ def test_expectation_value(fx_rng: Generator, backend: str) -> None:
     ],
 )
 def test_random_branch_selector(fx_rng: Generator, backend: str) -> None:
-    branch_selector = RandomBranchSelector(rng=fx_rng)
+    branch_selector = RandomBranchSelector()
     pattern = Pattern(cmds=[N(0), M(0)])
     for _ in range(NB_ROUNDS):
         measure_method = DefaultMeasureMethod()
-        pattern.simulate_pattern(backend, branch_selector=branch_selector, measure_method=measure_method)
+        pattern.simulate_pattern(backend, branch_selector=branch_selector, measure_method=measure_method, rng=fx_rng)
         assert measure_method.results[0] == 0
 
 
diff --git a/tests/test_tnsim.py b/tests/test_tnsim.py
index 14f6daa1e..cc669a590 100644
--- a/tests/test_tnsim.py
+++ b/tests/test_tnsim.py
@@ -34,7 +34,7 @@ def random_op(sites: int, dtype: type, rng: Generator) -> npt.NDArray:
 class TestTN:
     def test_add_node(self, fx_rng: Generator) -> None:
         node_index = fx_rng.integers(0, 1000)
-        tn = MBQCTensorNet(branch_selector=RandomBranchSelector(rng=fx_rng))
+        tn = MBQCTensorNet(branch_selector=RandomBranchSelector())
 
         tn.add_qubit(node_index)
 
@@ -43,7 +43,7 @@ def test_add_node(self, fx_rng: Generator) -> None:
 
     def test_add_nodes(self, fx_rng: Generator) -> None:
         node_index = set(fx_rng.integers(0, 1000, 20))
-        tn = MBQCTensorNet(branch_selector=RandomBranchSelector(rng=fx_rng))
+        tn = MBQCTensorNet(branch_selector=RandomBranchSelector())
 
         tn.graph_prep = "sequential"
         tn.add_qubits(node_index)

From b0be4ca5165efde82e7ba78d50bf38b4f39934b0 Mon Sep 17 00:00:00 2001
From: Thierry Martinez <thierry.martinez@inria.fr>
Date: Mon, 28 Jul 2025 22:40:24 +0200
Subject: [PATCH 17/19] Fix docstring syntax

---
 graphix/simulator.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/graphix/simulator.py b/graphix/simulator.py
index 222a71997..8dc363f7c 100644
--- a/graphix/simulator.py
+++ b/graphix/simulator.py
@@ -288,7 +288,7 @@ def run(self, input_state: Data = BasicStates.PLUS, rng: Generator | None = None
         input_state: Data, optional
             the output quantum state,
             in the representation depending on the backend used.
-            Default: |+>.
+            Default: ``|+>``.
         rng: Generator, optional
             Random-number generator for measurements.
             This generator is used only in case of random branch selection

From 58ac2177422b630586512bbd1aba3f19d4c36967 Mon Sep 17 00:00:00 2001
From: Thierry Martinez <thierry.martinez@inria.fr>
Date: Mon, 28 Jul 2025 22:44:15 +0200
Subject: [PATCH 18/19] Type annotations for pyright

---
 tests/test_branch_selector.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/tests/test_branch_selector.py b/tests/test_branch_selector.py
index e801b9e69..215c581c7 100644
--- a/tests/test_branch_selector.py
+++ b/tests/test_branch_selector.py
@@ -115,9 +115,9 @@ def test_random_branch_selector_without_pr_calc(backend: str) -> None:
 )
 @pytest.mark.parametrize("outcome", itertools.product([0, 1], repeat=3))
 def test_fixed_branch_selector(backend: str, outcome: list[Outcome]) -> None:
-    branch_selector = FixedBranchSelector(
-        results=dict(enumerate(outcome[:-1])), default=FixedBranchSelector({2: outcome[2]})
-    )
+    results1: dict[int, Outcome] = dict(enumerate(outcome[:-1]))
+    results2: dict[int, Outcome] = {2: outcome[2]}
+    branch_selector = FixedBranchSelector(results1, default=FixedBranchSelector(results2))
     pattern = Pattern(cmds=[cmd for qubit in range(3) for cmd in (N(qubit), M(qubit, angle=0.1))])
     measure_method = DefaultMeasureMethod()
     pattern.simulate_pattern(backend, branch_selector=branch_selector, measure_method=measure_method)

From 2dab4443e93b7bfb0fe5e86025e83bf58528b717 Mon Sep 17 00:00:00 2001
From: Thierry Martinez <thierry.martinez@inria.fr>
Date: Wed, 30 Jul 2025 17:50:49 +0200
Subject: [PATCH 19/19] Documentation update

---
 CHANGELOG.md              |  3 +++
 docs/source/simulator.rst | 29 +++++++++++++++++++++++++++++
 2 files changed, 32 insertions(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 542796872..09837a2e0 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -71,6 +71,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   The user can specify `pr_calc` in the constructor of
   `RandomBranchSelector` instead.
 
+- #300: `rng` is no longer stored in the backends; it is now passed as
+  an optional argument to each simulation method.
+
 - #261: Moved all device interface functionalities to an external
   library and removed their implementation from this library.
 
diff --git a/docs/source/simulator.rst b/docs/source/simulator.rst
index 78b103dc5..821a1be48 100644
--- a/docs/source/simulator.rst
+++ b/docs/source/simulator.rst
@@ -48,3 +48,32 @@ Density Matrix
 
 .. autoclass:: DensityMatrix
     :members:
+
+Branch Selection: :mod:`graphix.branch_selector` module
++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+.. currentmodule:: graphix.branch_selector
+
+Abstract Branch Selector
+------------------------
+
+.. autoclass:: BranchSelector
+    :members:
+
+Random Branch Selector
+----------------------
+
+.. autoclass:: RandomBranchSelector
+    :members:
+
+Fixed Branch Selector
+---------------------
+
+.. autoclass:: FixedBranchSelector
+    :members:
+
+Constant Branch Selector
+------------------------
+
+.. autoclass:: ConstBranchSelector
+    :members: