diff --git a/.gitignore b/.gitignore index 6c68c78..9e5ead9 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,5 @@ *~ -__pycache__/ \ No newline at end of file +__pycache__/ +*.egg-info/ +./*/_version.py +/graphix_symbolic/_version.py \ No newline at end of file diff --git a/graphix_symbolic/__init__.py b/graphix_symbolic/__init__.py index 4593f1b..7d0f92c 100644 --- a/graphix_symbolic/__init__.py +++ b/graphix_symbolic/__init__.py @@ -1,3 +1,5 @@ +from graphix_symbolic.density_matrix import DensityMatrix, DensityMatrixBackend +from graphix_symbolic.statevec import Statevec, StatevectorBackend from graphix_symbolic.sympy_parameter import SympyParameter -__all__ = ["SympyParameter"] +__all__ = ["DensityMatrix", "DensityMatrixBackend", "Statevec", "StatevectorBackend", "SympyParameter"] diff --git a/graphix_symbolic/_version.py b/graphix_symbolic/_version.py new file mode 100644 index 0000000..ab0a2ab --- /dev/null +++ b/graphix_symbolic/_version.py @@ -0,0 +1,24 @@ +# file generated by vcs-versioning +# don't change, don't track in version control +from __future__ import annotations + +__all__ = [ + "__version__", + "__version_tuple__", + "version", + "version_tuple", + "__commit_id__", + "commit_id", +] + +version: str +__version__: str +__version_tuple__: tuple[int | str, ...] +version_tuple: tuple[int | str, ...] +commit_id: str | None +__commit_id__: str | None + +__version__ = version = '0.1.dev12+ge4c27e99e.d20260601' +__version_tuple__ = version_tuple = (0, 1, 'dev12', 'ge4c27e99e.d20260601') + +__commit_id__ = commit_id = 'ge4c27e99e' diff --git a/graphix_symbolic/density_matrix.py b/graphix_symbolic/density_matrix.py new file mode 100644 index 0000000..060c10e --- /dev/null +++ b/graphix_symbolic/density_matrix.py @@ -0,0 +1,425 @@ +"""Density matrix simulator. + +Simulate MBQC with density matrix representation. +""" + +from __future__ import annotations + +import copy +import dataclasses +import math +from collections.abc import Collection, Iterable +from dataclasses import dataclass +from typing import TYPE_CHECKING + +import numpy as np +from graphix import linalg_validations as lv +from graphix import parameter +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.states import BasicStates, State +from typing_extensions import override + +from graphix_symbolic.statevec import CNOT_TENSOR, CZ_TENSOR, SWAP_TENSOR, Statevec + +if TYPE_CHECKING: + from collections.abc import Mapping, Sequence + from typing import SupportsComplex, SupportsFloat + + from graphix.noise_models.noise_model import Noise + from graphix.parameter import ExpressionOrSupportsFloat, Parameter + from graphix.sim.data import Data + + +class DensityMatrix(DenseState): + """DensityMatrix object.""" + + rho: Matrix + + def __init__( + self, + data: Data = BasicStates.PLUS, + nqubit: int | None = None, + ) -> None: + """Initialize density matrix objects. + + The behaviour builds on the one of *graphix.statevec.Statevec*. + `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 iterable of scalars (A *2**n x 2**n* numerical density matrix) + - a *graphix.statevec.DensityMatrix* object + - a *graphix.statevec.Statevector* object + + 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* or *graphix.statevec.DensityMatrix* is passed, returns a copy. + + + :param data: input data to prepare the state. Can be a classical description or a numerical input, defaults to graphix.states.BasicStates.PLUS + :type data: Data + :param nqubit: number of qubits to prepare, defaults to *None* + :type nqubit: int, optional + """ + if nqubit is not None and nqubit < 0: + raise ValueError("nqubit must be a non-negative integer.") + + def check_size_consistency(mat: Matrix) -> None: + if nqubit is not None and mat.shape != (2**nqubit, 2**nqubit): + raise ValueError( + f"Inconsistent parameters between nqubit = {nqubit} and the shape of the provided density matrix = {mat.shape}." + ) + + if isinstance(data, DensityMatrix): + check_size_consistency(data.rho) + # safe: https://numpy.org/doc/stable/reference/generated/numpy.ndarray.copy.html + self.rho = data.rho.copy() + return + if isinstance(data, Iterable): + input_list = list(data) + if len(input_list) != 0 and isinstance(input_list[0], Iterable): + + def cast_row( + item: Iterable[ExpressionOrSupportsComplex] | State | Expression | SupportsFloat | SupportsComplex, + ) -> list[ExpressionOrSupportsComplex]: + if isinstance(item, Iterable): + return list(item) + raise TypeError("Every row of a matrix should be iterable.") + + input_matrix: list[list[ExpressionOrSupportsComplex]] = [cast_row(item) for item in input_list] + self.rho = np.array(input_matrix) + if not lv.is_qubitop(self.rho): + raise ValueError("Cannot interpret the provided density matrix as a qubit operator.") + check_size_consistency(self.rho) + if self.rho.dtype != np.object_: + if not lv.is_unit_trace(self.rho): + raise ValueError("Density matrix must have unit trace.") + if not lv.is_psd(self.rho): + raise ValueError("Density matrix must be positive semi-definite.") + return + statevec = Statevec(data, nqubit) + # NOTE this works since np.outer flattens the inputs! + self.rho = outer(statevec.psi, statevec.psi.conj()) + + @property + def nqubit(self) -> int: + """Return the number of qubits.""" + # Circumvent typing bug with numpy>=2.3 + # `shape` field is typed `tuple[Any, ...]` instead of `tuple[int, ...]` + # See https://github.com/numpy/numpy/issues/29830 + nqubit: int = self.rho.shape[0].bit_length() - 1 + return nqubit + + def __str__(self) -> str: + """Return a string description.""" + return f"DensityMatrix object, with density matrix {self.rho} and shape {self.dims()}." + + @override + def add_nodes(self, nqubit: int, data: Data) -> None: + r""" + Add nodes (qubits) to the density matrix and initialize them in a specified state. + + Parameters + ---------- + nqubit : int + The number of qubits to add to the density matrix. + + data : Data, optional + 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` initializes the new nodes jointly. + - A density matrix must have shape :math:`2^n \times 2^n`, + and is used to jointly initialize the new nodes. + + Notes + ----- + Previously existing nodes remain unchanged. + """ + dm_to_add = DensityMatrix(nqubit=nqubit, data=data) + self.tensor(dm_to_add) + + @override + def evolve_single(self, op: Matrix, i: int) -> None: + """Single-qubit operation. + + Parameters + ---------- + op : np.ndarray + 2*2 matrix. + i : int + Index of qubit to apply operator. + """ + assert i >= 0 + assert i < 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)) + self.rho = rho_tensor.reshape((2**self.nqubit, 2**self.nqubit)) + + @override + def evolve(self, op: Matrix, qargs: Sequence[int]) -> None: + """Multi-qubit operation. + + Args: + op (np.array): 2^n*2^n matrix + qargs (list of ints): target qubits' indexes + """ + d = op.shape + # check it is a matrix. + if len(d) == 2: + # check it is square + if d[0] == d[1]: + pass + else: + raise ValueError(f"The provided operator has shape {op.shape} and is not a square matrix.") + else: + raise ValueError(f"The provided data has incorrect shape {op.shape}.") + + nqb_op = np.log2(len(op)) + if not np.isclose(nqb_op, int(nqb_op)): + raise ValueError("Incorrect operator dimension: not consistent with qubits.") + nqb_op = int(nqb_op) + + if nqb_op != len(qargs): + 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): + raise ValueError("Incorrect target indices.") + if len(set(qargs)) != nqb_op: + raise ValueError("A repeated target qubit index is not possible.") + + op_tensor = op.reshape((2,) * 2 * nqb_op) + + 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))), + 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)))), + ) + 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))], + ) + self.rho = rho_tensor.reshape((2**self.nqubit, 2**self.nqubit)) + + @override + def expectation_single(self, op: Matrix, loc: 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. + + 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 op.shape != (2, 2): + raise ValueError("op must be 2x2 matrix.") + + st1 = copy.copy(self) + st1.normalize() + + 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) + + # complex() needed with mypy strict mode (no-any-return) + return complex(np.trace(rho_tensor.reshape((2**nqubit, 2**nqubit)))) + + def dims(self) -> tuple[int, ...]: + """Return the dimensions of the density matrix.""" + return self.rho.shape + + def tensor(self, other: DensityMatrix) -> None: + r"""Tensor product state with other density matrix. + + Results in self :math:`\otimes` other. + + Parameters + ---------- + other : :class: `DensityMatrix` object + DensityMatrix object to be tensored with self. + """ + if not isinstance(other, DensityMatrix): + other = DensityMatrix(other) + self.rho = kron(self.rho, other.rho) + + def cnot(self, edge: tuple[int, int]) -> None: + """Apply CNOT gate to density matrix. + + Parameters + ---------- + edge : (int, int) or [int, int] + Edge to apply CNOT gate. + """ + self.evolve(CNOT_TENSOR.reshape(4, 4), edge) + + @override + def swap(self, qubits: tuple[int, int]) -> None: + """Swap qubits. + + Parameters + ---------- + qubits : (int, int) + (control, target) qubits indices. + """ + self.evolve(SWAP_TENSOR.reshape(4, 4), qubits) + + def entangle(self, edge: tuple[int, int]) -> None: + """Connect graph nodes. + + Parameters + ---------- + edge : (int, int) or [int, int] + (control, target) qubit indices. + """ + self.evolve(CZ_TENSOR.reshape(4, 4), edge) + + def normalize(self) -> None: + """Normalize density matrix.""" + # 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.rho.dtype == np.object_: + rho_o = self.rho.astype(np.object_, copy=False) + rho_o /= np.trace(rho_o) + else: + rho_c = self.rho.astype(np.complex128, copy=False) + rho_c /= np.trace(rho_c) + + @override + def remove_qubit(self, qarg: int) -> None: + """Remove a qubit.""" + self.ptrace(qarg) + self.normalize() + + def ptrace(self, qargs: Collection[int] | int) -> None: + """Partial trace. + + Parameters + ---------- + qargs : list of ints or int + Indices of qubit to trace out. + """ + n = int(np.log2(self.rho.shape[0])) + if isinstance(qargs, int): + qargs = [qargs] + assert isinstance(qargs, (list, tuple)) + qargs_num = len(qargs) + nqubit_after = n - qargs_num + assert n > 0 + assert all(qarg >= 0 and qarg < n for qarg in qargs) + + rho_res = self.rho.reshape((2,) * n * 2) + # ket, bra indices to trace out + trace_axes = list(qargs) + [n + qarg for qarg in qargs] + op: Matrix = np.eye(2**qargs_num).reshape((2,) * qargs_num * 2).astype(np.complex128) + rho_res = tensordot(op, rho_res, axes=(range(2 * qargs_num), trace_axes)) + + self.rho = rho_res.reshape((2**nqubit_after, 2**nqubit_after)) + + def fidelity(self, statevec: Statevec) -> ExpressionOrFloat: + """Calculate the fidelity against reference statevector. + + Parameters + ---------- + statevec : numpy array + statevector (flattened numpy array) to compare with + """ + result = vdot(statevec.psi, matmul(self.rho, statevec.psi)) + if isinstance(result, Expression): + return result + assert math.isclose(result.imag, 0) + return result.real + + def flatten(self) -> Matrix: + """Return flattened density matrix.""" + return self.rho.flatten() + + def apply_channel(self, channel: KrausChannel, qargs: Sequence[int]) -> None: + """Apply a channel to a density matrix. + + Parameters + ---------- + :rho: density matrix. + channel: :class:`graphix.channel.KrausChannel` object + KrausChannel to be applied to the density matrix + qargs: target qubit indices + + Returns + ------- + nothing + + Raises + ------ + ValueError + If the final density matrix is not normalized after application of the channel. + This shouldn't happen since :class:`graphix.channel.KrausChannel` objects are normalized by construction. + .... + """ + result_array = np.zeros((2**self.nqubit, 2**self.nqubit), dtype=np.complex128) + + if not isinstance(channel, KrausChannel): + raise TypeError("Can't apply a channel that is not a Channel object.") + + for k_op in channel: + dm = copy.copy(self) + dm.evolve(k_op.operator, qargs) + result_array += k_op.coef * np.conj(k_op.coef) * dm.rho + # reinitialize to input density matrix + + if not np.allclose(result_array.trace(), 1.0): + raise ValueError("The output density matrix is not normalized, check the channel definition.") + + self.rho = result_array + + @override + def apply_noise(self, qubits: Sequence[int], noise: Noise) -> None: + """Apply noise. + + Parameters + ---------- + qubits : sequence of ints. + Target qubits + noise : Noise + Noise to apply + """ + channel = noise.to_kraus_channel() + self.apply_channel(channel, qubits) + + def subs(self, variable: Parameter, substitute: ExpressionOrSupportsFloat) -> DensityMatrix: + """Return a copy of the density matrix where all occurrences of the given variable in measurement angles are substituted by the given value.""" + result = copy.copy(self) + result.rho = np.vectorize(lambda value: parameter.subs(value, variable, substitute))(self.rho) + return result + + def xreplace(self, assignment: Mapping[Parameter, ExpressionOrSupportsFloat]) -> DensityMatrix: + """Return a copy of the density matrix where all occurrences of the given keys in measurement angles are substituted by the given values in parallel.""" + result = copy.copy(self) + result.rho = np.vectorize(lambda value: parameter.xreplace(value, assignment))(self.rho) + return result + + +@dataclass(frozen=True) +class DensityMatrixBackend(DenseStateBackend[DensityMatrix]): + """MBQC simulator with density matrix method.""" + + state: DensityMatrix = dataclasses.field(init=False, default_factory=lambda: DensityMatrix(nqubit=0)) diff --git a/graphix_symbolic/statevec.py b/graphix_symbolic/statevec.py new file mode 100644 index 0000000..2ab4cf2 --- /dev/null +++ b/graphix_symbolic/statevec.py @@ -0,0 +1,449 @@ +"""MBQC state vector backend supporting symbolic computations.""" + +from __future__ import annotations + +import copy +import dataclasses +import functools +import math +from collections.abc import Iterable +from dataclasses import dataclass +from typing import TYPE_CHECKING, SupportsComplex, SupportsFloat + +import numpy as np +import numpy.typing as npt +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.states import BasicStates +from typing_extensions import override + +if TYPE_CHECKING: + from collections.abc import Mapping, Sequence + from typing import Any, Literal, TypeVar + + from graphix.parameter import ExpressionOrFloat, ExpressionOrSupportsFloat, Parameter + 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, +) + + +class Statevec(DenseState): + """Statevector object.""" + + psi: Matrix + + def __init__( + self, + data: Data = BasicStates.PLUS, + nqubit: int | None = None, + ) -> None: + """Initialize statevector objects. + + `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 + + 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. + + 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 + """ + if nqubit is not None and nqubit < 0: + raise ValueError("nqubit 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: + raise ValueError( + f"Inconsistent parameters between nqubit = {nqubit} and the inferred number of qubit = {len(data.flatten())}." + ) + self.psi = data.psi.copy() + return + + # The type + # list[states.State] | list[ExpressionOrSupportsComplex] | list[Iterable[ExpressionOrSupportsComplex]] + # would be more precise, but given a value X of type Iterable[A] | Iterable[B], + # mypy infers that list(X) has type list[A | B] instead of list[A] | list[B]. + input_list: list[states.State | ExpressionOrSupportsComplex | Iterable[ExpressionOrSupportsComplex]] + if isinstance(data, states.State): + if nqubit is None: + nqubit = 1 + input_list = [data] * nqubit + elif isinstance(data, Iterable): + input_list = list(data) + else: + 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) + + elif isinstance(input_list[0], states.State): + if nqubit is None: + nqubit = len(input_list) + elif nqubit != len(input_list): + raise ValueError("Mismatch between nqubit and length of input state.") + + def state_to_statevector( + s: states.State | ExpressionOrSupportsComplex | Iterable[ExpressionOrSupportsComplex], + ) -> npt.NDArray[np.complex128]: + if not isinstance(s, states.State): + 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] + + 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)): + 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) + else: + raise TypeError(f"First element of data has type {type(input_list[0])} whereas Number or State is expected") + + def __str__(self) -> str: + """Return a string description.""" + return f"Statevec object with statevector {self.psi} and length {self.dims()}." + + @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. + + Parameters + ---------- + nqubit : int + The number of qubits to add to the state vector. + + data : Data, optional + 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. + + Notes + ----- + Previously existing nodes remain unchanged. + """ + sv_to_add = Statevec(nqubit=nqubit, data=data) + self.tensor(sv_to_add) + + @override + def evolve_single(self, op: Matrix, qubit: int) -> None: + """Apply a single-qubit operator. + + Parameters + ---------- + op : Matrix + Matrix of shape :math:`(2, 2)` representing + the operator to apply. + qubit : int + Target qubit index. + qubit index + """ + psi = tensordot(op, self.psi, (1, qubit)) + self.psi = np.moveaxis(psi, 0, qubit) + + @override + def evolve(self, op: Matrix, qubits: Sequence[int]) -> None: + """Apply a multi-qubit operator. + + Parameters + ---------- + op : Matrix + 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. + """ + 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(qubits))), qubits), + ) + self.psi = np.moveaxis(psi, range(len(qubits)), qubits) + + def dims(self) -> tuple[int, ...]: + """Return the dimensions.""" + return self.psi.shape + + # Note that `@property` must appear before `@override` for pyright + @property + @override + def nqubit(self) -> int: + """Return the number of qubits.""" + return self.psi.ndim + + @override + def remove_qubit(self, qubit: int) -> None: + r"""Remove a separable qubit from the system and assemble a statevector for remaining qubits. + + This results in the same result as partial trace, if the qubit *qubit* is separable from the rest. + + For a statevector :math:`\ket{\psi} = \sum c_i \ket{i}` with sum taken over + :math:`i \in [ 0 \dots 00,\ 0\dots 01,\ \dots,\ + 1 \dots 11 ]`, this method returns + + .. math:: + \begin{align} + \ket{\psi}' =& + c_{0 \dots 0_{\mathrm{k-1}}0_{\mathrm{k}}0_{\mathrm{k+1}} \dots 00} + \ket{0 \dots 0_{\mathrm{k-1}}0_{\mathrm{k+1}} \dots 00} \\ + & + c_{0 \dots 0_{\mathrm{k-1}}0_{\mathrm{k}}0_{\mathrm{k+1}} \dots 01} + \ket{0 \dots 0_{\mathrm{k-1}}0_{\mathrm{k+1}} \dots 01} \\ + & + c_{0 \dots 0_{\mathrm{k-1}}0_{\mathrm{k}}0_{\mathrm{k+1}} \dots 10} + \ket{0 \dots 0_{\mathrm{k-1}}0_{\mathrm{k+1}} \dots 10} \\ + & + \dots \\ + & + c_{1 \dots 1_{\mathrm{k-1}}0_{\mathrm{k}}1_{\mathrm{k+1}} \dots 11} + \ket{1 \dots 1_{\mathrm{k-1}}1_{\mathrm{k+1}} \dots 11}, + \end{align} + + (after normalization) for :math:`k =` qubit. 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}}` . + + .. warning:: + This method assumes the qubit with index *qubit* 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. + + Parameters + ---------- + qubit : 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[qubit] = 0 + psi = self.psi[tuple(index)] + norm = _norm(psi) + if isinstance(norm, SupportsFloat) and math.isclose(norm, 0): + index[qubit] = 1 + psi = self.psi[tuple(index)] + self.psi = psi + self.normalize() + + @override + def entangle(self, qubits: tuple[int, int]) -> None: + """Apply a CZ gate on two qubits. + + Parameters + ---------- + qubits : tuple[int, int] + (control, target) qubit indices. + """ + # contraction: 2nd index - control index, and 3rd index - target index. + psi = tensordot(CZ_TENSOR, self.psi, ((2, 3), qubits)) + # sort back axes + self.psi = np.moveaxis(psi, (0, 1), qubits) + + def tensor(self, other: Statevec) -> None: + r"""Tensor product state with other qubits. + + 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() + + total_num = len(self.dims()) + len(other.dims()) + self.psi = kron(psi_self, psi_other).reshape((2,) * total_num) + + def cnot(self, qubits: tuple[int, int]) -> None: + """Apply CNOT. + + Parameters + ---------- + qubits : tuple of int + (control, target) qubit indices + """ + # 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) + + @override + def swap(self, qubits: tuple[int, int]) -> None: + """Apply SWAP gate between two qubits. + + Parameters + ---------- + 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 + + def flatten(self) -> Matrix: + """Return flattened statevector.""" + return self.psi.flatten() + + @override + def expectation_single(self, op: Matrix, qubit: int) -> complex: + """Return the expectation value of single-qubit operator. + + Parameters + ---------- + op : numpy.ndarray + 2*2 operator + qubit : int + target qubit index + + Returns + ------- + complex : expectation value. + """ + st1 = copy.copy(self) + st1.normalize() + st2 = copy.copy(st1) + st1.evolve_single(op, qubit) + return complex(np.dot(st2.psi.flatten().conjugate(), st1.psi.flatten())) + + def expectation_value(self, op: Matrix, qubits: Sequence[int]) -> complex: + """Return the expectation value of multi-qubit operator. + + Parameters + ---------- + op : numpy.ndarray + 2^n*2^n operator + qubits : list of int + target qubit indices + + Returns + ------- + complex : expectation value + """ + st2 = copy.copy(self) + st2.normalize() + st1 = copy.copy(st2) + st1.evolve(op, qubits) + return complex(np.dot(st2.psi.flatten().conjugate(), st1.psi.flatten())) + + 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 + + +@dataclass(frozen=True) +class StatevectorBackend(DenseStateBackend[Statevec]): + """MBQC simulator with statevector method.""" + + state: Statevec = dataclasses.field(init=False, default_factory=lambda: Statevec(nqubit=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))) + + +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) + + +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 _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 diff --git a/requirements.txt b/requirements.txt index ed72554..c0b24eb 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,2 +1,3 @@ -graphix @ git+https://github.com/TeamGraphix/graphix +graphix @ git+https://github.com/matulni/graphix@sv-backend +# graphix @ git+https://github.com/TeamGraphix/graphix sympy>=1.9 diff --git a/tests/conftest.py b/tests/conftest.py index 149bfd1..534878b 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,9 +1,17 @@ +from __future__ import annotations + import pytest from numpy.random import PCG64, Generator SEED = 25 +DEPTH = 1 -@pytest.fixture() +@pytest.fixture def fx_rng() -> Generator: return Generator(PCG64(SEED)) + + +@pytest.fixture +def fx_bg() -> PCG64: + return PCG64(SEED) diff --git a/tests/test_statevec_numeric.py b/tests/test_statevec_numeric.py new file mode 100644 index 0000000..feabe79 --- /dev/null +++ b/tests/test_statevec_numeric.py @@ -0,0 +1,251 @@ +from __future__ import annotations + +import math +from typing import TYPE_CHECKING + +import numpy as np +import numpy.typing as npt +import pytest +from graphix.clifford import Clifford +from graphix.random_objects import rand_circuit +from graphix.sim.statevec import Statevec as SVGraphix +from graphix.sim.statevec import StatevectorBackend as SBGraphix +from graphix.states import BasicStates +from numpy.random import Generator + +from graphix_symbolic import Statevec, StatevectorBackend + +if TYPE_CHECKING: + from graphix.states import State + from numpy.random import PCG64 + + +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 + + +class TestStatevec: + N_JUMPS = 3 + + @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) + + @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) + + @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) + 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) + + @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_remove_qubit(self, sv: Statevec, q: int, sv_ref: Statevec) -> None: + sv.remove_qubit(q) + assert np.allclose(sv.flatten(), sv_ref.flatten()) + + +class TestStatevecGraphix: + """Tests in this class compare the result against the existing statevector simulator in Graphix. They are not self-contained.""" + + N_JUMPS = 3 + + @pytest.mark.parametrize("jumps", range(1, N_JUMPS)) + def test_entangle(self, fx_bg: PCG64, jumps: int) -> None: + rng = Generator(fx_bg.jumped(jumps)) + nqubits = 5 + sv_test = Statevec(generate_rnd_data(rng, nqubits)) + sv_ref = SVGraphix(data=sv_test.flatten()) + edge: tuple[int, int] = tuple(rng.choice(range(nqubits), size=2, replace=False)) + for sv in [sv_test, sv_ref]: + sv.entangle(edge) + + assert sv_ref.isclose(SVGraphix(data=sv_test.flatten())) + + @pytest.mark.parametrize("jumps", range(1, N_JUMPS)) + def test_swap(self, fx_bg: PCG64, jumps: int) -> None: + rng = Generator(fx_bg.jumped(jumps)) + nqubits = 5 + sv_test = Statevec(generate_rnd_data(rng, nqubits)) + sv_ref = SVGraphix(data=sv_test.flatten()) + edge: tuple[int, int] = tuple(rng.choice(range(nqubits), size=2, replace=False)) + for sv in [sv_test, sv_ref]: + sv.swap(edge) + + assert sv_ref.isclose(SVGraphix(data=sv_test.flatten())) + + def test_evolve_single(self, fx_rng: Generator) -> None: + nqubits = 5 + for clifford in Clifford: + sv_test = Statevec(generate_rnd_data(fx_rng, nqubits)) + sv_ref = SVGraphix(data=sv_test.flatten()) + qubit = int(fx_rng.integers(0, nqubits)) + for sv in [sv_test, sv_ref]: + sv.evolve_single(clifford.matrix, qubit) + assert sv_ref.isclose(SVGraphix(data=sv_test.flatten())) + + def test_expectation_single(self, fx_rng: Generator) -> None: + nqubits = 5 + for clifford in Clifford: + sv_test = Statevec(generate_rnd_data(fx_rng, nqubits)) + sv_ref = SVGraphix(data=sv_test.flatten()) + qubit = int(fx_rng.integers(0, nqubits)) + + val_test = sv_test.expectation_single(clifford.matrix, qubit) + val_ref = sv_ref.expectation_single(clifford.matrix, qubit) + + assert math.isclose(val_test.real, val_ref.real, abs_tol=1e-12) + assert math.isclose(val_test.imag, val_ref.imag, abs_tol=1e-12) + + def test_add_nodes(self, fx_rng: Generator) -> None: + + max_qubits = 5 + sv_test = Statevec(nqubit=0) + sv_ref = SVGraphix(nqubit=0) + + for _ in range(max_qubits): # Add a node at each iteration + data = generate_rnd_data(fx_rng, nqubits=1) + sv_test.add_nodes(1, data) + sv_ref.add_nodes(1, data) + + assert sv_ref.isclose(SVGraphix(data=sv_test.flatten())) + + @pytest.mark.parametrize( + "projector", [np.array([[1, 0], [0, 0]], dtype=np.complex128), np.array([[0, 0], [0, 1]], dtype=np.complex128)] + ) + def test_remove_nodes(self, fx_rng: Generator, projector: npt.NDArray[np.complex128]) -> None: + + nqubits = 5 + sv_test = Statevec(generate_rnd_data(fx_rng, nqubits)) + sv_ref = SVGraphix(data=sv_test.flatten()) + q = 0 + for _ in range(nqubits - 1): # Remove a node at each iteration + sv_test.evolve_single(projector, q) + sv_test.remove_qubit(q) + sv_ref.evolve_single(projector, q) + sv_ref.remove_qubit(q) + + assert sv_ref.isclose(SVGraphix(data=sv_test.flatten())) + + +@pytest.mark.parametrize("jumps", range(1, 6)) +def test_pattern_simulator(fx_bg: PCG64, jumps: int) -> None: + rng = Generator(fx_bg.jumped(jumps)) + + nqubits = 5 + + pattern = rand_circuit(nqubits, depth=5, rng=rng).transpile().pattern + pattern.remove_pauli_measurements() + + sv_test = pattern.simulate_pattern(backend=StatevectorBackend(), rng=rng) + sv_ref = pattern.simulate_pattern(backend=SBGraphix(), rng=rng) + + assert sv_ref.isclose(SVGraphix(data=sv_test.flatten())) diff --git a/tests/test_sympy_parameter.py b/tests/test_sympy_parameter.py index dd2204b..89cd36d 100644 --- a/tests/test_sympy_parameter.py +++ b/tests/test_sympy_parameter.py @@ -6,10 +6,11 @@ from graphix.branch_selector import RandomBranchSelector from numpy.random import Generator -from graphix_symbolic import SympyParameter +from graphix_symbolic import DensityMatrixBackend, StatevectorBackend, SympyParameter if TYPE_CHECKING: from graphix.parameter import Parameter + from graphix.sim.base_backend import DenseStateBackend def test_parameter_circuit_simulation(fx_rng: Generator) -> None: @@ -18,8 +19,10 @@ def test_parameter_circuit_simulation(fx_rng: Generator) -> None: circuit.rz(0, alpha) result_subs_then_simulate = circuit.subs(alpha, 0.5).simulate_statevector().statevec assert result_subs_then_simulate.psi.dtype == np.complex128 - result_simulate_then_subs = circuit.simulate_statevector().statevec.subs(alpha, 0.5) - assert np.allclose(result_subs_then_simulate.psi, result_simulate_then_subs.psi) + result_simulate_then_subs = circuit.simulate_statevector( + backend=StatevectorBackend(branch_selector=RandomBranchSelector(pr_calc=False), symbolic=True) + ).statevec.subs(alpha, 0.5) + assert np.allclose(result_subs_then_simulate.flatten(), result_simulate_then_subs.psi) def test_parameter_parallel_substitution(fx_rng: Generator) -> None: @@ -30,8 +33,10 @@ def test_parameter_parallel_substitution(fx_rng: Generator) -> None: circuit.rz(1, beta) mapping: dict[Parameter, float] = {alpha: 0.5, beta: 0.4} result_subs_then_simulate = circuit.xreplace(mapping).simulate_statevector().statevec - result_simulate_then_subs = circuit.simulate_statevector().statevec.xreplace(mapping) - assert np.allclose(result_subs_then_simulate.psi, result_simulate_then_subs.psi) + result_simulate_then_subs = circuit.simulate_statevector( + backend=StatevectorBackend(branch_selector=RandomBranchSelector(pr_calc=False), symbolic=True) + ).statevec.xreplace(mapping) + assert np.allclose(result_subs_then_simulate.flatten(), result_simulate_then_subs.flatten()) @pytest.mark.parametrize("backend", ["statevector", "densitymatrix"]) @@ -43,10 +48,14 @@ def test_parameter_pattern_simulation(backend, fx_rng: Generator) -> None: result_subs_then_simulate = pattern.subs(alpha, 0.5).simulate_pattern(backend, rng=fx_rng) # We cannot compute probabilities on symbolic states; we explore # one arbitrary branch. - result_simulate_then_subs = pattern.simulate_pattern( - backend, branch_selector=RandomBranchSelector(pr_calc=False), rng=fx_rng, symbolic=True - ).subs(alpha, 0.5) + symb_backend: DenseStateBackend if backend == "statevector": - assert np.allclose(result_subs_then_simulate.psi, result_simulate_then_subs.psi) + symb_backend = StatevectorBackend(branch_selector=RandomBranchSelector(pr_calc=False), symbolic=True) + elif backend == "densitymatrix": + symb_backend = DensityMatrixBackend(branch_selector=RandomBranchSelector(pr_calc=False), symbolic=True) + + result_simulate_then_subs = pattern.simulate_pattern(backend=symb_backend, rng=fx_rng).subs(alpha, 0.5) + if backend == "statevector": + assert np.allclose(result_subs_then_simulate.flatten(), result_simulate_then_subs.flatten()) elif backend == "densitymatrix": assert np.allclose(result_subs_then_simulate.rho, result_simulate_then_subs.rho)