diff --git a/docs/source/tutorial.rst b/docs/source/tutorial.rst index a828719e7..c5e75a9b0 100644 --- a/docs/source/tutorial.rst +++ b/docs/source/tutorial.rst @@ -304,7 +304,7 @@ With the noise model written, we can simulate it. dm_result = simulator.run() ->>> print(dm_result.fidelity(out_state.psi.flatten())) +>>> print(dm_result.fidelity(out_state.flatten())) 0.9718678141724848 We can plot the results from the model, @@ -318,7 +318,7 @@ We can plot the results from the model, for i in range(10): simulator = PatternSimulator(pattern, backend="densitymatrix", noise_model=NoisyGraphState(p_z=err_arr[i])) dm_result = simulator.run() - fidelity[i] = dm_result.fidelity(out_state.psi.flatten()) + fidelity[i] = dm_result.fidelity(out_state.flatten()) plt.semilogx(err_arr, fidelity, "o:") plt.xlabel("dephasing error of qubit preparation") diff --git a/examples/deutsch_jozsa.py b/examples/deutsch_jozsa.py index 68a032846..aaf7c085b 100644 --- a/examples/deutsch_jozsa.py +++ b/examples/deutsch_jozsa.py @@ -91,6 +91,6 @@ out_state = pattern.simulate_pattern() state = circuit.simulate_statevector().statevec -print("overlap of states: ", np.abs(np.dot(state.psi.flatten().conjugate(), out_state.psi.flatten()))) +print("overlap of states: ", np.abs(np.dot(state.flatten().conjugate(), out_state.flatten()))) # %% diff --git a/examples/qaoa.py b/examples/qaoa.py index f0f2c1a42..a4a54ac67 100644 --- a/examples/qaoa.py +++ b/examples/qaoa.py @@ -48,7 +48,7 @@ out_state = pattern.simulate_pattern() state = circuit.simulate_statevector().statevec -print("overlap of states: ", np.abs(np.dot(state.psi.flatten().conjugate(), out_state.psi.flatten()))) +print("overlap of states: ", np.abs(np.dot(state.flatten().conjugate(), out_state.flatten()))) # sphinx_gallery_thumbnail_number = 2 # %% diff --git a/examples/rotation.py b/examples/rotation.py index 0de287855..44a1c1621 100644 --- a/examples/rotation.py +++ b/examples/rotation.py @@ -77,7 +77,7 @@ state = Statevec(nqubit=2, data=BasicStates.ZERO) # starts with |0> states state.evolve_single(Ops.rx(theta[0]), 0) state.evolve_single(Ops.rx(theta[1]), 1) -print("overlap of states: ", np.abs(np.dot(state.psi.flatten().conjugate(), out_state.psi.flatten()))) +print("overlap of states: ", np.abs(np.dot(state.flatten().conjugate(), out_state.flatten()))) # %% # Now let us compile more complex pattern and inspect the graph using the visualization tool. diff --git a/graphix/sim/base_backend.py b/graphix/sim/base_backend.py index ef1977b3c..64fdcbb13 100644 --- a/graphix/sim/base_backend.py +++ b/graphix/sim/base_backend.py @@ -6,7 +6,7 @@ import math from abc import ABC, abstractmethod from dataclasses import dataclass -from typing import TYPE_CHECKING, Generic, SupportsFloat, TypeAlias, TypeVar +from typing import TYPE_CHECKING, Generic, SupportsFloat, TypeAlias, TypedDict, TypeVar import numpy as np import numpy.typing as npt @@ -364,15 +364,14 @@ def flatten(self) -> Matrix: @abstractmethod def add_nodes(self, nqubit: int, data: Data) -> None: - """ - Add nodes (qubits) to the state and initialize them in a specified state. + """Add nodes (qubits) to the state and initialize them in a specified state. Parameters ---------- nqubit : int The number of qubits to add to the state. - data : Data, optional + data : Data The state in which to initialize the newly added nodes. The supported forms of state specification depend on the backend implementation. @@ -380,67 +379,94 @@ def add_nodes(self, nqubit: int, data: Data) -> None: """ @abstractmethod - def entangle(self, edge: tuple[int, int]) -> None: - """Connect graph nodes. + def entangle(self, qubits: tuple[int, int]) -> None: + """Apply a CZ gate on two qubits. Parameters ---------- - edge : tuple of int - (control, target) qubit indices + qubits : tuple[int, int] + (control, target) qubit indices. """ @abstractmethod - def evolve(self, op: Matrix, qargs: Sequence[int]) -> None: - """Apply a multi-qubit operation. + def evolve(self, op: Matrix, qubits: Sequence[int]) -> None: + """Apply a multi-qubit operator. Parameters ---------- - op : numpy.ndarray - 2^n*2^n matrix - qargs : list of int - target qubits' indices + op : Matrix + Matrix of shape :math:`(2^n, 2^n)` representing + the operator to apply. + qubits : Sequence[int] + Target qubit indices. """ @abstractmethod - def evolve_single(self, op: Matrix, i: int) -> None: + def evolve_single(self, op: Matrix, qubit: int) -> None: """Apply a single-qubit operation. Parameters ---------- - op : numpy.ndarray - 2*2 matrix - i : int - qubit index + op : Matrix + Matrix of shape :math:`(2, 2)` representing + the operator to apply. + qubit : int + Target qubit index. """ @abstractmethod - def expectation_single(self, op: Matrix, loc: int) -> complex: - """Return the expectation value of single-qubit operator. + def expectation_single(self, op: Matrix, qubit: int) -> complex: + """Return the expectation value of a single-qubit operator. Parameters ---------- - op : numpy.ndarray - 2*2 operator - loc : int - target qubit index + op : Matrix + Matrix of shape :math:`(2, 2)` representing + the operator to measure. + qubit : int + Target qubit index. Returns ------- - complex : expectation value. + complex + Expectation value. + """ @abstractmethod - def remove_qubit(self, qarg: int) -> None: - """Remove a separable qubit from the system.""" + def remove_qubit(self, qubit: int) -> None: + """Remove a separable qubit from the system. + + Parameters + ---------- + qubit : int + Target qubit index. + """ + + def project_qubit(self, op: Matrix, qubit: int) -> None: + r"""Project out a qubit from the system and assemble the statevector of the remaining qubits. + + This method combines :meth:`evolve_single` and :meth:`remove_qubit`. It evolves the statevector with ``op`` and removes ``qubit``. It assumes that after the application of ``op``, ``qubit`` is a separable qubit. + + Parameters + ---------- + op : npt.NDArray[np.complex128] + Complex-valued matrix of shape :math:`(2, 2)` representing + the projector to apply. + qubit : int + Target qubit index. + """ + self.evolve_single(op, qubit) + self.remove_qubit(qubit) @abstractmethod def swap(self, qubits: tuple[int, int]) -> None: - """Swap qubits. + """Apply SWAP gate between two qubits. Parameters ---------- - qubits : tuple of int - (control, target) qubit indices + qubits : tuple[int, int] + (control, target) qubit indices. """ def apply_noise(self, qubits: Sequence[int], noise: Noise) -> None: # noqa: ARG002 @@ -682,6 +708,14 @@ def measure( _DenseStateT_co = TypeVar("_DenseStateT_co", bound="DenseState", covariant=True) +class DenseStateBackendKwargs(TypedDict, total=False): + """Keyword arguments for initializing a `DenseStateBackend`.""" + + node_index: NodeIndex + branch_selector: BranchSelector + symbolic: bool + + @dataclass(frozen=True) class DenseStateBackend(Backend[_DenseStateT_co], Generic[_DenseStateT_co]): """ @@ -709,11 +743,14 @@ class DenseStateBackend(Backend[_DenseStateT_co], Generic[_DenseStateT_co]): symbolic : bool, optional If True, support arbitrary objects (typically, symbolic expressions) in matrices. + All parameters are key-word only. + See Also -------- :class:`StatevecBackend`, :class:`DensityMatrixBackend`, :class:`TensorNetworkBackend` """ + _: dataclasses.KW_ONLY node_index: NodeIndex = dataclasses.field(default_factory=NodeIndex) branch_selector: BranchSelector = dataclasses.field(default_factory=RandomBranchSelector) symbolic: bool = False @@ -756,13 +793,23 @@ def entangle_nodes(self, edge: tuple[int, int]) -> None: def measure( self, node: int, measurement: Measurement, rng: Generator | None = None, *, stacklevel: int = 1 ) -> Outcome: - """Perform measurement of a node and trace out the qubit. + """Measure a node and trace out the corresponding qubit. Parameters ---------- - node: int - measurement: Measurement - rng: Generator, optional + node : int + Index of the node to measure. + measurement : Measurement + Measurement specification defining the measurement plane and angle. + rng : Generator, optional + Random number generator used for probabilistic outcome sampling. + stacklevel : int, default=1 + Stack level passed to the branch selector for warning reporting. + + Returns + ------- + Outcome + Measurement outcome. """ loc = self.node_index.index(node) bloch = measurement.to_bloch() @@ -784,9 +831,8 @@ def f_expectation0() -> float: outcome = self.branch_selector.measure(node, f_expectation0, rng, stacklevel=stacklevel + 1) op_mat = _outcome_to_operator_matrix(vec, 1, symbolic=self.symbolic) if outcome else compute_op_mat0() - self.state.evolve_single(op_mat, loc) + self.state.project_qubit(op_mat, loc) self.node_index.remove(node) - self.state.remove_qubit(loc) return outcome @override @@ -810,7 +856,7 @@ def apply_noise(self, cmd: ApplyNoise) -> None: def apply_single(self, node: int, op: Matrix) -> None: """Apply a single gate to the state.""" index = self.node_index.index(node) - self.state.evolve_single(op=op, i=index) + self.state.evolve_single(op=op, qubit=index) @override def apply_clifford(self, node: int, clifford: Clifford) -> None: diff --git a/graphix/sim/density_matrix.py b/graphix/sim/density_matrix.py index cb5570e2a..1f04d0859 100644 --- a/graphix/sim/density_matrix.py +++ b/graphix/sim/density_matrix.py @@ -20,7 +20,7 @@ from graphix.channels import KrausChannel from graphix.parameter import Expression, ExpressionOrFloat, ExpressionOrSupportsComplex from graphix.sim.base_backend import DenseState, DenseStateBackend, Matrix, kron, matmul, outer, tensordot, vdot -from graphix.sim.statevec import CNOT_TENSOR, CZ_TENSOR, SWAP_TENSOR, Statevec +from graphix.sim.statevec import Statevec from graphix.states import BasicStates, State if TYPE_CHECKING: @@ -31,6 +31,19 @@ from graphix.parameter import ExpressionOrSupportsFloat, Parameter from graphix.sim.data import Data +CZ_TENSOR = np.array( + [[[[1, 0], [0, 0]], [[0, 1], [0, 0]]], [[[0, 0], [1, 0]], [[0, 0], [0, -1]]]], + dtype=np.complex128, +) +CNOT_TENSOR = np.array( + [[[[1, 0], [0, 0]], [[0, 1], [0, 0]]], [[[0, 0], [0, 1]], [[0, 0], [1, 0]]]], + dtype=np.complex128, +) +SWAP_TENSOR = np.array( + [[[[1, 0], [0, 0]], [[0, 0], [1, 0]]], [[[0, 1], [0, 0]], [[0, 0], [0, 1]]]], + dtype=np.complex128, +) + class DensityMatrix(DenseState): """DensityMatrix object.""" @@ -146,33 +159,33 @@ def add_nodes(self, nqubit: int, data: Data) -> None: self.tensor(dm_to_add) @override - def evolve_single(self, op: Matrix, i: int) -> None: + def evolve_single(self, op: Matrix, qubit: int) -> None: """Single-qubit operation. Parameters ---------- op : np.ndarray 2*2 matrix. - i : int + qubit : int Index of qubit to apply operator. """ - assert i >= 0 - assert i < self.nqubit + assert qubit >= 0 + assert qubit < self.nqubit if op.shape != (2, 2): raise ValueError("op must be 2*2 matrix.") rho_tensor = self.rho.reshape((2,) * self.nqubit * 2) - rho_tensor = tensordot(tensordot(op, rho_tensor, axes=(1, i)), op.conj().T, axes=(i + self.nqubit, 0)) - rho_tensor = np.moveaxis(rho_tensor, (0, -1), (i, i + self.nqubit)) + rho_tensor = tensordot(tensordot(op, rho_tensor, axes=(1, qubit)), op.conj().T, axes=(qubit + self.nqubit, 0)) + rho_tensor = np.moveaxis(rho_tensor, (0, -1), (qubit, qubit + self.nqubit)) self.rho = rho_tensor.reshape((2**self.nqubit, 2**self.nqubit)) @override - def evolve(self, op: Matrix, qargs: Sequence[int]) -> None: + def evolve(self, op: Matrix, qubits: Sequence[int]) -> None: """Multi-qubit operation. Args: op (np.array): 2^n*2^n matrix - qargs (list of ints): target qubits' indexes + qubits (list of ints): target qubits' indexes """ d = op.shape # check it is a matrix. @@ -190,12 +203,12 @@ def evolve(self, op: Matrix, qargs: Sequence[int]) -> None: raise ValueError("Incorrect operator dimension: not consistent with qubits.") nqb_op = int(nqb_op) - if nqb_op != len(qargs): + if nqb_op != len(qubits): raise ValueError("The dimension of the operator doesn't match the number of targets.") - if not all(0 <= i < self.nqubit for i in qargs): + if not all(0 <= i < self.nqubit for i in qubits): raise ValueError("Incorrect target indices.") - if len(set(qargs)) != nqb_op: + if len(set(qubits)) != nqb_op: raise ValueError("A repeated target qubit index is not possible.") op_tensor = op.reshape((2,) * 2 * nqb_op) @@ -203,31 +216,31 @@ def evolve(self, op: Matrix, qargs: Sequence[int]) -> None: rho_tensor = self.rho.reshape((2,) * self.nqubit * 2) rho_tensor = tensordot( - tensordot(op_tensor, rho_tensor, axes=(tuple(nqb_op + i for i in range(len(qargs))), tuple(qargs))), + tensordot(op_tensor, rho_tensor, axes=(tuple(nqb_op + i for i in range(len(qubits))), tuple(qubits))), op.conj().T.reshape((2,) * 2 * nqb_op), - axes=(tuple(i + self.nqubit for i in qargs), tuple(i for i in range(len(qargs)))), + axes=(tuple(i + self.nqubit for i in qubits), tuple(i for i in range(len(qubits)))), ) rho_tensor = np.moveaxis( rho_tensor, - list(range(len(qargs))) + [-i for i in range(1, len(qargs) + 1)], - list(qargs) + [i + self.nqubit for i in reversed(list(qargs))], + list(range(len(qubits))) + [-i for i in range(1, len(qubits) + 1)], + list(qubits) + [i + self.nqubit for i in reversed(list(qubits))], ) self.rho = rho_tensor.reshape((2**self.nqubit, 2**self.nqubit)) @override - def expectation_single(self, op: Matrix, loc: int) -> complex: + def expectation_single(self, op: Matrix, qubit: int) -> complex: """Return the expectation value of single-qubit operator. Args: op (np.array): 2*2 Hermite operator - loc (int): Index of qubit on which to apply operator. + qubit (int): Index of qubit on which to apply operator. Returns ------- complex: expectation value (real for hermitian ops!). """ - if not (0 <= loc < self.nqubit): - raise ValueError(f"Wrong target qubit {loc}. Must between 0 and {self.nqubit - 1}.") + if not (0 <= qubit < self.nqubit): + raise ValueError(f"Wrong target qubit {qubit}. Must between 0 and {self.nqubit - 1}.") if op.shape != (2, 2): raise ValueError("op must be 2x2 matrix.") @@ -237,8 +250,8 @@ def expectation_single(self, op: Matrix, loc: int) -> complex: nqubit = self.nqubit rho_tensor: Matrix = st1.rho.reshape((2,) * nqubit * 2) - rho_tensor = tensordot(op, rho_tensor, axes=((1,), (loc,))) - rho_tensor = np.moveaxis(rho_tensor, 0, loc) + rho_tensor = tensordot(op, rho_tensor, axes=((1,), (qubit,))) + rho_tensor = np.moveaxis(rho_tensor, 0, qubit) # complex() needed with mypy strict mode (no-any-return) return complex(np.trace(rho_tensor.reshape((2**nqubit, 2**nqubit)))) @@ -282,15 +295,15 @@ def swap(self, qubits: tuple[int, int]) -> None: """ self.evolve(SWAP_TENSOR.reshape(4, 4), qubits) - def entangle(self, edge: tuple[int, int]) -> None: + def entangle(self, qubits: tuple[int, int]) -> None: """Connect graph nodes. Parameters ---------- - edge : (int, int) or [int, int] + qubits : (int, int) or [int, int] (control, target) qubit indices. """ - self.evolve(CZ_TENSOR.reshape(4, 4), edge) + self.evolve(CZ_TENSOR.reshape(4, 4), qubits) def normalize(self) -> None: """Normalize density matrix.""" @@ -306,9 +319,9 @@ def normalize(self) -> None: rho_c /= np.trace(rho_c) @override - def remove_qubit(self, qarg: int) -> None: + def remove_qubit(self, qubit: int) -> None: """Remove a qubit.""" - self.ptrace(qarg) + self.ptrace(qubit) self.normalize() def ptrace(self, qargs: Collection[int] | int) -> None: diff --git a/graphix/sim/statevec.py b/graphix/sim/statevec.py index 1f7dd23e6..63238ca95 100644 --- a/graphix/sim/statevec.py +++ b/graphix/sim/statevec.py @@ -1,90 +1,143 @@ -"""MBQC state vector backend.""" +"""Numba JIT-compiled MBQC statevector backend simulator based on Ref. [1]. + +[1] McGuffin, M. J., Robert J-M., and Ikeda K. "How to Write a Simulator for Quantum Circuits from Scratch: A Tutorial.", 2025 (arXiv:2506.08142). +""" from __future__ import annotations -import copy import dataclasses import functools import math from collections.abc import Iterable +from copy import deepcopy from dataclasses import dataclass from typing import TYPE_CHECKING, SupportsComplex, SupportsFloat +import numba as nb import numpy as np import numpy.typing as npt from typing_extensions import override -from graphix import parameter, states -from graphix.parameter import Expression, ExpressionOrSupportsComplex, check_expression_or_float -from graphix.sim.base_backend import DenseState, DenseStateBackend, Matrix, kron, tensordot +from graphix import states +from graphix.sim.base_backend import ( + DenseState, + DenseStateBackend, + DenseStateBackendKwargs, + Matrix, +) from graphix.states import BasicStates if TYPE_CHECKING: - from collections.abc import Callable, Mapping, Sequence - from typing import Any, Literal, TypeVar + from collections.abc import Callable, Sequence + from typing import Any, Literal, Self, TypeVar + + # Unpack introduced in Python 3.12 + from typing_extensions import Unpack - from graphix.parameter import ExpressionOrFloat, ExpressionOrSupportsFloat, Parameter + from graphix.parameter import ExpressionOrSupportsComplex from graphix.sim.data import Data _ENCODING = Literal["LSB", "MSB"] _ScalarT = TypeVar("_ScalarT", bound=np.generic[Any]) -CZ_TENSOR = np.array( - [[[[1, 0], [0, 0]], [[0, 1], [0, 0]]], [[[0, 0], [1, 0]], [[0, 0], [0, -1]]]], - dtype=np.complex128, -) -CNOT_TENSOR = np.array( - [[[[1, 0], [0, 0]], [[0, 1], [0, 0]]], [[[0, 0], [0, 1]], [[0, 0], [1, 0]]]], - dtype=np.complex128, -) -SWAP_TENSOR = np.array( - [[[[1, 0], [0, 0]], [[0, 0], [1, 0]]], [[[0, 1], [0, 0]], [[0, 0], [0, 1]]]], - dtype=np.complex128, -) +NUM_QUBIT_PARALLEL = 15 +"""This compilation constant determines the number of qubits above which matrix operations +are multi-threaded. For lower counts, the overhead does not compensate parallelization. +This number was determined empirically and it may be platform dependent.""" class Statevec(DenseState): - """Statevector object.""" - - psi: Matrix - - def __init__( - self, - data: Data = BasicStates.PLUS, - nqubit: int | None = None, - ) -> None: - """Initialize statevector objects. + """Statevector object. + + Attributes + ---------- + _psi : npt.NDArray[np.complex128] + Complex-valued 1-dimensional array representing the quantum statevector. + Only the first ``2**nqubit`` complex values have meaning. + + _nqubit : int + Number of active qubits at any given time. + + _max_qubits : int + Maximum Hilbert space size allowed for internal computations. It determines + the size of ``self._psi``. For circuit simulations, it corresponds to the number + of qubits, while for pattern simulations it corresponds to the pattern's + maximum space. The method :meth:`Statevec.ensure_capacity` allows to increase + this number. + + Notes + ----- + The internal representation of the quantum state is guaranteed to be + normalized after initialization, and it is assumed to remain normalized + thereafter. + + Using :meth:`evolve_single`, :meth:`expectation_single`, :meth:`evolve`, + or :meth:`expectation_value` with non-unitary operators does not preserve + the norm of the statevector and may lead to unexpected behavior. + + In pattern simulation, node measurements call :meth:`evolve_single` with + a projector (a non-unitary operator). However, the measured qubit is + immediately removed via :meth:`remove_qubit`, which restores the unit + norm of the internal quantum state. + See :meth:`graphix.sim.base_backend.DenseStateBackend.measure` for additional + details. + """ + + _psi: npt.NDArray[np.complex128] + _nqubit: int + _max_qubits: int + + def __init__(self, data: Data = BasicStates.PLUS, nqubit: int | None = None, max_qubits: int | None = None) -> None: + """Initialize a statevector object. `data` can be: - a single :class:`graphix.states.State` (classical description of a quantum state) - an iterable of :class:`graphix.states.State` objects - - an iterable of scalars (A 2**n numerical statevector) - - a *graphix.statevec.Statevec* object + - an iterable of scalars (a :math:`2^n` numerical statevector) + - a single :class:`graphix.sim.statevec.Statevec` - If *nqubit* is not provided, the number of qubit is inferred from *data* and checked for consistency. - If only one :class:`graphix.states.State` is provided and nqubit is a valid integer, initialize the statevector - in the tensor product state. - If both *nqubit* and *data* are provided, consistency of the dimensions is checked. - If a *graphix.statevec.Statevec* is passed, returns a copy. + If ``nqubit`` is not provided, it is inferred from ``data``. + If ``max_qubits`` is not provided, it is set to match the provided or inferred ``nqubit``. + If only one :class:`graphix.states.State` is provided and ``nqubit`` is a valid integer, the statevector is initialized in the tensor product state. + If a class:`graphix.sim.statevec.Statevec` is provided, a copy is returned. + Consistency between provided ``nqubit``, ``max_qubits`` and ``data`` is checked. Parameters ---------- data : Data, optional - input data to prepare the state. Can be a classical description or a numerical input, defaults to graphix.states.BasicStates.PLUS - nqubit : int, optional - number of qubits to prepare, defaults to None + Input data to prepare the state. Can be a classical description or a numerical input, defaults to :class:`graphix.states.BasicStates.PLUS`. + nqubit : int | None, optional + Number of qubits to prepare. If ``None`` (default), it's inferred from ``data``. + max_qubits : int | None, optional. + Maximum Hilbert space size for array preallocation. If ``None`` (default), it's set equal to ``nqubit``. + + Raises + ------ + ValueError + If ``nqubit``, ``max_qubits`` or ``data`` are not consistent with each other. """ if nqubit is not None and nqubit < 0: - raise ValueError("nqubit must be a non-negative integer.") + raise ValueError("`nqubit` must be a non-negative integer.") + + if max_qubits is not None and max_qubits < 0: + raise ValueError("`max_qubits` must be a non-negative integer.") if isinstance(data, Statevec): - # assert nqubit is None or len(state.flatten()) == 2**nqubit - if nqubit is not None and len(data.flatten()) != 2**nqubit: + if nqubit is not None and len(data.flatten()) != 1 << nqubit: raise ValueError( f"Inconsistent parameters between nqubit = {nqubit} and the inferred number of qubit = {len(data.flatten())}." ) - self.psi = data.psi.copy() + self._psi = data._psi.copy() + self._max_qubits = data.max_qubits + self._nqubit = data.nqubit + + if max_qubits is not None: + if max_qubits < data.max_qubits: + raise ValueError( + f"`max_qubits` can't be smaller than the capacity of input state: {max_qubits} < {data.max_qubits}." + ) + self.ensure_capacity(max_qubits) return # The type @@ -99,19 +152,20 @@ def __init__( elif isinstance(data, Iterable): input_list = list(data) else: - raise TypeError(f"Incorrect type for data: {type(data)}") + raise TypeError(f"Incorrect type for data: {type(data)}.") if len(input_list) == 0: if nqubit is not None and nqubit != 0: - raise ValueError("nqubit is not null but input state is empty.") - - self.psi = np.array(1, dtype=np.complex128) + raise ValueError("`nqubit` is not null but input state is empty.") + nqubit = 0 + psi = np.array([1], dtype=np.complex128) elif isinstance(input_list[0], states.State): + length = len(input_list) if nqubit is None: - nqubit = len(input_list) - elif nqubit != len(input_list): - raise ValueError("Mismatch between nqubit and length of input state.") + nqubit = length + elif nqubit != length: + raise ValueError(f"Mismatch between nqubit and length of input state: {nqubit} != {length}.") def state_to_statevector( s: states.State | ExpressionOrSupportsComplex | Iterable[ExpressionOrSupportsComplex], @@ -120,113 +174,260 @@ def state_to_statevector( raise TypeError("Data should be an homogeneous sequence of states.") return s.to_statevector() - list_of_sv = [state_to_statevector(s) for s in input_list] + psi = functools.reduce( + lambda m0, m1: np.kron(m0, m1).astype(np.complex128, copy=False), + (state_to_statevector(s) for s in input_list), + ) - tmp_psi = functools.reduce(lambda m0, m1: np.kron(m0, m1).astype(np.complex128), list_of_sv) - # reshape - self.psi = tmp_psi.reshape((2,) * nqubit) # `SupportsFloat` is needed because `numpy.float64` is not an instance of `SupportsComplex`! - elif isinstance(input_list[0], (Expression, SupportsComplex, SupportsFloat)): + elif isinstance(input_list[0], (SupportsComplex, SupportsFloat)): + length = len(input_list) + inferred_nqubit = length.bit_length() - 1 if nqubit is None: - length = len(input_list) if length & (length - 1): - raise ValueError("Length is not a power of two") - nqubit = length.bit_length() - 1 - elif nqubit != len(input_list).bit_length() - 1: - raise ValueError("Mismatch between nqubit and length of input state") - psi = np.array(input_list) - # check only if the matrix is not symbolic - if psi.dtype != "O" and not np.allclose(np.sqrt(np.sum(np.abs(psi) ** 2)), 1): - raise ValueError("Input state is not normalized") - self.psi = psi.reshape((2,) * nqubit) + raise ValueError(f"Length of input data is not a power of two: {length}.") + nqubit = inferred_nqubit + elif nqubit != inferred_nqubit: + raise ValueError(f"Mismatch between nqubit and inferred nqubit: {nqubit} != {inferred_nqubit}.") + psi = np.array(input_list, dtype=np.complex128) + if not np.isclose(np.linalg.norm(psi), 1.0): + raise ValueError("Input state is not normalized.") + else: - raise TypeError(f"First element of data has type {type(input_list[0])} whereas Number or State is expected") + raise TypeError( + f"First element of data has type {type(input_list[0])} whereas Number or State is expected." + ) + + if max_qubits is not None: + if max_qubits < nqubit: + raise ValueError( + f"`max_qubits` can't be smaller than the length of input state: {max_qubits} < {nqubit}." + ) + else: + max_qubits = nqubit + + self._psi = psi + self._max_qubits = nqubit # bootstrap for self.ensure_capacity + self._nqubit = nqubit + self.ensure_capacity(max_qubits) # may extend both self._psi and self._max_qubits def __str__(self) -> str: """Return a string description.""" - return f"Statevec object with statevector {self.psi} and length {self.dims()}." + sv = self.psi + return f"Statevec object with statevector {sv} and length {len(sv)}." + + @property + def psi(self) -> npt.NDArray[np.complex128]: + """Return a view of the meaningful elements in ``self._psi``. + + These are the first ``2**self.nqubit`` elements. + """ + size_valid_psi = 1 << self.nqubit # 2**self.nqubit + return self._psi[:size_valid_psi] + + # Note that `@property` must appear before `@override` for pyright + @property + @override + def nqubit(self) -> int: + """Return the number of qubits.""" + return self._nqubit + + @property + def max_qubits(self) -> int: + """Return the preallocated number of qubits.""" + return self._max_qubits + + def ensure_capacity(self, required_qubits: int) -> None: + """Extend the state vector if the required qubit capacity exceeds the current one. + + It copies the full vector state if ``required_qubits > self.max_qubits``, + otherwise does nothing. + + Parameters + ---------- + required_qubits : int + Minimum number of qubits the state vector must support. If expansion + is needed, ``self._psi`` is extended to size ``2**required_qubits``. + """ + if required_qubits > self.max_qubits: + offset = (1 << required_qubits) - len(self._psi) + self._psi = np.concatenate([self._psi, np.empty(offset, dtype=self._psi.dtype)]) + self._max_qubits = required_qubits + + @override + def flatten(self) -> Matrix: + """Return flattened state. + + A view of only the first ``2**self.nqubit`` elements of ``self._psi`` is returned. + """ + return self.psi @override def add_nodes(self, nqubit: int, data: Data) -> None: - r""" - Add nodes (qubits) to the state vector and initialize them in a specified state. + r"""Add nodes (qubits) to the state vector and initialize them in a specified state. + + Previously existing nodes remain unchanged. Parameters ---------- nqubit : int The number of qubits to add to the state vector. - data : Data, optional + data : Data The state in which to initialize the newly added nodes. - If a single basic state is provided, all new nodes are initialized in that state. - If a list of basic states is provided, it must match the length of ``nodes``, and each node is initialized with its corresponding state. - A single-qubit state vector will be broadcast to all nodes. - - A multi-qubit state vector of dimension :math:`2^n`, where :math:`n = \mathrm{len}(nodes)`, - initializes the new nodes jointly. + - A multi-qubit state vector of dimension :math:`2^n`, where :math:`n = \mathrm{len}(nodes)`, initializes the new nodes jointly. Notes ----- - Previously existing nodes remain unchanged. + This method can extend the size of ``self._psi`` for convenience, but this requires allocating a full new array. """ - sv_to_add = Statevec(nqubit=nqubit, data=data) - self.tensor(sv_to_add) + self.ensure_capacity(required_qubits=self.nqubit + nqubit) + if nqubit == 1 and data is BasicStates.PLUS: + # Simulating standard N commands falls in this branch. + _add_default_node_jit(self._psi, self.nqubit) + self._nqubit += 1 + else: + sv_to_add = Statevec(nqubit=nqubit, data=data) + self.tensor(sv_to_add) @override - def evolve_single(self, op: Matrix, i: int) -> None: - """Apply a single-qubit operation. + def entangle(self, qubits: tuple[int, int]) -> None: + """Apply a CZ gate on two qubits. Parameters ---------- - op : numpy.ndarray - 2*2 matrix - i : int - qubit index + qubits : tuple[int, int] + (control, target) qubit indices. """ - psi = tensordot(op, self.psi, (1, i)) - self.psi = np.moveaxis(psi, 0, i) + # `_entangle_jit` is not unsafe if called on out-of-bound indices but + # we check them for robustness. + for qubit in qubits: + self._check_bounds(qubit) + _entangle_jit(self._psi, self.nqubit, *qubits) @override - def evolve(self, op: Matrix, qargs: Sequence[int]) -> None: - """Apply a multi-qubit operation. + def evolve_single(self, op: Matrix, qubit: int) -> None: + """Apply a single-qubit operator. Parameters ---------- - op : numpy.ndarray - 2^n*2^n matrix - qargs : list of int - target qubits' indices + op : npt.NDArray[np.complex128] + Complex-valued matrix of shape :math:`(2, 2)` representing + the operator to apply. + qubit : int + Target qubit index. """ - op_dim = int(np.log2(len(op))) - # TODO shape = (2,)* 2 * op_dim - shape = [2 for _ in range(2 * op_dim)] - op_tensor = op.reshape(shape) - psi = tensordot( - op_tensor, - self.psi, - (tuple(op_dim + i for i in range(len(qargs))), qargs), - ) - self.psi = np.moveaxis(psi, range(len(qargs)), qargs) + self._check_bounds(qubit) + # Downcast from Matrix to npt.NDArray[np.complex128] to match numba signature. + op_as_complex = _cast_op(op) + _evolve_single_jit(self._psi, op_as_complex, self.nqubit, qubit) - def dims(self) -> tuple[int, ...]: - """Return the dimensions.""" - return self.psi.shape + @override + def expectation_single(self, op: Matrix, qubit: int) -> complex: + """Return the expectation value of a single-qubit operator. + + Parameters + ---------- + op : npt.NDArray[np.complex128] + Complex-valued matrix of shape :math:`(2, 2)` representing + the operator to measure. + qubit : int + Target qubit index. + + Returns + ------- + complex + Expectation value. + + Notes + ----- + This method assumes that quantum state represented by ``self.psi`` is normalized. + See the class docstring for details. + """ + self._check_bounds(qubit) + # Downcast from Matrix to npt.NDArray[np.complex128] to match numba signature. + op_as_complex = _cast_op(op) + return _expectation_single_jit(self._psi, op_as_complex, self.nqubit, qubit) - # Note that `@property` must appear before `@override` for pyright - @property @override - def nqubit(self) -> int: - """Return the number of qubits.""" - return self.psi.ndim + def evolve(self, op: Matrix, qubits: Sequence[int]) -> None: + r"""Apply a multi-qubit operator. + + Parameters + ---------- + op : npt.NDArray[np.complex128] + Complex-valued matrix of shape :math:`(2^n, 2^n)` representing + the operator to apply. + qubits : Sequence[int] + Target qubit indices. + + Notes + ----- + This method is a fallback for circuit simulation and it's not required + for pattern simulation. It does not have an efficient JIT-compiled + implementation. + """ + nq = len(qubits) + # treat op as a tensor with nq output + nq input legs + op_t = op.reshape((2,) * (nq * 2)).astype(np.complex128, copy=False) + psi_t = self.flatten().reshape((2,) * self.nqubit).astype(np.complex128, copy=False) + + psi_idx = np.array(range(self.nqubit)) + out_idx = np.array(range(self.nqubit, self.nqubit + nq)) # fresh labels + + op_idx = np.concatenate((out_idx, qubits)) + + # result subscripts: same as psi but modified indices (qargs) replaced by out labels + res_idx = psi_idx.copy() + for i, s in enumerate(qubits): + res_idx[s] = out_idx[i] + + self._psi[: len(self.psi)] = np.einsum(op_t, op_idx, psi_t, psi_idx, res_idx).reshape(1 << self.nqubit) # type: ignore[arg-type] # https://github.com/numpy/numpy/issues/31513 + + def expectation_value(self, op: Matrix, qubits: Sequence[int]) -> complex: + """Return the expectation value of a multi-qubit operator. + + Parameters + ---------- + op : npt.NDArray[np.complex128] + Complex-valued matrix of shape :math:`(2^n, 2^n)` representing + the operator to measure. + qubits : Sequence[int] + Target qubit indices. + + Notes + ----- + This method assumes that quantum state represented by ``self.psi`` is normalized. + See the class docstring for details. + """ + sv = deepcopy(self) + sv.evolve(op, qubits) + return complex(np.dot(self.flatten().conjugate(), sv.flatten())) @override - def remove_qubit(self, qarg: int) -> None: - r"""Remove a separable qubit from the system and assemble a statevector for remaining qubits. + def remove_qubit(self, qubit: int) -> None: + r"""Remove a separable qubit from the system and assemble the statevector of the remaining qubits. + + This method is deprecated. See :meth:`self.project_qubit`. - This results in the same result as partial trace, if the qubit *qarg* is separable from the rest. + Parameters + ---------- + qubit : int + Target qubit index. + """ + raise NotImplementedError("This method is deprecated. See :meth:`self.project_qubit`.") + + def project_qubit(self, op: Matrix, qubit: int) -> None: + r"""Project out a qubit from the system and assemble the statevector of the remaining qubits. + + This method evolves the statevector with ``op`` and removes ``qubit``. It assumes that after the application of ``op``, ``qubit`` is a separable qubit. - For a statevector :math:`\ket{\psi} = \sum c_i \ket{i}` with sum taken over + Let :math:`P \ket{\psi} = \sum c_i \ket{i}` be the statevector after the application of the projector :math:`P` (``op``), with sum taken over :math:`i \in [ 0 \dots 00,\ 0\dots 01,\ \dots,\ 1 \dots 11 ]`, this method returns @@ -244,153 +445,97 @@ def remove_qubit(self, qarg: int) -> None: \ket{1 \dots 1_{\mathrm{k-1}}1_{\mathrm{k+1}} \dots 11}, \end{align} - (after normalization) for :math:`k =` qarg. If the :math:`k` th qubit is in :math:`\ket{1}` state, - above will return zero amplitudes; in such a case the returned state will be the one above with - :math:`0_{\mathrm{k}}` replaced with :math:`1_{\mathrm{k}}` . + (after normalization) for :math:`k =` ``qubit``. If the :math:`k` th qubit is in the :math:`\ket{1}` state, all the amplitudes above will be zero. + In that case the returned state will be the one above with + :math:`0_{\mathrm{k}}` replaced with :math:`1_{\mathrm{k}}`. .. warning:: - This method assumes the qubit with index *qarg* to be separable from the rest, - and is implemented as a significantly faster alternative for partial trace to - be used after single-qubit measurements. - Care needs to be taken when using this method. - Checks for separability will be implemented soon as an option. + This method assumes that ``op`` is a projector, and therefore that ``qubit`` is a separable ``qubit`` after applying ``op``. It is + implemented as a significantly faster alternative for partial trace to + be used after single-qubit measurements. Separability is not checked. Parameters ---------- - qarg : int - qubit index - """ - norm = _norm(self.psi) - if isinstance(norm, SupportsFloat): - assert not np.isclose(norm, 0) - index: list[slice[int] | int] = [slice(None)] * self.psi.ndim - index[qarg] = 0 - psi = self.psi[tuple(index)] - norm = _norm(psi) - if isinstance(norm, SupportsFloat) and math.isclose(norm, 0): - index[qarg] = 1 - psi = self.psi[tuple(index)] - self.psi = psi - self.normalize() - - @override - def entangle(self, edge: tuple[int, int]) -> None: - """Connect graph nodes. + op : npt.NDArray[np.complex128] + Complex-valued matrix of shape :math:`(2, 2)` representing + the projector to apply. + qubit : int + Target qubit index. - Parameters - ---------- - edge : tuple of int - (control, target) qubit indices - """ - # contraction: 2nd index - control index, and 3rd index - target index. - psi = tensordot(CZ_TENSOR, self.psi, ((2, 3), edge)) - # sort back axes - self.psi = np.moveaxis(psi, (0, 1), edge) + Notes + ----- + This method combines :meth:`evolve_single` and :meth:`remove_qubit`. - def tensor(self, other: Statevec) -> None: - r"""Tensor product state with other qubits. + First, the relevant entries of ``self.psi`` are partitioned into two + disjoint subvectors. The partition depends on the value of ``q``. For a + three-qubit statevector ``[a0, a1, a2, a3, a4, a5, a6, a7]``: - Results in self :math:`\otimes` other. + :: - Parameters - ---------- - other : :class:`graphix.sim.statevec.Statevec` - statevector to be tensored with self - """ - psi_self = self.psi.flatten() - psi_other = other.psi.flatten() + q = 0: sv0 = [a0, a1, a2, a3], sv1 = [a4, a5, a6, a7] + q = 1: sv0 = [a0, a1, a4, a5], sv1 = [a2, a3, a6, a7] + q = 2: sv0 = [a0, a2, a4, a6], sv1 = [a1, a3, a5, a7] - total_num = len(self.dims()) + len(other.dims()) - self.psi = kron(psi_self, psi_other).reshape((2,) * total_num) + The statevector is then iterated over as in :meth:`evolve_single`, + updating ``self.psi`` in place while computing the norm of each + partition separately. - def cnot(self, qubits: tuple[int, int]) -> None: - """Apply CNOT. + If both resulting subvectors have nonzero norm, they are guaranteed to + be identical. The subvector ``sv0`` is normalized and written to the first + ``2**(nqubit - 1)`` entries of ``self.psi``. - Parameters - ---------- - qubits : tuple of int - (control, target) qubit indices + If exactly one subvector has nonzero norm, that subvector is selected. + If both have zero norm, a ``RuntimeError`` is raised. """ - # contraction: 2nd index - control index, and 3rd index - target index. - psi = tensordot(CNOT_TENSOR, self.psi, ((2, 3), qubits)) - # sort back axes - self.psi = np.moveaxis(psi, (0, 1), qubits) + self._check_bounds(qubit) + op_as_complex = _cast_op(op) + self._nqubit = _project_qubit_jit(self._psi, op_as_complex, self.nqubit, qubit, atol=1e-10) @override def swap(self, qubits: tuple[int, int]) -> None: - """Swap qubits. + """Apply SWAP gate between two qubits. Parameters ---------- - qubits : tuple of int - (control, target) qubit indices + qubits : tuple[int, int] + (control, target) qubit indices. """ - # contraction: 2nd index - control index, and 3rd index - target index. - psi = tensordot(SWAP_TENSOR, self.psi, ((2, 3), qubits)) - # sort back axes - self.psi = np.moveaxis(psi, (0, 1), qubits) - - def normalize(self) -> None: - """Normalize the state in-place.""" - # Note that the following calls to `astype` are guaranteed to - # return the original NumPy array itself, since `copy=False` and - # the `dtype` matches. This is important because the array is - # then modified in place. - if self.psi.dtype == np.object_: - psi_o = self.psi.astype(np.object_, copy=False) - norm_o = _norm_symbolic(psi_o) - psi_o /= norm_o - self.psi = psi_o - else: - psi_c = self.psi.astype(np.complex128, copy=False) - norm_c = _norm_numeric(psi_c) - psi_c /= norm_c - self.psi = psi_c + _swap_jit(self._psi, self.nqubit, *qubits) - def flatten(self) -> Matrix: - """Return flattened statevector.""" - return self.psi.flatten() + def tensor(self, other: Statevec) -> None: + r"""Tensor product state with other qubits. - @override - def expectation_single(self, op: Matrix, loc: int) -> complex: - """Return the expectation value of single-qubit operator. + Results in ``self`` :math:`\otimes` ``other``. Parameters ---------- - op : numpy.ndarray - 2*2 operator - loc : int - target qubit index + other : :class:`graphix.sim.statevec.Statevec` + Statevector to be tensored with ``self``. - Returns - ------- - complex : expectation value. + Notes + ----- + This method is used internally by :meth:`add_nodes`. """ - st1 = copy.copy(self) - st1.normalize() - st2 = copy.copy(st1) - st1.evolve_single(op, loc) - return complex(np.dot(st2.psi.flatten().conjugate(), st1.psi.flatten())) + _tensor_jit(self._psi, other.psi, self.nqubit, other.nqubit) + self._nqubit += other.nqubit + + def _check_bounds(self, qubit: int) -> None: + """Check if qubit index is valid. - def expectation_value(self, op: Matrix, qargs: Sequence[int]) -> complex: - """Return the expectation value of multi-qubit operator. + This check is necessary because there is no bounds checking in Numba. See + https://numba.pydata.org/numba-doc/dev/reference/pysemantics.html#bounds-checking Parameters ---------- - op : numpy.ndarray - 2^n*2^n operator - qargs : list of int - target qubit indices + qubit : int + Target qubit index. - Returns - ------- - complex : expectation value + Raises + ------ + IndexError """ - st2 = copy.copy(self) - st2.normalize() - st1 = copy.copy(st2) - st1.evolve(op, qargs) - return complex(np.dot(st2.psi.flatten().conjugate(), st1.psi.flatten())) + if not 0 <= qubit < self.nqubit: + raise IndexError(f"Qubit index {qubit} out of range [0, {self.nqubit - 1}]") def fidelity(self, other: Statevec) -> float: r"""Calculate the fidelity against another statevector. @@ -399,8 +544,8 @@ def fidelity(self, other: Statevec) -> float: Parameters ---------- - other : :class:`Statevec` - statevector to compare with + other : :class:`graphix.sim.statevec.Statevec` + Statevector to compare with. Returns ------- @@ -417,12 +562,12 @@ def isclose(self, other: Statevec, *, rtol: float = 1e-09, atol: float = 0.0) -> Parameters ---------- - other : :class:`Statevec` - statevector to compare with + other : :class:`graphix.sim.statevec.Statevec` + Statevector to compare with. rtol : float - relative tolerance for :func:`math.isclose` + Relative tolerance for :func:`math.isclose`. atol : float - absolute tolerance for :func:`math.isclose` + Absolute tolerance for :func:`math.isclose`. Returns ------- @@ -534,51 +679,452 @@ def _to_dict_map( return {_format_encoding(self.nqubit, i, encoding): amp for i, amp in zip(i_vals, amp_vals, strict=True)} - def subs(self, variable: Parameter, substitute: ExpressionOrSupportsFloat) -> Statevec: - """Return a copy of the state vector where all occurrences of the given variable in measurement angles are substituted by the given value.""" - result = Statevec() - result.psi = np.vectorize(lambda value: parameter.subs(value, variable, substitute))(self.psi) - return result - def xreplace(self, assignment: Mapping[Parameter, ExpressionOrSupportsFloat]) -> Statevec: - """Return a copy of the state vector where all occurrences of the given keys in measurement angles are substituted by the given values in parallel.""" - result = Statevec() - result.psi = np.vectorize(lambda value: parameter.xreplace(value, assignment))(self.psi) - return result +def _format_encoding(nqubit: int, i: int, encoding: _ENCODING) -> str: + """Format the i-th basis vector as a ket. + + See :meth:`Statevec.to_dict` for additional details. + """ + display_width = nqubit + output = f"{i:0{display_width}b}" + if encoding == "LSB": + return output[::-1] + return output @dataclass(frozen=True) class StatevectorBackend(DenseStateBackend[Statevec]): - """MBQC simulator with statevector method.""" + """Numba JIT-compiled MBQC statevector backend simulator based on Ref. [1]. - state: Statevec = dataclasses.field(init=False, default_factory=lambda: Statevec(nqubit=0)) + See Also + -------- + graphix.sim.base_backend.DenseStateBackend + Base class describing available parameters and shared behavior. + Notes + ----- + By default, the backend is initialized with a 0-dimensional statevector + (a scalar ``1``) and ``max_qubits = 0``. -def _norm_symbolic(psi: npt.NDArray[np.object_]) -> ExpressionOrFloat: - """Return norm of the state.""" - flat = psi.flatten() - return check_expression_or_float(np.sqrt(np.sum(flat.conj() * flat))) + The internal state representation can be expanded using + ``StatevectorBackend.add_nodes``, but this is inefficient since it + requires copying the full quantum state array. + To preallocate memory for a fixed system size, use + :meth:`StatevectorBackend.with_capacity`. -def _norm_numeric(psi: npt.NDArray[np.complex128]) -> float: - flat = psi.flatten() - norm_sq = np.sum(flat.conj() * flat) - assert math.isclose(norm_sq.imag, 0, abs_tol=1e-15) - return math.sqrt(norm_sq.real) + References + ---------- + [1] McGuffin, M. J., Robert J-M., and Ikeda K. "How to Write a Simulator for Quantum Circuits from Scratch: A Tutorial.", 2025 (arXiv:2506.08142). + """ + state: Statevec = dataclasses.field(init=True, default_factory=lambda: Statevec(nqubit=0)) -def _norm(psi: Matrix) -> ExpressionOrFloat: - """Return norm of the state.""" - # Narrow psi to concrete dtype - if psi.dtype == np.object_: - return _norm_symbolic(psi.astype(np.object_, copy=False)) - return _norm_numeric(psi.astype(np.complex128, copy=False)) + def __post_init__(self) -> None: + """Validate backend configuration. + Raises + ------ + ValueError + If ``symbolic`` is ``True``, since the statevector backend + does not support symbolic simulation. + """ + if self.symbolic: + raise ValueError( + "Statevector backend does not support `symbolic` simulation. Consider using backend in `graphix-symbolic` plugin." + ) -def _format_encoding(nqubit: int, i: int, encoding: _ENCODING) -> str: - """Format the i-th basis vector as a ket. See :meth:`Statevec.to_dict` for additional details.""" - display_width = nqubit - output = f"{i:0{display_width}b}" - if encoding == "LSB": - return output[::-1] - return output + @classmethod + def with_capacity( + cls, max_qubits: int, state: Statevec | None = None, **kwargs: Unpack[DenseStateBackendKwargs] + ) -> Self: + """Initialize the backend with preallocated statevector capacity. + + Parameters + ---------- + max_qubits : int + Maximum number of qubits supported by the statevector. For pattern simulation this corresponds to ``Pattern.max_space()``. + state: Statevec | None = None + Initial backend state. If ``None``, the backend is initialized + with a 0-dimensional statevector (scalar ``1``). + **kwargs + Options for :class:`graphix.sim.base_backend.DenseStateBackend`. See + :class:`graphix.sim.base_backend.DenseStateBackendKwargs`. + + Returns + ------- + Self + Backend instance with capacity for up to ``max_qubits`` qubits. + """ + state_init = ( + Statevec(nqubit=0, max_qubits=max_qubits) if state is None else Statevec(state, max_qubits=max_qubits) + ) + return cls(state_init, **kwargs) + + +@nb.njit("(c16[::1], c16[::1], int32, int32)") +def _tensor_jit( + psi: npt.NDArray[np.complex128], + psi_other: npt.NDArray[np.complex128], + nqubit: int, + nqubit_other: int, +) -> None: + """Tensor in-place two state vectors. + + This function assumes that ``psi`` has enough capacity to contain the tensor product. + """ + size_psi = 1 << nqubit + size_other = 1 << nqubit_other + # We update the elements of `psi` in-place. + # This requires starting the update for the last element of the new psi, `size_psi * size_other - 1` + k = size_psi * size_other - 1 + max_psi_idx = size_psi - 1 + max_other_idx = size_other - 1 + + for i in range(size_psi): + alpha_old = psi[max_psi_idx - i] + for j in range(size_other): + psi[k] = alpha_old * psi_other[max_other_idx - j] + k -= 1 + + +@nb.njit("(c16[::1], int32)") +def _add_default_node_jit( + psi: npt.NDArray[np.complex128], + nqubit: int, +) -> None: + r"""Tensor in-place a one-qubit |+> state. + + This function follows the same logic as :func:`_tensor_jit` but is specialized to ``psi_other`` being a :math:`(1, 1)/\sqrt{2}` state. + """ + read_idx = (1 << nqubit) - 1 # 2**nqubit - 1 + write_idx = (1 << nqubit + 1) - 1 # 2**(nqubit + 1) -1 + amp = 1 / np.sqrt(2) + + while read_idx >= 0: + v = amp * psi[read_idx] + + psi[write_idx] = v + psi[write_idx - 1] = v + + read_idx -= 1 + write_idx -= 2 + + +@nb.njit("(c16[::1], int32, int32, int32)") +def _swap_jit(psi: npt.NDArray[np.complex128], nqubit: int, q1: int, q2: int) -> None: + """Swap two qubits. + + This function is inspired from Ref. [1]. + """ + if q1 == q2: + return + size_sv = 1 << nqubit + mask_1 = 1 << nqubit - 1 - q1 + mask_2 = 1 << nqubit - 1 - q2 + mask = mask_1 | mask_2 + # `mask` is an integer number whose binary representation has 1s at positions `q1` and `q2` and 0s elsewhere. + + for i in range(size_sv): + # i & mask_1 = 2^(nqubit - 1 - q_1) if the binary representation of `i` has a 1 at position `q1` and 0 otherwise. + i_has_1_at_q1 = bool(i & mask_1) + i_has_1_at_q2 = bool(i & mask_2) + if i_has_1_at_q1 != i_has_1_at_q2: + # `j` has the same binary representation as `i` except for bits `q1` and `q2` which are flipped. + j = i ^ mask + if j > i: # Ensure we don't swap the same indices twice. + psi[j], psi[i] = psi[i], psi[j] + + +@nb.njit("(c16[::1], c16[:,:], int32, int32, int32)") +def _compute_op_psi( + psi: npt.NDArray[np.complex128], op: npt.NDArray[np.complex128], b0: int, j: int, size_half_block: int +) -> None: + """Update ``psi`` in place in step of :func:`_evolve_single_jit`.""" + i1 = b0 + j + i2 = i1 + size_half_block + psi1 = psi[i1] + psi2 = psi[i2] + psi[i1] = op[0, 0] * psi1 + op[0, 1] * psi2 + psi[i2] = op[1, 0] * psi1 + op[1, 1] * psi2 + + +@nb.njit("(c16[::1], c16[:,:], int32, int32)", parallel=True) +def _evolve_single_jit(psi: npt.NDArray[np.complex128], op: npt.NDArray[np.complex128], nqubit: int, q: int) -> None: + """Apply a single-qubit operator. + + This function is inspired from Ref. [1]. + + The kernel switches to a parallel implementation when the + number of qubits exceeds ``NUM_QUBIT_PARALLEL`` (module constant). + """ + nblocks = 1 << q + size_block = 1 << nqubit - q # 2**(nqubit - q) + size_half_block = ( + size_block >> 1 + ) # Left-to-right tensor product encoding (first qubit corresponds to most significant bit). For right-to-left encoding use `size_half_block = 1 << i` + + if nqubit > NUM_QUBIT_PARALLEL: + if nblocks > 1: + for b in nb.prange(nblocks): + # WARNING: setting `b0 += size_block` may result in a race condition if parallel loop. + b0 = size_block * b + for j in range(size_half_block): + _compute_op_psi(psi, op, b0, j, size_half_block) + else: + for j in nb.prange(size_half_block): + _compute_op_psi(psi, op, 0, j, size_half_block) + else: + b0 = 0 + for _ in range(nblocks): + for j in range(size_half_block): + _compute_op_psi(psi, op, b0, j, size_half_block) + b0 += size_block + + +@nb.njit("c16(c16[::1], c16[:,:], int32, int32, int32)") +def _compute_psic_op_psi( + psi: npt.NDArray[np.complex128], op: npt.NDArray[np.complex128], b0: int, j: int, size_half_block: int +) -> complex: + i1 = b0 + j + i2 = i1 + size_half_block + psi1 = psi[i1] + psi2 = psi[i2] + b1 = op[0, 0] * psi1 + op[0, 1] * psi2 + b2 = op[1, 0] * psi1 + op[1, 1] * psi2 + return psi1.conjugate() * b1 + psi2.conjugate() * b2 # type: ignore[no-any-return] + + +@nb.njit("c16(c16[::1], c16[:,:], int32, int32)", parallel=True) +def _expectation_single_jit( + psi: npt.NDArray[np.complex128], op: npt.NDArray[np.complex128], nqubit: int, q: int +) -> complex: + """Compute expectation value of single-qubit operator. + + This function applies ``op`` on ``psi`` in the same way as :func:`_evolve_single`. + + The kernel switches to a parallel implementation when the + number of qubits exceeds ``NUM_QUBIT_PARALLEL`` (module constant). + """ + nblocks = 1 << q + size_block = 1 << nqubit - q + size_half_block = size_block >> 1 + + result = 0.0 + 0.0j + + if nqubit > NUM_QUBIT_PARALLEL: + if nblocks > 1: + for b in nb.prange(nblocks): + # WARNING: setting `b0 += size_block` may result in a race condition if parallel loop. + b0 = size_block * b + for j in range(size_half_block): + result += _compute_psic_op_psi(psi, op, b0, j, size_half_block) + else: + for j in nb.prange(size_half_block): + result += _compute_psic_op_psi(psi, op, 0, j, size_half_block) + else: + b0 = 0 + for _ in range(nblocks): + for j in range(size_half_block): + result += _compute_psic_op_psi(psi, op, b0, j, size_half_block) + b0 += size_block + + return result + + +@nb.njit("(c16[::1], int32, int32, int32)", parallel=True) +def _entangle_jit(psi: npt.NDArray[np.complex128], nqubit: int, control: int, target: int) -> None: + """Apply CZ gate on two qubits. + + This function is inspired from Ref. [1]. + + The kernel switches to a parallel implementation when the + number of qubits exceeds ``NUM_QUBIT_PARALLEL`` (module constant). + """ + size_sv = 1 << nqubit + mask_control = 1 << nqubit - 1 - control + mask_target = 1 << nqubit - 1 - target + mask = mask_control | mask_target + # `mask` is an integer number whose binary representation has 1s at positions `control` and `target` and 0s elsewhere. + + if nqubit > NUM_QUBIT_PARALLEL: + for i in nb.prange(size_sv): + if mask & i == mask: + psi[i] = -psi[i] + else: + for i in range(size_sv): + if mask & i == mask: + psi[i] = -psi[i] + + +@nb.njit("types.UniTuple(f8, 2)(c16[::1], c16[:,:], int32, int32, int32)") +def _compute_op_psi_a2( + psi: npt.NDArray[np.complex128], op: npt.NDArray[np.complex128], b0: int, j: int, size_half_block: int +) -> tuple[float, float]: + """Update ``psi`` in place in step of :func:`_evolve_single_and_compute_norm_jit` and compute ``abs(psi)**2`` of updated elements.""" + # WARNING: Updating in-place a two-dimensional np.array for the the norm, e.g., + # `norm2_array[0] += norm2`` + # may result in a race condition. Return a tuple instead!! + i1 = b0 + j + i2 = i1 + size_half_block + psi1 = psi[i1] + psi2 = psi[i2] + a1 = op[0, 0] * psi1 + op[0, 1] * psi2 + a2 = op[1, 0] * psi1 + op[1, 1] * psi2 + psi[i1] = a1 + psi[i2] = a2 + + a1_re = a1.real + a1_im = a1.imag + norm2_1 = a1_re * a1_re + a1_im * a1_im + + a2_re = a2.real + a2_im = a2.imag + norm2_2 = a2_re * a2_re + a2_im * a2_im + + return norm2_1, norm2_2 + + +@nb.njit("types.UniTuple(f8, 2)(c16[::1], c16[:,:], int32, int32, int32, int32)", parallel=True) +def _evolve_single_and_compute_norm( + psi: npt.NDArray[np.complex128], + op: npt.NDArray[np.complex128], + nqubit: int, + n_blocks: int, + size_block: int, + size_half_block: int, +) -> tuple[float, float]: + """Evolve the full statevector with ``op`` and compute the norm of ``psi[b0+j]``. + + The kernel switches to a parallel implementation when the + number of qubits exceeds ``NUM_QUBIT_PARALLEL`` (module constant). + """ + norm2_0 = 0.0 + norm2_1 = 0.0 + if nqubit > NUM_QUBIT_PARALLEL: + if n_blocks > 1: + for b in nb.prange(n_blocks): + b0 = b * size_block + for j in range(size_half_block): + norm2 = _compute_op_psi_a2(psi, op, b0, j, size_half_block) + norm2_0 += norm2[0] + norm2_1 += norm2[1] + else: + for j in nb.prange(size_half_block): + norm2 = _compute_op_psi_a2(psi, op, 0, j, size_half_block) + norm2_0 += norm2[0] + norm2_1 += norm2[1] + else: + for b in range(n_blocks): + b0 = b * size_block + for j in range(size_half_block): + norm2 = _compute_op_psi_a2(psi, op, b0, j, size_half_block) + norm2_0 += norm2[0] + norm2_1 += norm2[1] + return norm2_0, norm2_1 + + +@nb.njit("void(c16[::1], int32, int32, int32, int32, f8)") +def _scale_psi_kernel( + psi: npt.NDArray[np.complex128], + n_blocks: int, + size_block: int, + size_half_block: int, + shift: int, + inv_norm: float, +) -> None: + """Update ``psi`` with selected and normalized elements. + + The implementation of this function does not support a parallelized + kernel because data is read and written on the same array. + """ + b0 = shift + k = 0 + for _ in range(n_blocks): + for j in range(size_half_block): + psi[k] = psi[b0 + j] * inv_norm + k += 1 + b0 += size_block + + +@nb.njit("void(c16[::1], int32, int32, int32, int32, f8)") +def _scale_psi( + psi: npt.NDArray[np.complex128], + n_blocks: int, + size_block: int, + size_half_block: int, + shift: int, + inv_norm: float, +) -> None: + """Update ``psi`` with selected and normalized elements.""" + # If the inner loop in `_scale_psi_kernel` has too few elements + # it introduces some overhead and it's best to unroll it. + # Numba can do that only if ``size_half_block`` is a constant at + # compile time, so we dispatch the relevant cases. + # We observe speed improvements of up to 10%. + if size_block == 2: + _scale_psi_kernel(psi, n_blocks, 2, 1, shift, inv_norm) + elif size_block == 4: + _scale_psi_kernel(psi, n_blocks, 4, 2, shift, inv_norm) + elif size_block == 8: + _scale_psi_kernel(psi, n_blocks, 8, 4, shift, inv_norm) + else: + _scale_psi_kernel(psi, n_blocks, size_block, size_half_block, shift, inv_norm) + + +@nb.njit("int32(c16[::1], c16[:,:], int32, int32, f8)") +def _project_qubit_jit( + psi: npt.NDArray[np.complex128], + op: npt.NDArray[np.complex128], + nqubit: int, + q: int, + atol: float, +) -> int: + """Evolve the statevector with ``op`` and remove qubit ``q``. + + Argument ``atol`` controls the tolerance below which norm of statevector is 0. + See :meth:`Statevec.remove_qubit` for additional details on the implementation. + + This implementation benefited from Jérôme Richard's advice. + https://stackoverflow.com/questions/79948374/improving-efficiency-of-numba-jit-function + """ + new_nqubit = nqubit - 1 + + n_blocks = 1 << q + size_block = 1 << nqubit - q # 2**(nqubits - q) + size_half_block = size_block >> 1 + + # Apply `op` to the full statevector and compute norm of branch 0 + shift = 0 + + norm2, norm2_1 = _evolve_single_and_compute_norm(psi, op, nqubit, n_blocks, size_block, size_half_block) + + # If norm of branch 0 is 0, we pick norm of branch 1 and set shift to branch 1. + if norm2 <= atol: + shift = size_half_block + norm2 = norm2_1 + + if norm2 <= atol: + raise RuntimeError(f"Attempted to remove qubit {q} from 0-norm statevector.") + + inv_norm = 1.0 / math.sqrt(norm2) + _scale_psi(psi, n_blocks, size_block, size_half_block, shift, inv_norm) + + return new_nqubit + + +def _cast_op(op: Matrix) -> npt.NDArray[np.complex128]: + if op.dtype == np.object_: + raise TypeError( + "Statevector backend does not support symbolic operators. Consider using backend in `graphix-symbolic` plugin." + ) + # By default, the numba signature c16[:,:] assumes a writeable array. + # Arrays obtained from, e.g., `graphix.clifford.Clifford.H.matrix` or + # `graphix.ops.Ops.X` are not writtable and therefore do not match the + # numba signature. + # Since ``op`` is a 2-by-2 matrix it's probably not worth it to dispatch + # multiple jit kernels (with the appropriate numba signature) + # depending on its WRITEABLE flag to avoid copying. + # https://github.com/numba/numba/issues/4511#issuecomment-527350694 + # https://numba.pydata.org/numba-doc/0.17.0/reference/types.html#numba.types.Array + return op.astype(np.complex128, copy=True) diff --git a/graphix/simulator.py b/graphix/simulator.py index e81439ebc..12a64b9a7 100644 --- a/graphix/simulator.py +++ b/graphix/simulator.py @@ -356,7 +356,7 @@ def __init__( 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. + [Density matrix backend only] If True, support arbitrary objects (typically, symbolic expressions) in measurement angles. stacklevel : int, optional Stack level to use for warnings. Defaults to 1, meaning that warnings are reported at this function's call site. @@ -531,7 +531,7 @@ def _initialize_backend( 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. + [Density matrix backend only] If True, support arbitrary objects (typically, symbolic expressions) in measurement angles. Returns ------- @@ -562,7 +562,12 @@ def _initialize_backend( case "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 symbolic: + raise ValueError( + "Statevector backend does not support `symbolic` simulation. Consider using backend in `graphix-symbolic` plugin." + ) + nqubits = pattern.max_space() + return StatevectorBackend.with_capacity(nqubits, branch_selector=branch_selector) case "densitymatrix": if noise_model is None: warnings.warn( diff --git a/graphix/transpiler.py b/graphix/transpiler.py index 7859d3f0f..e4fd98ee0 100644 --- a/graphix/transpiler.py +++ b/graphix/transpiler.py @@ -1019,7 +1019,7 @@ def simulate_statevector( result : :class:`SimulateResult` output state of the statevector simulation and results of classical measures. """ - _backend = _initialize_backend(backend, branch_selector) + _backend = _initialize_backend(backend, branch_selector, self.width) if input_state is None: _backend.add_nodes(range(self.width)) @@ -1207,6 +1207,7 @@ def transpile_swaps(circuit: Circuit) -> TranspileSwapsResult: def _initialize_backend( backend: StatevectorBackend | Literal["statevector"], branch_selector: BranchSelector | None, + width: int, ) -> StatevectorBackend: ... @@ -1214,6 +1215,7 @@ def _initialize_backend( def _initialize_backend( backend: DensityMatrixBackend | Literal["densitymatrix"], branch_selector: BranchSelector | None, + width: int, ) -> DensityMatrixBackend: ... @@ -1221,12 +1223,14 @@ def _initialize_backend( def _initialize_backend( backend: DenseStateBackend[_DenseStateT_co], branch_selector: BranchSelector | None, + width: int, ) -> DenseStateBackend[_DenseStateT_co]: ... def _initialize_backend( backend: DenseStateBackend[_DenseStateT_co] | _DenseStateBackendLiteral, branch_selector: BranchSelector | None, + width: int, ) -> _BuiltinDenseStateBackend | DenseStateBackend[_DenseStateT_co]: """Initialize backend for circuit simulation. @@ -1236,6 +1240,9 @@ def _initialize_backend( Simulation backend branch_selector: :class:`BranchSelector` Branch selector used for measurements. Can only be specified if ``backend`` is not an already instantiated :class:`Backend` object. If ``None``, it defaults to :class:`RandomBranchSelector`. + width : int + Number of qubits in circuit. It is required to initialize the :class:`StatevectorBackend` with the appropriate + capacity. Returns ------- @@ -1252,7 +1259,7 @@ def _initialize_backend( match backend: case "statevector": - return StatevectorBackend(branch_selector=branch_selector) + return StatevectorBackend.with_capacity(width, branch_selector=branch_selector) case "densitymatrix": return DensityMatrixBackend(branch_selector=branch_selector) case _: diff --git a/noxfile.py b/noxfile.py index 65419b640..641d9c6fe 100644 --- a/noxfile.py +++ b/noxfile.py @@ -113,7 +113,7 @@ class ReverseDependency: ReverseDependency( "https://github.com/thierry-martinez/graphix-stim-backend", branch="fix/graphix_498_remove_pauli" ), - ReverseDependency("https://github.com/TeamGraphix/graphix-symbolic"), + ReverseDependency("https://github.com/matulni/graphix-symbolic", branch="backend"), ReverseDependency("https://github.com/TeamGraphix/graphix-qasm-parser"), ReverseDependency( "https://github.com/thierry-martinez/veriphix", diff --git a/tests/test_density_matrix.py b/tests/test_density_matrix.py index 7c67909c5..44db6fc92 100644 --- a/tests/test_density_matrix.py +++ b/tests/test_density_matrix.py @@ -15,8 +15,8 @@ from graphix.channels import KrausChannel, dephasing_channel, depolarising_channel from graphix.fundamentals import ANGLE_PI, Plane from graphix.ops import Ops -from graphix.sim.density_matrix import DensityMatrix, DensityMatrixBackend -from graphix.sim.statevec import CNOT_TENSOR, CZ_TENSOR, SWAP_TENSOR, Statevec +from graphix.sim.density_matrix import CNOT_TENSOR, CZ_TENSOR, SWAP_TENSOR, DensityMatrix, DensityMatrixBackend +from graphix.sim.statevec import Statevec from graphix.simulator import DefaultMeasureMethod from graphix.states import BasicStates, PlanarState from graphix.transpiler import Circuit diff --git a/tests/test_graphsim.py b/tests/test_graphsim.py index 180c38929..a10c636e3 100644 --- a/tests/test_graphsim.py +++ b/tests/test_graphsim.py @@ -83,23 +83,17 @@ def test_fig2(self) -> None: g = GraphState(nodes=np.arange(nqubit), edges=edges) gstate = graph_state_to_statevec(g) g.measure_x(0) - gstate.evolve_single(meas_op(0), 0) # x meas - gstate.normalize() - gstate.remove_qubit(0) + gstate.project_qubit(meas_op(0), 0) # x meas gstate2 = graph_state_to_statevec(g) assert gstate.isclose(gstate2) g.measure_y(1, choice=0) - gstate.evolve_single(meas_op(0.5 * ANGLE_PI), 0) # y meas - gstate.normalize() - gstate.remove_qubit(0) + gstate.project_qubit(meas_op(0.5 * ANGLE_PI), 0) # y meas gstate2 = graph_state_to_statevec(g) assert gstate.isclose(gstate2) g.measure_z(3) - gstate.evolve_single(meas_op(0.5 * ANGLE_PI, plane=Plane.YZ), 1) # z meas - gstate.normalize() - gstate.remove_qubit(1) + gstate.project_qubit(meas_op(0.5 * ANGLE_PI, plane=Plane.YZ), 1) # z meas gstate2 = graph_state_to_statevec(g) assert gstate.isclose(gstate2) diff --git a/tests/test_parameter.py b/tests/test_parameter.py index 88f68d07e..d8b2160dd 100644 --- a/tests/test_parameter.py +++ b/tests/test_parameter.py @@ -12,7 +12,6 @@ from graphix.pattern import DrawPatternAnnotations, Pattern from graphix.random_objects import rand_circuit from graphix.sim.density_matrix import DensityMatrix -from graphix.sim.statevec import Statevec if TYPE_CHECKING: from numpy.random import PCG64, Generator @@ -132,19 +131,6 @@ def test_parallel_substitution_with_zero() -> None: assert not pattern23.is_parameterized() -def test_statevec_subs() -> None: - alpha = Placeholder("alpha") - statevec = Statevec([alpha]) - assert np.allclose(statevec.subs(alpha, 1).psi, np.array([1])) - - -def test_statevec_xreplace() -> None: - alpha = Placeholder("alpha") - beta = Placeholder("beta") - statevec = Statevec([alpha, beta]) - assert np.allclose(statevec.xreplace({alpha: 1, beta: 2}).psi, np.array([1, 2])) - - def test_density_matrix_subs() -> None: alpha = Placeholder("alpha") dm = DensityMatrix([[alpha]]) diff --git a/tests/test_pattern.py b/tests/test_pattern.py index 6a87c4016..57c0589cc 100644 --- a/tests/test_pattern.py +++ b/tests/test_pattern.py @@ -170,7 +170,7 @@ def simulate_and_measure() -> int: sim.run(rng=fx_rng) state = sim.backend.state if isinstance(state, Statevec): - assert state.dims() == () + assert state.nqubit == 0 elif isinstance(state, DensityMatrix): assert state.dims() == (1, 1) elif isinstance(state, MBQCTensorNet): diff --git a/tests/test_statevec.py b/tests/test_statevec.py index e7ce8d925..3030250bc 100644 --- a/tests/test_statevec.py +++ b/tests/test_statevec.py @@ -4,213 +4,398 @@ from typing import TYPE_CHECKING import numpy as np +import numpy.typing as npt import pytest - -from graphix.fundamentals import ANGLE_PI, Plane -from graphix.pattern import Pattern -from graphix.sim.statevec import Statevec, _norm_numeric -from graphix.states import BasicStates, PlanarState +from numpy.random import Generator + +from graphix.branch_selector import ConstBranchSelector +from graphix.clifford import Clifford +from graphix.instruction import Instruction +from graphix.measurements import Measurement +from graphix.ops import Ops +from graphix.pauli import Pauli +from graphix.random_objects import rand_unit +from graphix.sim.statevec import NUM_QUBIT_PARALLEL, Statevec, StatevectorBackend +from graphix.states import BasicStates +from graphix.transpiler import Circuit if TYPE_CHECKING: from collections.abc import Mapping from typing import Literal - from numpy.random import Generator + from numpy.random import PCG64 _ENCODING = Literal["LSB", "MSB"] + from graphix.states import PlanarState, State -class TestStatevec: - """Test for Statevec class. Particularly new constructor.""" - - # test initializing one qubit in plus state - def test_default_success(self) -> None: - vec = Statevec(nqubit=1) - assert np.allclose(vec.psi, np.array([1, 1] / np.sqrt(2))) - assert len(vec.dims()) == 1 - - def test_basicstates_success(self) -> None: - # minus - vec = Statevec(nqubit=1, data=BasicStates.MINUS) - assert np.allclose(vec.psi, np.array([1, -1] / np.sqrt(2))) - assert len(vec.dims()) == 1 - - # zero - vec = Statevec(nqubit=1, data=BasicStates.ZERO) - assert np.allclose(vec.psi, np.array([1, 0]), rtol=0, atol=1e-15) - assert len(vec.dims()) == 1 - - # one - vec = Statevec(nqubit=1, data=BasicStates.ONE) - assert np.allclose(vec.psi, np.array([0, 1]), rtol=0, atol=1e-15) - assert len(vec.dims()) == 1 - - # plus_i - vec = Statevec(nqubit=1, data=BasicStates.PLUS_I) - assert np.allclose(vec.psi, np.array([1, 1j] / np.sqrt(2))) - assert len(vec.dims()) == 1 - - # minus_i - vec = Statevec(nqubit=1, data=BasicStates.MINUS_I) - assert np.allclose(vec.psi, np.array([1, -1j] / np.sqrt(2))) - assert len(vec.dims()) == 1 - - # even more tests? - def test_default_tensor_success(self, fx_rng: Generator) -> None: - nqb = int(fx_rng.integers(2, 5)) - print(f"nqb is {nqb}") - vec = Statevec(nqubit=nqb) - assert np.allclose(vec.psi, np.ones((2,) * nqb) / (np.sqrt(2)) ** nqb) - assert len(vec.dims()) == nqb - - vec = Statevec(nqubit=nqb, data=BasicStates.MINUS_I) - sv_list = [BasicStates.MINUS_I.to_statevector() for _ in range(nqb)] - sv = functools.reduce(lambda a, b: np.kron(a, b).astype(np.complex128, copy=False), sv_list) - assert np.allclose(vec.psi, sv.reshape((2,) * nqb)) - assert len(vec.dims()) == nqb - - # tensor of same state - rand_angle = fx_rng.random() * 2 * ANGLE_PI - rand_plane = fx_rng.choice(np.array(Plane)) - state = PlanarState(rand_plane, rand_angle) - vec = Statevec(nqubit=nqb, data=state) - sv_list = [state.to_statevector() for _ in range(nqb)] - sv = functools.reduce(lambda a, b: np.kron(a, b).astype(np.complex128, copy=False), sv_list) - assert np.allclose(vec.psi, sv.reshape((2,) * nqb)) - assert len(vec.dims()) == nqb - - # tensor of different states - rand_angles = fx_rng.random(nqb) * 2 * ANGLE_PI - rand_planes = fx_rng.choice(np.array(Plane), nqb) - states = [PlanarState(plane=i, angle=j) for i, j in zip(rand_planes, rand_angles, strict=True)] - vec = Statevec(nqubit=nqb, data=states) - sv_list = [state.to_statevector() for state in states] - sv = functools.reduce(lambda a, b: np.kron(a, b).astype(np.complex128, copy=False), sv_list) - assert np.allclose(vec.psi, sv.reshape((2,) * nqb)) - assert len(vec.dims()) == nqb - - def test_data_success(self, fx_rng: Generator) -> None: - nqb = fx_rng.integers(2, 5) - length = 2**nqb - rand_vec = fx_rng.random(length) + 1j * fx_rng.random(length) - rand_vec /= np.sqrt(np.sum(np.abs(rand_vec) ** 2)) - vec = Statevec(data=rand_vec) - assert np.allclose(vec.psi, rand_vec.reshape((2,) * nqb)) - assert len(vec.dims()) == nqb +N_JUMPS = 3 + + +def generate_rnd_data(rng: Generator, nqubits: int) -> npt.NDArray[np.complex128]: + length = 1 << nqubits + data = rng.random(length) + 1j * rng.random(length) + data /= np.sqrt(np.sum(np.abs(data) ** 2)) + return data - # fail: incorrect len - def test_data_dim_fail(self, fx_rng: Generator) -> None: - length = 5 + +class TestStatevec: + @pytest.mark.parametrize( + ("state", "data_ref"), + [ + (BasicStates.PLUS, np.array([1, 1] / np.sqrt(2))), + (BasicStates.MINUS, np.array([1, -1] / np.sqrt(2))), + (BasicStates.ZERO, np.array([1, 0])), + (BasicStates.ONE, np.array([0, 1])), + (BasicStates.PLUS_I, np.array([1, 1j] / np.sqrt(2))), + (BasicStates.MINUS_I, np.array([1, -1j] / np.sqrt(2))), + ], + ) + def test_init_basic_states(self, state: State, data_ref: npt.NDArray[np.complex128]) -> None: + sv = Statevec(data=state) + assert np.allclose(sv.flatten(), data_ref) + + @pytest.mark.parametrize("nqubit", range(5)) + def test_init_random_state(self, fx_rng: Generator, nqubit: int) -> None: + data = generate_rnd_data(fx_rng, nqubit) + sv = Statevec(data) + assert np.allclose(sv.flatten(), data) + + def test_init_preallocation(self) -> None: + nqubit = 2 + max_qubits = 5 + sv = Statevec(data=[BasicStates.PLUS, BasicStates.ZERO], nqubit=nqubit, max_qubits=max_qubits) + + assert np.allclose(sv.flatten(), np.array([1, 0, 1, 0]) / np.sqrt(2)) + assert sv.nqubit == nqubit + assert sv.max_qubits == max_qubits + assert len(sv._psi) == 2**max_qubits + + @pytest.mark.parametrize("nqubit", range(5)) + def test_init_statevec(self, fx_rng: Generator, nqubit: int) -> None: + data = generate_rnd_data(fx_rng, nqubit) + sv_1 = Statevec(data) + sv_2 = Statevec(sv_1) + assert np.allclose(sv_1.psi, sv_2.psi) + assert sv_1.nqubit == sv_2.nqubit + assert sv_1.max_qubits == sv_2.max_qubits + + @pytest.mark.parametrize("nqubit", range(5)) + def test_init_statevec_max_qubits(self, fx_rng: Generator, nqubit: int) -> None: + data = generate_rnd_data(fx_rng, nqubit) + sv_1 = Statevec(data) + sv_2 = Statevec(data=sv_1, max_qubits=nqubit + 1) + assert np.allclose(sv_1.flatten(), sv_2.flatten()) + assert sv_1.nqubit == sv_2.nqubit + assert sv_1.max_qubits + 1 == sv_2.max_qubits + + # fail: length data is not a power of 2 + @pytest.mark.parametrize("length", [3, 5, 6, 7]) + def test_init_dim_fail(self, fx_rng: Generator, length: int) -> None: rand_vec = fx_rng.random(length) + 1j * fx_rng.random(length) rand_vec /= np.sqrt(np.sum(np.abs(rand_vec) ** 2)) with pytest.raises(ValueError): - _vec = Statevec(data=rand_vec) + Statevec(data=rand_vec) - # fail: with less qubit than number of qubits inferred from a correct state vect - def test_data_dim_fail_mismatch(self, fx_rng: Generator) -> None: - nqb = 3 - rand_vec = fx_rng.random(2**nqb) + 1j * fx_rng.random(2**nqb) - rand_vec /= np.sqrt(np.sum(np.abs(rand_vec) ** 2)) + # fail: different qubit than number of qubits inferred from data + @pytest.mark.parametrize("nqubit", [2, 4]) + def test_init_dim_mismatch_fail(self, fx_rng: Generator, nqubit: int) -> None: + data = generate_rnd_data(fx_rng, nqubit) with pytest.raises(ValueError): - _vec = Statevec(nqubit=2, data=rand_vec) + Statevec(nqubit=3, data=data) # fail: not normalized - def test_data_norm_fail(self, fx_rng: Generator) -> None: - nqb = fx_rng.integers(2, 5) - length = 2**nqb - rand_vec = fx_rng.random(length) + 1j * fx_rng.random(length) + def test_init_norm_fail(self, fx_rng: Generator) -> None: + data = 5 * generate_rnd_data(fx_rng, 3) with pytest.raises(ValueError): - _vec = Statevec(data=rand_vec) + Statevec(data=data) - def test_defaults_to_one(self) -> None: - vec = Statevec() - assert len(vec.dims()) == 1 + # fail: max qubits smaller than number of qubits + def test_init_max_qubits_fail(self) -> None: + nqubit = 4 + max_qubits = 3 + with pytest.raises(ValueError): + Statevec(nqubit=nqubit, max_qubits=max_qubits) - # try copying Statevec input - def test_copy_success(self, fx_rng: Generator) -> None: - nqb = fx_rng.integers(2, 5) - length = 2**nqb - rand_vec = fx_rng.random(length) + 1j * fx_rng.random(length) - rand_vec /= np.sqrt(np.sum(np.abs(rand_vec) ** 2)) - test_vec = Statevec(data=rand_vec) - # try to copy it - vec = Statevec(data=test_vec) + # fail: incorrect number of qubits + @pytest.mark.parametrize("nqubit", range(5)) + def test_init_statevec_fail(self, fx_rng: Generator, nqubit: int) -> None: + data = generate_rnd_data(fx_rng, nqubit) + sv_1 = Statevec(data) + with pytest.raises(ValueError): + Statevec(data=sv_1, nqubit=nqubit + 1) - assert np.allclose(vec.psi, test_vec.psi) - assert len(vec.dims()) == len(test_vec.dims()) + # fail: incorrect max_qubits qubits + @pytest.mark.parametrize("nqubit", range(5)) + def test_init_statevec_max_qubits_fail(self, fx_rng: Generator, nqubit: int) -> None: + data = generate_rnd_data(fx_rng, nqubit) + sv_1 = Statevec(data=data, max_qubits=6) + with pytest.raises(ValueError): + Statevec(data=sv_1, max_qubits=nqubit + 1) - # try calling with incorrect number of qubits compared to inferred one - def test_copy_fail(self, fx_rng: Generator) -> None: - nqb = int(fx_rng.integers(2, 5)) - length = 1 << nqb - rand_vec = fx_rng.random(length) + 1j * fx_rng.random(length) - rand_vec /= np.sqrt(np.sum(np.abs(rand_vec) ** 2)) - test_vec = Statevec(data=rand_vec) + @pytest.mark.parametrize( + ("sv", "edge", "data_ref"), + [ + (Statevec(data=BasicStates.ZERO, nqubit=2), (0, 1), np.array([1, 0, 0, 0])), + (Statevec(data=[BasicStates.PLUS, BasicStates.PLUS]), (0, 1), np.array([1, 1, 1, -1]) / 2), + (Statevec(data=[BasicStates.ONE, BasicStates.MINUS]), (0, 1), np.array([0, 0, 1, 1]) / np.sqrt(2)), + ( + Statevec(data=np.array([1, 0, 0, 0, 0, 0, 0, 1]) / np.sqrt(2)), + (0, 2), + np.array([1, 0, 0, 0, 0, 0, 0, -1]) / np.sqrt(2), + ), + ], + ) + def test_entangle(self, sv: Statevec, edge: tuple[int, int], data_ref: npt.NDArray[np.complex128]) -> None: + sv.entangle(edge) + assert np.allclose(sv.flatten(), data_ref) - with pytest.raises(ValueError): - _vec = Statevec(nqubit=length - 1, data=test_vec) - - def test_nqubits(self) -> None: - for i in [1, 2, 5]: - sv = Statevec(nqubit=i) - assert sv.nqubit == i - - def test_nqubits_pattern(self) -> None: - p = Pattern(input_nodes=[0, 1, 2]) - sv = p.simulate_pattern(backend="statevector") - assert sv.nqubit == 3 - - -class TestFidelityIsclose: - def test_fidelity_same_state(self) -> None: - state = Statevec(data=BasicStates.PLUS) - assert state.fidelity(state) == pytest.approx(1) - - def test_fidelity_orthogonal(self) -> None: - zero = Statevec(data=BasicStates.ZERO) - one = Statevec(data=BasicStates.ONE) - assert zero.fidelity(one) == pytest.approx(0) - - def test_fidelity_known_value(self) -> None: - # F(|0>, |+>) = 0.5 - zero = Statevec(data=BasicStates.ZERO) - plus = Statevec(data=BasicStates.PLUS) - assert zero.fidelity(plus) == pytest.approx(0.5) - - def test_fidelity_global_phase(self) -> None: - plus = Statevec(data=BasicStates.PLUS) - plus_rotated = Statevec(data=np.array([1, 1]) / np.sqrt(2) * 1j) - assert plus.fidelity(plus_rotated) == pytest.approx(1) - - def test_fidelity_symmetry(self, fx_rng: Generator) -> None: - length = 4 - vec_a = fx_rng.random(length) + 1j * fx_rng.random(length) - vec_a /= np.sqrt(np.sum(np.abs(vec_a) ** 2)) - vec_b = fx_rng.random(length) + 1j * fx_rng.random(length) - vec_b /= np.sqrt(np.sum(np.abs(vec_b) ** 2)) - a = Statevec(data=vec_a) - b = Statevec(data=vec_b) - assert a.fidelity(b) == pytest.approx(b.fidelity(a)) - - def test_isclose_same_state(self) -> None: - state = Statevec(data=BasicStates.PLUS) - assert state.isclose(state) - - def test_isclose_orthogonal(self) -> None: - zero = Statevec(data=BasicStates.ZERO) - one = Statevec(data=BasicStates.ONE) - assert not zero.isclose(one) - - def test_isclose_global_phase(self) -> None: - plus = Statevec(data=BasicStates.PLUS) - rotated = Statevec(data=np.array([1, 1]) / np.sqrt(2) * np.exp(1j * 0.7)) - assert plus.isclose(rotated) - - def test_isclose_tolerance(self) -> None: - zero = Statevec(data=BasicStates.ZERO) - almost = Statevec(data=np.array([np.sqrt(1 - 1e-8), np.sqrt(1e-8)])) - assert not zero.isclose(almost) - assert zero.isclose(almost, atol=1e-6) + @pytest.mark.parametrize( + ("sv", "q", "op", "data_ref"), + [ + (Statevec(data=BasicStates.ZERO, nqubit=2), 0, Clifford.X.matrix, np.array([0, 0, 1, 0])), + ( + Statevec(data=[BasicStates.PLUS, BasicStates.PLUS]), + 1, + Clifford.H.matrix, + np.array([1, 0, 1, 0]) / np.sqrt(2), + ), + ( + Statevec(data=[BasicStates.PLUS, BasicStates.MINUS]), + 0, + np.array([[1, 0], [0, np.exp(0.25j * np.pi)]]), + np.array([1, -1, np.exp(0.25j * np.pi), -np.exp(0.25j * np.pi)]) / 2, + ), + ( + Statevec(data=np.array([1, 0, 0, 0, 0, 0, 0, 1]) / np.sqrt(2)), + 1, + Clifford.Z.matrix, + np.array([1, 0, 0, 0, 0, 0, 0, -1]) / np.sqrt(2), + ), + ], + ) + def test_evolve_single( + self, sv: Statevec, q: int, op: npt.NDArray[np.complex128], data_ref: npt.NDArray[np.complex128] + ) -> None: + sv.evolve_single(op, q) + assert np.allclose(sv.flatten(), data_ref) + + @pytest.mark.parametrize( + ("sv", "q", "op", "exp_ref"), + [ + (Statevec(data=BasicStates.ZERO, nqubit=2), 0, Clifford.X.matrix, 0), + (Statevec(data=[BasicStates.PLUS, BasicStates.PLUS]), 1, Clifford.H.matrix, 1 / np.sqrt(2)), + ( + Statevec(data=[BasicStates.PLUS, BasicStates.MINUS]), + 0, + np.array([[1, 0], [0, np.exp(0.25j * np.pi)]]), + (1 + np.exp(0.25j * np.pi)) / 2, + ), + ( + Statevec(data=np.array([1, 0, 0, 0, 0, 0, 0, 1]) / np.sqrt(2)), + 1, + Clifford.Z.matrix, + 0, + ), + ], + ) + def test_expectation_single( + self, sv: Statevec, q: int, op: npt.NDArray[np.complex128], exp_ref: np.complex128 + ) -> None: + assert np.isclose(sv.expectation_single(op, q), exp_ref) + + def test_add_nodes(self, fx_rng: Generator) -> None: + max_qubits = 5 + sv_test = Statevec(nqubit=0, max_qubits=max_qubits) + psi_ref = np.array([1.0 + 0.0j]) + + for _ in range(max_qubits): # Add a node at each iteration + data = generate_rnd_data(fx_rng, nqubits=1) + psi_ref = np.kron(psi_ref, data) + sv_test.add_nodes(1, data) + assert np.allclose(sv_test.flatten(), psi_ref) + + def test_add_nodes_beyond_max_qubits(self, fx_rng: Generator) -> None: + max_qubits = 3 + sv_test = Statevec(nqubit=max_qubits, max_qubits=max_qubits) + + # We create a reference to `psi` to ensure that array is not extended with np.ndarray.resize + psi_ref = sv_test.psi # noqa: F841 + + nqubits_new = 2 + data = generate_rnd_data(fx_rng, nqubits=nqubits_new) + psi = np.kron(sv_test.psi, data) + + sv_test.add_nodes(nqubits_new, data) + assert np.allclose(sv_test.flatten(), psi) + + @pytest.mark.parametrize( + ("sv", "q", "sv_ref"), + [ + (Statevec(data=BasicStates.ZERO, nqubit=2), 0, Statevec(data=BasicStates.ZERO, nqubit=1)), + (Statevec(data=[BasicStates.PLUS, BasicStates.PLUS]), 1, Statevec(data=BasicStates.PLUS, nqubit=1)), + (Statevec(data=[BasicStates.PLUS, BasicStates.MINUS]), 0, Statevec(data=BasicStates.MINUS, nqubit=1)), + (Statevec(data=[BasicStates.ZERO, BasicStates.ONE]), 0, Statevec(data=BasicStates.ONE, nqubit=1)), + # In previous testcase, branch 1 is 0 (psi_10 == psi_11 == 0), and first element of branch 0 is 0 too (psi_00 == 0)! + ( + Statevec(data=[BasicStates.PLUS_I, BasicStates.ONE, BasicStates.PLUS]), + 1, + Statevec(data=[BasicStates.PLUS_I, BasicStates.PLUS], nqubit=2), + ), + ], + ) + def test_project_qubit(self, sv: Statevec, q: int, sv_ref: Statevec) -> None: + # This test mimics the behavior of former `remove_qubit`. + op = np.eye(2, dtype=np.complex128) + sv.project_qubit(op, q) + assert np.linalg.norm(sv.psi) == pytest.approx(1) + assert np.allclose(sv.flatten(), sv_ref.flatten()) + + @pytest.mark.parametrize( + "state", + [ + BasicStates.PLUS, + BasicStates.MINUS, + BasicStates.ZERO, + BasicStates.ONE, + BasicStates.PLUS_I, + BasicStates.MINUS_I, + ], + ) + def test_measurement_into_each_xyz_basis(self, state: PlanarState) -> None: + n = 3 + k = 0 + statevector = state.to_statevector() + m_op = np.outer(statevector, statevector.T.conjugate()) + sv = Statevec(nqubit=n) + sv.evolve_single(m_op, k) + + if state is BasicStates.MINUS: + # Measurement into |-> results in a 0-norm vector + with pytest.raises(RuntimeError): + sv.project_qubit(m_op, k) + else: + sv.project_qubit(m_op, k) + sv2 = Statevec(nqubit=n - 1) + assert sv.isclose(sv2) + + @pytest.mark.parametrize( + ("sv", "qargs", "op", "data_ref"), + [ + (Statevec(data=BasicStates.ZERO, nqubit=2), (0,), Clifford.X.matrix, np.array([0, 0, 1, 0])), + ( + Statevec(data=[BasicStates.PLUS, BasicStates.PLUS]), + (1,), + Clifford.H.matrix, + np.array([1, 0, 1, 0]) / np.sqrt(2), + ), + ( + Statevec(data=[BasicStates.PLUS, BasicStates.MINUS]), + (0,), + np.array([[1, 0], [0, np.exp(0.25j * np.pi)]]), + np.array([1, -1, np.exp(0.25j * np.pi), -np.exp(0.25j * np.pi)]) / 2, + ), + ( + Statevec(data=np.array([1, 0, 0, 0, 0, 0, 0, 1]) / np.sqrt(2)), + (1,), + Clifford.Z.matrix, + np.array([1, 0, 0, 0, 0, 0, 0, -1]) / np.sqrt(2), + ), + ( + Statevec(data=np.array([1, 0, 0, 0, 0, 0, 0, 1]) / np.sqrt(2)), + (0, 1), + Ops.CNOT, + np.array([1, 0, 0, 0, 0, 1, 0, 0]) / np.sqrt(2), + ), + ( + Statevec(data=np.array([1, 0, 0, 0, 0, 0, 0, 1]) / np.sqrt(2)), + (0, 2), + Ops.CNOT, + np.array([1, 0, 0, 0, 0, 0, 1, 0]) / np.sqrt(2), + ), + ( + Statevec(data=np.array([1, 0, 0, 0, 0, 0, 0, 1]) / np.sqrt(2), max_qubits=5), + (1, 2), + Ops.CNOT, + np.array([1, 0, 0, 0, 0, 0, 1, 0]) / np.sqrt(2), + ), + ( + Statevec(data=np.array([1, 0, 0, 0, 0, 0, 0, 1]) / np.sqrt(2), max_qubits=5), + (0, 1, 2), + Ops.CCX, + np.array([1, 0, 0, 0, 0, 0, 1, 0]) / np.sqrt(2), + ), + ], + ) + def test_evolve( + self, sv: Statevec, qargs: tuple[int, ...], op: npt.NDArray[np.complex128], data_ref: npt.NDArray[np.complex128] + ) -> None: + sv.evolve(op, qargs) + assert np.allclose(sv.flatten(), data_ref) + + @pytest.mark.parametrize("jumps", range(1, N_JUMPS)) + def test_evolve_rnd(self, fx_bg: PCG64, jumps: int) -> None: + rng = Generator(fx_bg.jumped(jumps)) + nqubits = 4 + data = generate_rnd_data(rng, nqubits) + sv = Statevec(data=data) + op = rand_unit(8, rng) + + sv.evolve(op, (1, 2, 3)) + data_ref = np.kron(np.eye(2), op) @ data + sv_ref = Statevec(data_ref) + + assert sv.isclose(sv_ref) + + @pytest.mark.parametrize("jumps", range(1, N_JUMPS)) + def test_expectation_value(self, fx_bg: PCG64, jumps: int) -> None: + rng = Generator(fx_bg.jumped(jumps)) + nqubits = 4 + data = generate_rnd_data(rng, nqubits) + sv = Statevec(data=data) + op = rand_unit(4, rng) + + val_test = sv.expectation_value(op, (1, 2)) + val_ref = np.conjugate(data) @ functools.reduce(np.kron, (np.eye(2), op, np.eye(2))) @ data + assert val_test == pytest.approx(val_ref) + + @pytest.mark.parametrize( + ("sv1", "sv2", "fidelity"), + [ + (Statevec(data=BasicStates.PLUS), Statevec(data=BasicStates.PLUS), 1), + (Statevec(data=BasicStates.ZERO), Statevec(data=BasicStates.ONE), 0), + (Statevec(data=BasicStates.ZERO), Statevec(data=BasicStates.PLUS), 0.5), + (Statevec(data=BasicStates.PLUS), Statevec(data=np.array([1, 1]) / np.sqrt(2) * 1j), 1), + ], + ) + def test_fidelity(self, sv1: Statevec, sv2: Statevec, fidelity: float) -> None: + assert sv1.fidelity(sv2) == pytest.approx(fidelity) + + @pytest.mark.parametrize( + ("sv1", "sv2", "isclose", "atol"), + [ + (Statevec(data=BasicStates.PLUS), Statevec(data=BasicStates.PLUS), True, 0.0), + (Statevec(data=BasicStates.ZERO), Statevec(data=BasicStates.ONE), False, 0.0), + ( + Statevec(data=BasicStates.PLUS), + Statevec(data=np.array([1, 1]) / np.sqrt(2) * np.exp(1j * 0.7)), + True, + 0.0, + ), + (Statevec(data=BasicStates.ZERO), Statevec(data=np.array([np.sqrt(1 - 1e-8), np.sqrt(1e-8)])), False, 0.0), + ( + Statevec(data=BasicStates.ZERO), + Statevec(data=np.array([np.sqrt(1 - 1e-8), np.sqrt(1e-8)])), + True, + 1.0e-6, + ), + ], + ) + def test_isclose(self, sv1: Statevec, sv2: Statevec, isclose: bool, atol: float) -> None: + if isclose: + assert sv1.isclose(sv2, atol=atol) + else: + assert not sv1.isclose(sv2, atol=atol) @pytest.mark.parametrize( ("encoding", "dict_ref"), @@ -238,8 +423,201 @@ def test_to_prob_dict(self, encoding: _ENCODING, dict_ref: Mapping[str, float]) assert np.isclose(dict_ref[ket], amp2.real) assert np.isclose(0, amp2.imag) + # Run with `nqubits` larger and smaller than ``graphix.sim.statevec.NUM_QUBIT_PARALLEL`` + # to test parallelized and non-parallelized kernels. + # This usually allows to detect race conditions. + @pytest.mark.parametrize("nqubits", [3, NUM_QUBIT_PARALLEL + 1]) + def test_simulation_identity(self, fx_rng: Generator, nqubits: int) -> None: + qc = Circuit(nqubits) + angle_1 = 2 * fx_rng.random() + angle_2 = 2 * fx_rng.random() + qc.extend( + [ + Instruction.H(0), + Instruction.CNOT(0, 1), + Instruction.CNOT(1, 2), + Instruction.RZZ(0, 1, angle_1), + Instruction.RY(2, angle_2), + ] + ) + qc.extend( + [ + Instruction.RZZ(0, 1, -angle_1), + Instruction.RY(2, -angle_2), + Instruction.CNOT(1, 2), + Instruction.CNOT(0, 1), + Instruction.H(0), + ] + ) + pattern = qc.transpile().pattern + data = generate_rnd_data(fx_rng, nqubits) + sv = Statevec(data=data) + sv_test = pattern.simulate_pattern(backend="statevector", input_state=sv, rng=fx_rng) + + assert sv_test.isclose(sv) + + +class TestStatevectorBackend: + @pytest.mark.parametrize( + ("state", "data_ref"), + [ + (BasicStates.PLUS, np.array([1, 1] / np.sqrt(2))), + (BasicStates.MINUS, np.array([1, -1] / np.sqrt(2))), + (BasicStates.ZERO, np.array([1, 0])), + (BasicStates.ONE, np.array([0, 1])), + (BasicStates.PLUS_I, np.array([1, 1j] / np.sqrt(2))), + (BasicStates.MINUS_I, np.array([1, -1j] / np.sqrt(2))), + ], + ) + def test_init_basic_states(self, state: State, data_ref: npt.NDArray[np.complex128]) -> None: + backend = StatevectorBackend() + backend.add_nodes([0], data=state) + sv = Statevec(data=data_ref) + assert backend.state.isclose(sv) + + @pytest.mark.parametrize( + ("n_nodes"), + range(5), + ) + def test_init_capacity(self, n_nodes: int, fx_rng: Generator) -> None: + capacity = 3 + data = generate_rnd_data(fx_rng, n_nodes) + backend = StatevectorBackend.with_capacity(capacity) + backend.add_nodes(range(n_nodes), data=data) + sv = Statevec(data=data) + assert backend.state.isclose(sv) + assert backend.state.max_qubits == max(n_nodes, capacity) + + def test_init_branch_selector(self) -> None: + max_qubits = 5 + sv = Statevec(nqubit=3) + bs = ConstBranchSelector(result=0) + + backend = StatevectorBackend.with_capacity(max_qubits, sv, branch_selector=bs) + + assert backend.branch_selector is bs + assert backend.state.isclose(sv) + assert backend.state.max_qubits == 5 + + def test_init_branch_selector_kw(self) -> None: + sv = Statevec(nqubit=3) + bs = ConstBranchSelector(result=0) -def test_normalize() -> None: - statevec = Statevec(nqubit=1, data=BasicStates.PLUS) - statevec.remove_qubit(0) - assert _norm_numeric(statevec.psi.astype(np.complex128, copy=False)) == 1 + with pytest.raises(TypeError): + # Branch selector is keyword only (enforced by dataclass sentinel) + StatevectorBackend(sv, bs) # type: ignore[misc, arg-type] + + backend = StatevectorBackend(sv, branch_selector=bs) + + assert backend.branch_selector is bs + assert backend.state.isclose(sv) + assert backend.state.max_qubits == 3 + + @pytest.mark.parametrize( + ("clifford"), + Clifford, + ) + def test_clifford(self, clifford: Clifford, fx_rng: Generator) -> None: + nqubits = 4 + data = generate_rnd_data(fx_rng, nqubits=nqubits) + + backend = StatevectorBackend() + backend.add_nodes(nodes=range(nqubits), data=data) + backend.apply_clifford(node=0, clifford=clifford) + + vec = Statevec(nqubit=nqubits, data=data) + vec.evolve_single(clifford.matrix, 0) + + assert backend.state.isclose(vec) + + @pytest.mark.parametrize("jumps", range(1, N_JUMPS)) + def test_deterministic_measure_0(self, fx_bg: PCG64, jumps: int) -> None: + rng = Generator(fx_bg.jumped(jumps)) + # plus state & zero state (default), but with tossed coins + + backend = StatevectorBackend() + coins = [rng.choice([0, 1]), rng.choice([0, 1])] + expected_result = sum(coins) % 2 + states = [ + Pauli.X.eigenstate(coins[0]), + Pauli.Z.eigenstate(coins[1]), + ] + nodes = range(len(states)) + backend.add_nodes(nodes=nodes, data=states) + backend.entangle_nodes(edge=(nodes[0], nodes[1])) + measurement = Measurement.X + node_to_measure = backend.node_index[0] + result = backend.measure(node=node_to_measure, measurement=measurement, rng=rng) + assert result == expected_result + + @pytest.mark.parametrize("jumps", range(1, N_JUMPS)) + def test_deterministic_measure_1(self, fx_bg: PCG64, jumps: int) -> None: + """Entangle |+> state with N |0> states, the (XY,0) measurement yields the outcome 0 with probability 1.""" + rng = Generator(fx_bg.jumped(jumps)) + n_nodes = 11 + backend = StatevectorBackend.with_capacity(n_nodes) + states = [BasicStates.PLUS, *(BasicStates.ZERO for _ in range(n_nodes - 1))] + backend.add_nodes(nodes=range(n_nodes), data=states) + for i in range(1, n_nodes): + backend.entangle_nodes(edge=(0, i)) + measurement = Measurement.X + node_to_measure = backend.node_index[0] + result = backend.measure(node=node_to_measure, measurement=measurement, rng=rng) + assert result == 0 + assert list(backend.node_index) == list(range(1, n_nodes)) + + @pytest.mark.parametrize("jumps", range(1, N_JUMPS)) + def test_deterministic_measure_many(self, fx_bg: PCG64, jumps: int) -> None: + """Entangle |+> state with N |0> states, the (XY,0) measurement yields the outcome 0 with probability 1.""" + rng = Generator(fx_bg.jumped(jumps)) + # plus state (default) + backend = StatevectorBackend() + n_traps = 5 + n_neighbors = 5 + n_others = 5 + traps = [Pauli.X.eigenstate() for _ in range(n_traps)] + dummies = [Pauli.Z.eigenstate() for _ in range(n_neighbors)] + others = [Pauli.I.eigenstate() for _ in range(n_others)] + states = traps + dummies + others + nodes = range(len(states)) + backend.add_nodes(nodes=nodes, data=states) + + for dummy in nodes[n_traps : n_traps + n_neighbors]: + for trap in nodes[:n_traps]: + backend.entangle_nodes(edge=(trap, dummy)) + for other in nodes[n_traps + n_neighbors :]: + backend.entangle_nodes(edge=(other, dummy)) + + # Same measurement for all traps + measurement = Measurement.X + + for trap in nodes[:n_traps]: + node_to_measure = trap + result = backend.measure(node=node_to_measure, measurement=measurement, rng=rng) + assert result == 0 + + assert list(backend.node_index) == list(range(n_traps, n_neighbors + n_traps + n_others)) + + @pytest.mark.parametrize("jumps", range(1, N_JUMPS)) + def test_deterministic_measure_with_coin(self, fx_bg: PCG64, jumps: int) -> None: + """Entangle |+> state with N |0> states, the (XY,0) measurement yields the outcome 0 with probability 1. + + We add coin toss to that. + """ + rng = Generator(fx_bg.jumped(jumps)) + # plus state (default) + backend = StatevectorBackend() + n_neighbors = 10 + coins = [rng.choice([0, 1])] + [rng.choice([0, 1]) for _ in range(n_neighbors)] + expected_result = sum(coins) % 2 + states = [Pauli.X.eigenstate(coins[0])] + [Pauli.Z.eigenstate(coins[i + 1]) for i in range(n_neighbors)] + nodes = range(len(states)) + backend.add_nodes(nodes=nodes, data=states) + + for i in range(1, n_neighbors + 1): + backend.entangle_nodes(edge=(nodes[0], i)) + measurement = Measurement.X + node_to_measure = backend.node_index[0] + result = backend.measure(node=node_to_measure, measurement=measurement, rng=rng) + assert result == expected_result + assert list(backend.node_index) == list(range(1, n_neighbors + 1)) diff --git a/tests/test_statevec_backend.py b/tests/test_statevec_backend.py deleted file mode 100644 index acab2c527..000000000 --- a/tests/test_statevec_backend.py +++ /dev/null @@ -1,189 +0,0 @@ -from __future__ import annotations - -from typing import TYPE_CHECKING - -import numpy as np -import pytest - -from graphix.clifford import Clifford -from graphix.fundamentals import ANGLE_PI, Plane -from graphix.measurements import Measurement -from graphix.pauli import Pauli -from graphix.sim.statevec import Statevec, StatevectorBackend -from graphix.states import BasicStates, PlanarState - -if TYPE_CHECKING: - from numpy.random import Generator - - from graphix import Pattern - - -class TestStatevec: - @pytest.mark.parametrize( - "state", [BasicStates.PLUS, BasicStates.ZERO, BasicStates.ONE, BasicStates.PLUS_I, BasicStates.MINUS_I] - ) - def test_measurement_into_each_xyz_basis(self, state: PlanarState) -> None: - n = 3 - k = 0 - # for measurement into |-> returns [[0, 0], ..., [0, 0]] (whose norm is zero) - statevector = state.to_statevector() - m_op = np.outer(statevector, statevector.T.conjugate()) - sv = Statevec(nqubit=n) - sv.evolve(m_op.astype(np.complex128, copy=False), [k]) - sv.remove_qubit(k) - - sv2 = Statevec(nqubit=n - 1) - assert sv.isclose(sv2) - - def test_measurement_into_minus_state(self) -> None: - n = 3 - k = 0 - m_op = np.outer(BasicStates.MINUS.to_statevector(), BasicStates.MINUS.to_statevector().T.conjugate()) - sv = Statevec(nqubit=n) - sv.evolve(m_op.astype(np.complex128, copy=False), [k]) - with pytest.raises(AssertionError): - sv.remove_qubit(k) - - -class TestStatevecNew: - # test initialization only - def test_init_success(self, hadamardpattern: Pattern, fx_rng: Generator) -> None: - # plus state (default) - backend = StatevectorBackend() - backend.add_nodes(hadamardpattern.input_nodes) - vec = Statevec(nqubit=1) - assert np.allclose(vec.psi, backend.state.psi) - assert len(backend.state.dims()) == 1 - - # minus state - backend = StatevectorBackend() - backend.add_nodes(hadamardpattern.input_nodes, data=BasicStates.MINUS) - vec = Statevec(nqubit=1, data=BasicStates.MINUS) - assert np.allclose(vec.psi, backend.state.psi) - assert len(backend.state.dims()) == 1 - - # random planar state - rand_angle = fx_rng.random() * 2 * ANGLE_PI - rand_plane = fx_rng.choice(np.array(Plane)) - state = PlanarState(rand_plane, rand_angle) - backend = StatevectorBackend() - backend.add_nodes(hadamardpattern.input_nodes, data=state) - vec = Statevec(nqubit=1, data=state) - assert np.allclose(vec.psi, backend.state.psi) - # assert backend.state.nqubit == 1 - assert len(backend.state.dims()) == 1 - - # data input and Statevec input - - def test_init_fail(self, hadamardpattern: Pattern, fx_rng: Generator) -> None: - rand_angle = fx_rng.random(2) * 2 * ANGLE_PI - rand_plane = fx_rng.choice(np.array(Plane), 2) - - state = PlanarState(rand_plane[0], rand_angle[0]) - state2 = PlanarState(rand_plane[1], rand_angle[1]) - with pytest.raises(ValueError): - StatevectorBackend().add_nodes(hadamardpattern.input_nodes, data=[state, state2]) - - def test_clifford(self) -> None: - for clifford in Clifford: - state = BasicStates.PLUS - vec = Statevec(nqubit=1, data=state) - backend = StatevectorBackend() - backend.add_nodes(nodes=[0], data=state) - # Applies clifford gate "Z" - vec.evolve_single(clifford.matrix, 0) - backend.apply_clifford(node=0, clifford=clifford) - np.testing.assert_allclose( - vec.psi.astype(np.complex128, copy=False), backend.state.psi.astype(np.complex128, copy=False) - ) - - def test_deterministic_measure_one(self, fx_rng: Generator) -> None: - # plus state & zero state (default), but with tossed coins - for _ in range(10): - backend = StatevectorBackend() - coins = [fx_rng.choice([0, 1]), fx_rng.choice([0, 1])] - expected_result = sum(coins) % 2 - states = [ - Pauli.X.eigenstate(coins[0]), - Pauli.Z.eigenstate(coins[1]), - ] - nodes = range(len(states)) - backend.add_nodes(nodes=nodes, data=states) - - backend.entangle_nodes(edge=(nodes[0], nodes[1])) - measurement = Measurement.X - node_to_measure = backend.node_index[0] - result = backend.measure(node=node_to_measure, measurement=measurement, rng=fx_rng) - assert result == expected_result - - def test_deterministic_measure(self, fx_rng: Generator) -> None: - """Entangle |+> state with N |0> states, the (XY,0) measurement yields the outcome 0 with probability 1.""" - for _ in range(10): - # plus state (default) - backend = StatevectorBackend() - n_neighbors = 10 - states = [Pauli.X.eigenstate()] + [Pauli.Z.eigenstate() for i in range(n_neighbors)] - nodes = range(len(states)) - backend.add_nodes(nodes=nodes, data=states) - - for i in range(1, n_neighbors + 1): - backend.entangle_nodes(edge=(nodes[0], i)) - measurement = Measurement.X - node_to_measure = backend.node_index[0] - result = backend.measure(node=node_to_measure, measurement=measurement, rng=fx_rng) - assert result == 0 - assert list(backend.node_index) == list(range(1, n_neighbors + 1)) - - def test_deterministic_measure_many(self, fx_rng: Generator) -> None: - """Entangle |+> state with N |0> states, the (XY,0) measurement yields the outcome 0 with probability 1.""" - for _ in range(10): - # plus state (default) - backend = StatevectorBackend() - n_traps = 5 - n_neighbors = 5 - n_whatever = 5 - traps = [Pauli.X.eigenstate() for _ in range(n_traps)] - dummies = [Pauli.Z.eigenstate() for _ in range(n_neighbors)] - others = [Pauli.I.eigenstate() for _ in range(n_whatever)] - states = traps + dummies + others - nodes = range(len(states)) - backend.add_nodes(nodes=nodes, data=states) - - for dummy in nodes[n_traps : n_traps + n_neighbors]: - for trap in nodes[:n_traps]: - backend.entangle_nodes(edge=(trap, dummy)) - for other in nodes[n_traps + n_neighbors :]: - backend.entangle_nodes(edge=(other, dummy)) - - # Same measurement for all traps - measurement = Measurement.X - - for trap in nodes[:n_traps]: - node_to_measure = trap - result = backend.measure(node=node_to_measure, measurement=measurement, rng=fx_rng) - assert result == 0 - - assert list(backend.node_index) == list(range(n_traps, n_neighbors + n_traps + n_whatever)) - - def test_deterministic_measure_with_coin(self, fx_rng: Generator) -> None: - """Entangle |+> state with N |0> states, the (XY,0) measurement yields the outcome 0 with probability 1. - - We add coin toss to that. - """ - for _ in range(10): - # plus state (default) - backend = StatevectorBackend() - n_neighbors = 10 - coins = [fx_rng.choice([0, 1])] + [fx_rng.choice([0, 1]) for _ in range(n_neighbors)] - expected_result = sum(coins) % 2 - states = [Pauli.X.eigenstate(coins[0])] + [Pauli.Z.eigenstate(coins[i + 1]) for i in range(n_neighbors)] - nodes = range(len(states)) - backend.add_nodes(nodes=nodes, data=states) - - for i in range(1, n_neighbors + 1): - backend.entangle_nodes(edge=(nodes[0], i)) - measurement = Measurement.X - node_to_measure = backend.node_index[0] - result = backend.measure(node=node_to_measure, measurement=measurement, rng=fx_rng) - assert result == expected_result - assert list(backend.node_index) == list(range(1, n_neighbors + 1)) diff --git a/tests/test_tnsim.py b/tests/test_tnsim.py index c347e492b..bd4d23e13 100644 --- a/tests/test_tnsim.py +++ b/tests/test_tnsim.py @@ -12,7 +12,7 @@ from graphix.command import E from graphix.fundamentals import ANGLE_PI from graphix.ops import Ops -from graphix.random_objects import rand_circuit +from graphix.random_objects import rand_circuit, rand_unit from graphix.sim.statevec import Statevec from graphix.sim.tensornet import MBQCTensorNet, gen_str from graphix.states import BasicStates @@ -399,8 +399,8 @@ def test_evolve(self, fx_bg: PCG64, jumps: int, fx_rng: Generator) -> None: pattern.remove_pauli_measurements() state = circuit.simulate_statevector().statevec tn_mbqc = pattern.simulate_pattern(backend="tensornetwork", rng=fx_rng) - random_op3 = random_op(3, rng) - random_op3_exp = random_op(3, rng) + random_op3 = rand_unit(2**3, rng) + random_op3_exp = rand_unit(2**3, rng) state.evolve(random_op3, [0, 1, 2]) tn_mbqc.evolve(random_op3, [0, 1, 2], decompose=False) diff --git a/tests/test_transpiler.py b/tests/test_transpiler.py index f848a5ed2..219c2240f 100644 --- a/tests/test_transpiler.py +++ b/tests/test_transpiler.py @@ -304,8 +304,9 @@ def test_transpile_swaps(fx_bg: PCG64, jumps: int) -> None: for qubit in transpiled_swaps.qubits: assert qubit is not None qubits.append(qubit) - state2.psi = np.transpose(state2.psi, qubits) - assert state.isclose(state2) + psi_t = np.transpose(state2.flatten().reshape((2,) * nqubits), qubits) + state2_test = Statevec(psi_t.reshape(1 << nqubits)) + assert state.isclose(state2_test) @pytest.mark.parametrize("jumps", range(1, 11))