diff --git a/CHANGELOG.md b/CHANGELOG.md index 727372ff0..09837a2e0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -16,6 +16,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`. +- #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 + deriving the abstract class `BranchSelector`. + - #312: The separation between `TensorNetworkBackend` and backends that operate on a full-state representation, such as `StatevecBackend` and `DensityMatrixBackend`, is now clearer with @@ -61,6 +67,13 @@ 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 removed in back-end initializers. + 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 da35237a1..821a1be48 100644 --- a/docs/source/simulator.rst +++ b/docs/source/simulator.rst @@ -12,7 +12,6 @@ Pattern Simulation .. automethod:: run - Simulator backends ++++++++++++++++++ @@ -49,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: diff --git a/graphix/branch_selector.py b/graphix/branch_selector.py new file mode 100644 index 000000000..8e7e54f9f --- /dev/null +++ b/graphix/branch_selector.py @@ -0,0 +1,149 @@ +"""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 collections.abc import Mapping +from dataclasses import dataclass +from typing import TYPE_CHECKING, Generic, TypeVar + +from typing_extensions import override + +from graphix.measurements import Outcome, outcome +from graphix.rng import ensure_rng + +if TYPE_CHECKING: + from collections.abc import Callable + + 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, f_expectation0: Callable[[], float], rng: Generator | None = None) -> Outcome: + """Return the measurement outcome of ``qubit``. + + Parameters + ---------- + qubit : int + Index of qubit to measure + + 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 + 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``. + """ + + +@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``. + """ + + pr_calc: bool = True + + @override + def measure(self, qubit: int, f_expectation0: Callable[[], float], rng: Generator | None = None) -> Outcome: + """ + 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. + """ + rng = ensure_rng(rng) + if self.pr_calc: + prob_0 = f_expectation0() + return outcome(rng.random() > prob_0) + result: Outcome = rng.choice([0, 1]) + return result + + +_T = TypeVar("_T", bound=Mapping[int, Outcome]) + + +@dataclass +class FixedBranchSelector(BranchSelector, Generic[_T]): + """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: _T + default: BranchSelector | None = None + + @override + def measure(self, qubit: int, f_expectation0: Callable[[], float], rng: Generator | None = None) -> Outcome: + """ + 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, f_expectation0) + return result + + +@dataclass +class ConstBranchSelector(BranchSelector): + """Branch selector with a constant measurement outcome. + + The value ``result`` is returned for every qubit. + + Parameters + ---------- + result : Outcome + The fixed measurement outcome for all qubits. + """ + + result: Outcome + + @override + 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/measurements.py b/graphix/measurements.py index e9a7b9dff..fa382705e 100644 --- a/graphix/measurements.py +++ b/graphix/measurements.py @@ -22,6 +22,11 @@ def outcome(b: bool) -> Outcome: return 1 if b else 0 +def toggle_outcome(outcome: Outcome) -> Outcome: + """Toggle outcome.""" + return 1 if outcome == 0 else 0 + + @dataclasses.dataclass class Domains: """Represent `X^sZ^t` where s and t are XOR of results from given sets of indices.""" diff --git a/graphix/pattern.py b/graphix/pattern.py index 370c50f05..dfaae7e88 100644 --- a/graphix/pattern.py +++ b/graphix/pattern.py @@ -23,7 +23,7 @@ from graphix.fundamentals import Axis, Plane, Sign from graphix.gflow import find_flow, find_gflow, get_layers from graphix.graphsim import GraphState -from graphix.measurements import Outcome, PauliMeasurement +from graphix.measurements import Outcome, PauliMeasurement, toggle_outcome from graphix.pretty_print import OutputFormat, pattern_to_str from graphix.simulator import PatternSimulator from graphix.states import BasicStates @@ -32,7 +32,9 @@ if TYPE_CHECKING: from collections.abc import Container, Iterator, Mapping from collections.abc import Set as AbstractSet - from typing import Any, Literal + 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: @@ -1887,14 +1897,7 @@ def extract_signal(plane: Plane, s_domain: set[int], t_domain: set[int]) -> Extr assert_never(plane) -def toggle_outcome(outcome: Literal[0, 1]) -> Literal[0, 1]: - """Toggle outcome.""" - if outcome == 0: - return 1 - return 0 - - -def shift_outcomes(outcomes: dict[int, Literal[0, 1]], signal_dict: dict[int, set[int]]) -> dict[int, Literal[0, 1]]: +def shift_outcomes(outcomes: dict[int, Outcome], signal_dict: dict[int, set[int]]) -> dict[int, Outcome]: """Update outcomes with shifted signals. Shifted signals (as returned by the method diff --git a/graphix/sim/base_backend.py b/graphix/sim/base_backend.py index 85de43aa4..9ab263137 100644 --- a/graphix/sim/base_backend.py +++ b/graphix/sim/base_backend.py @@ -3,6 +3,7 @@ from __future__ import annotations import dataclasses +import math import sys from abc import ABC, abstractmethod from dataclasses import dataclass @@ -15,12 +16,11 @@ # override introduced in Python 3.12 from typing_extensions import TypeAlias, override +from graphix.branch_selector import BranchSelector, RandomBranchSelector from graphix.clifford import Clifford from graphix.command import CommandKind -from graphix.measurements import outcome from graphix.ops import Ops from graphix.parameter import check_expression_or_complex -from graphix.rng import ensure_rng from graphix.states import BasicStates if TYPE_CHECKING: @@ -521,27 +521,35 @@ def _op_mat_from_result( def perform_measure( - qubit: int, + qubit_node: int, + qubit_loc: int, plane: Plane, angle: ExpressionOrFloat, state: DenseState, - rng: Generator, - pr_calc: bool = True, + branch_selector: BranchSelector, + rng: Generator | None = None, symbolic: bool = False, ) -> Outcome: """Perform measurement of a qubit.""" vec = plane.polar(angle) - if pr_calc: - op_mat = _op_mat_from_result(vec, 0, symbolic=symbolic) - prob_0 = state.expectation_single(op_mat, qubit) - result = outcome(rng.random() > abs(prob_0)) - if result: - op_mat = _op_mat_from_result(vec, 1, 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_mat0 may contain the matrix operator associated with the outcome 0, + # but the value is computed lazily, i.e., only if needed. + op_mat0 = None + + 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 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, 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 @@ -702,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`). """ @@ -737,11 +749,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. @@ -751,8 +760,7 @@ class DenseStateBackend(Backend[_DenseStateT_co], Generic[_DenseStateT_co]): """ node_index: NodeIndex = dataclasses.field(default_factory=NodeIndex) - pr_calc: bool = True - rng: Generator = dataclasses.field(default_factory=ensure_rng) + branch_selector: BranchSelector = dataclasses.field(default_factory=RandomBranchSelector) symbolic: bool = False @override @@ -790,22 +798,24 @@ 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( + node, loc, measurement.plane, measurement.angle, self.state, - rng=self.rng, - pr_calc=self.pr_calc, + 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 fb6284dfe..8892f62f6 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 @@ -53,7 +53,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 +82,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. @@ -206,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. @@ -231,7 +232,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, rng=rng) # Basis state to be projected if isinstance(basis, np.ndarray): if outcome is not None: @@ -537,7 +538,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 +582,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 +610,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 +627,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 +652,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 +664,7 @@ def __init__( pattern, graph_prep, input_state, - rng, + branch_selector, pattern.output_nodes, results, decomposed_cz, @@ -735,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 @@ -752,12 +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.rng.choice([0, 1], p=probs) + result: Outcome = self.branch_selector.measure(node, lambda: probs[0], rng=rng) 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, rng=rng) self.results[node] = result buffer = 2**0.5 if isinstance(measurement.angle, Expression): @@ -766,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 48e982f95..8dc363f7c 100644 --- a/graphix/simulator.py +++ b/graphix/simulator.py @@ -13,6 +13,7 @@ import numpy as np from graphix import command +from graphix.branch_selector import BranchSelector, RandomBranchSelector from graphix.clifford import Clifford from graphix.command import BaseM, CommandKind, MeasureUpdate from graphix.measurements import Measurement, Outcome @@ -24,7 +25,8 @@ if TYPE_CHECKING: from collections.abc import Mapping - from typing import Any + + from numpy.random import Generator from graphix.noise_models.noise_model import NoiseModel from graphix.pattern import Pattern @@ -42,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) @@ -181,47 +189,75 @@ def __init__( backend: Backend[BackendState] | str = "statevector", measure_method: MeasureMethod | None = None, noise_model: NoiseModel | None = None, - **kwargs: Any, + branch_selector: BranchSelector | None = None, + graph_prep: str | None = None, + symbolic: bool = False, ) -> None: """ Construct a pattern simulator. 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'. - noise_model: - kwargs: keyword args for specified backend. + measure_method: :class:`MeasureMethod`, optional + Measure method used by the simulator. Default is :class:`DefaultMeasureMethod`. + noise_model: :class:`NoiseModel`, optional + [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`. + graph_prep: str, optional + [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. .. seealso:: :class:`graphix.sim.statevec.StatevectorBackend`\ :class:`graphix.sim.tensornet.TensorNetworkBackend`\ :class:`graphix.sim.density_matrix.DensityMatrixBackend`\ """ - if isinstance(backend, Backend): - assert kwargs == {} - self.backend = backend - elif backend == "statevector": - self.backend = StatevectorBackend(**kwargs) - elif backend == "densitymatrix": - if noise_model is None: - self.noise_model = None - self.backend = DensityMatrixBackend(**kwargs) - 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, - ) - else: - self.backend = DensityMatrixBackend(pr_calc=True, **kwargs) - self.set_noise_model(noise_model) - elif backend in {"tensornetwork", "mps"} and noise_model is None: - self.noise_model = None - self.backend = TensorNetworkBackend(pattern, **kwargs) - else: - raise ValueError("Unknown backend.") - self.set_noise_model(noise_model) + + def initialize_backend() -> Backend[BackendState]: + 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 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() + 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) + 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.") + 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, symbolic=symbolic) + raise ValueError(f"Unknown backend {backend}.") + + self.backend = initialize_backend() + self.noise_model = noise_model self.__pattern = pattern if measure_method is None: measure_method = DefaultMeasureMethod(pattern.results) @@ -244,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) @@ -262,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) @@ -284,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 2b3c675a6..0c288cb1d 100644 --- a/graphix/transpiler.py +++ b/graphix/transpiler.py @@ -13,13 +13,13 @@ from typing_extensions import assert_never from graphix import command, instruction, parameter +from graphix.branch_selector import BranchSelector, RandomBranchSelector from graphix.command import E, M, N, X, Z from graphix.fundamentals import Plane from graphix.instruction import Instruction, InstructionKind from graphix.ops import Ops from graphix.parameter import ExpressionOrFloat, Parameter from graphix.pattern import Pattern -from graphix.rng import ensure_rng from graphix.sim import Data, Statevec, base_backend if TYPE_CHECKING: @@ -876,14 +876,23 @@ def _ccx_command( ) return ancilla[17], ancilla[15], ancilla[13], seq - def simulate_statevector(self, input_state: Data | None = None, rng: Generator | 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 : Data - rng : Generator - Random number generator used to sample measurement outcomes. + branch_selector: :class:`graphix.branch_selector.BranchSelector` + 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 ------- @@ -891,7 +900,9 @@ def simulate_statevector(self, input_state: Data | None = None, rng: Generator | output state of the statevector simulation and results of classical measures. """ symbolic = self.is_parameterized() - rng = ensure_rng(rng) + if branch_selector is None: + branch_selector = RandomBranchSelector() + state = Statevec(nqubit=self.width) if input_state is None else Statevec(nqubit=self.width, data=input_state) classical_measures = [] @@ -926,7 +937,14 @@ def simulate_statevector(self, input_state: Data | None = None, rng: Generator | 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.plane, instr.angle * np.pi, state, rng, 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/noxfile.py b/noxfile.py index b64c9a64a..4f19886e6 100644 --- a/noxfile.py +++ b/noxfile.py @@ -62,6 +62,8 @@ 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") + # 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"): session.run("pytest", "--doctest-modules") diff --git a/tests/test_branch_selector.py b/tests/test_branch_selector.py new file mode 100644 index 000000000..215c581c7 --- /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 + +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 Callable, 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.""" + + expected: Mapping[int, float] = dataclasses.field(default_factory=dict) + + @override + 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]) + return super().measure(qubit, lambda: expectation0) + + +@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(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.") +@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() + 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, rng=fx_rng) + 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 / 5 + + +@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: + 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) + 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: + 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): + 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 diff --git a/tests/test_density_matrix.py b/tests/test_density_matrix.py index 35b376b7a..ca0210c00 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 @@ -921,7 +922,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 4f452894a..f3c82e750 100644 --- a/tests/test_pattern.py +++ b/tests/test_pattern.py @@ -10,10 +10,11 @@ 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, Command, CommandKind, E, M, N, X, Z from graphix.fundamentals import Plane -from graphix.measurements import PauliMeasurement +from graphix.measurements import Outcome, PauliMeasurement from graphix.pattern import Pattern, shift_outcomes from graphix.random_objects import rand_circuit, rand_gate from graphix.sim.density_matrix import DensityMatrix @@ -24,9 +25,7 @@ from graphix.transpiler import Circuit if TYPE_CHECKING: - import collections.abc from collections.abc import Sequence - from typing import Literal from graphix.sim.base_backend import BackendState @@ -39,17 +38,6 @@ def compare_backend_result_with_statevec(backend_state: BackendState, statevec: raise NotImplementedError(backend_state) -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() @@ -225,7 +213,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(0)) assert np.abs(np.dot(state.flatten().conjugate(), state_ref.flatten())) == pytest.approx(1) @pytest.mark.parametrize("jumps", range(1, 11)) @@ -344,14 +332,14 @@ def test_shift_signals_plane(self, plane: Plane, method: str) -> None: pattern.standardize() signal_dict = pattern.shift_signals(method=method) # Test for every possible outcome of each measure - zero_one: list[Literal[0, 1]] = [0, 1] - outcomes_ref: tuple[Literal[0, 1], ...] - for outcomes_ref in itertools.product(*([zero_one] * 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))) - ) + zero_one: list[Outcome] = [0, 1] + for outcomes_ref_list in itertools.product(*([zero_one] * 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)) diff --git a/tests/test_tnsim.py b/tests/test_tnsim.py index 0410f7f71..cc669a590 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()) 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()) tn.graph_prep = "sequential" tn.add_qubits(node_index)