diff --git a/.gitignore b/.gitignore index 7395f9c65..5bc66443c 100644 --- a/.gitignore +++ b/.gitignore @@ -14,4 +14,7 @@ graphix/sim/graphix.code-workspace graphix/graphix.code-workspace *~ .vscode/settings.json +graphix/sim/graphix.code-workspace +.vscode/settings.json +.pre-commit-config.yaml graphix/_version.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 506b6f280..bdcf446e9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -24,6 +24,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Pauli-flow finding algorithm (#117) - workflow for isort, codecov (#148, #147) +- Allow arbitrary states for initializing input nodes in state vector + and density matrix backends. + ### Fixed - Fix output node order sorting bug in Pauli preprocessing `measure_pauli` (#145) @@ -40,6 +43,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 empty output set. - Completely migrated to pytest, no `unittest` usage remains (#134) +- Basic states are now defined in `states.BasicStates` and no longer + in `ops.States`. + ## [0.2.11] - 2024-03-16 ### Added diff --git a/examples/MBQCvqe.py b/examples/MBQCvqe.py index 05eefdbc2..37c0a5e7b 100644 --- a/examples/MBQCvqe.py +++ b/examples/MBQCvqe.py @@ -22,9 +22,11 @@ import networkx as nx import numpy as np from scipy.optimize import minimize + from graphix import Circuit from graphix.simulator import PatternSimulator + # %% # Define the Hamiltonian for the VQE problem (Example: H = Z0Z1 + X0 + X1) def create_hamiltonian(): @@ -33,6 +35,7 @@ def create_hamiltonian(): H = np.kron(Z, Z) + np.kron(X, np.eye(2)) + np.kron(np.eye(2), X) return H + # %% # Function to build the VQE circuit def build_vqe_circuit(n_qubits, params): @@ -45,6 +48,7 @@ def build_vqe_circuit(n_qubits, params): circuit.cnot(i, i + 1) return circuit + # %% class MBQCVQE: def __init__(self, n_qubits, hamiltonian): @@ -85,6 +89,7 @@ def compute_energy(self, params): energy = tn.expectation_value(self.hamiltonian, qubit_indices=range(self.n_qubits)) return energy + # %% # Set parameters for VQE n_qubits = 2 @@ -94,18 +99,20 @@ def compute_energy(self, params): # Instantiate the MBQCVQE class mbqc_vqe = MBQCVQE(n_qubits, hamiltonian) + # %% # Define the cost function def cost_function(params): return mbqc_vqe.compute_energy(params) + # %% # Random initial parameters initial_params = np.random.rand(n_qubits * 3) # %% # Perform the optimization using COBYLA -result = minimize(cost_function, initial_params, method='COBYLA', options={'maxiter': 100}) +result = minimize(cost_function, initial_params, method="COBYLA", options={"maxiter": 100}) print(f"Optimized parameters: {result.x}") print(f"Optimized energy: {result.fun}") diff --git a/graphix/linalg_validations.py b/graphix/linalg_validations.py index f452504c8..a9d3fb61e 100644 --- a/graphix/linalg_validations.py +++ b/graphix/linalg_validations.py @@ -18,16 +18,21 @@ def check_square(matrix: np.ndarray) -> bool: return True +def truncate(s: str, max_length: int = 80, ellipsis: str = "...") -> str: + "Auxilliary function to truncate a long string for formatting error messages." + if len(s) <= max_length: + return s + return s[: max_length - len(ellipsis)] + ellipsis + + def check_psd(matrix: np.ndarray, tol: float = 1e-15) -> bool: """ check if a density matrix is positive semidefinite by diagonalizing. - After check_square and check_hermitian (osef) so that it already is square with power of 2 dimension. - Parameters ---------- matrix : np.ndarray - matrix to check. Normally already square and 2**n x 2**n + matrix to check tol : float tolerance on the small negatives. Default 1e-15. """ @@ -35,7 +40,7 @@ def check_psd(matrix: np.ndarray, tol: float = 1e-15) -> bool: evals = np.linalg.eigvalsh(matrix) if not all(evals >= -tol): - raise ValueError("The matrix is not positive semi-definite.") + raise ValueError("The matrix {truncate(str(matrix))} is not positive semi-definite.") return True @@ -70,7 +75,6 @@ def check_data_normalization(data: Union[list, tuple, np.ndarray]) -> bool: def check_data_dims(data: Union[list, tuple, np.ndarray]) -> bool: - # convert to set to remove duplicates dims = set([i["operator"].shape for i in data]) @@ -85,7 +89,6 @@ def check_data_dims(data: Union[list, tuple, np.ndarray]) -> bool: def check_data_values_type(data: Union[list, tuple, np.ndarray]) -> bool: - if not all( isinstance(i, dict) for i in data ): # ni liste ni ensemble mais iterable (lazy) pas stocké, executé au besoin diff --git a/graphix/ops.py b/graphix/ops.py index 33d6d507d..9404ebe10 100644 --- a/graphix/ops.py +++ b/graphix/ops.py @@ -8,16 +8,6 @@ import numpy as np -class States: - plus = np.array([1.0 / np.sqrt(2), 1.0 / np.sqrt(2)]) # plus - minus = np.array([1.0 / np.sqrt(2), -1.0 / np.sqrt(2)]) # minus - zero = np.array([1.0, 0.0]) # zero - one = np.array([0.0, 1.0]) # one - iplus = np.array([1.0 / np.sqrt(2), 1.0j / np.sqrt(2)]) # +1 eigenstate of Pauli Y - iminus = np.array([1.0 / np.sqrt(2), -1.0j / np.sqrt(2)]) # -1 eigenstate of Pauli Y - vec = [plus, minus, zero, one, iplus, iminus] - - class Ops: """Basic single- and two-qubits operators""" diff --git a/graphix/pattern.py b/graphix/pattern.py index fc4efc46a..293d48bd3 100644 --- a/graphix/pattern.py +++ b/graphix/pattern.py @@ -100,7 +100,7 @@ def add(self, cmd): cmd : list MBQC command. """ - assert type(cmd) == list + assert isinstance(cmd, list) assert cmd[0] in ["N", "E", "M", "X", "Z", "S", "C"] if cmd[0] == "N": if cmd[1] in self.__output_nodes: @@ -775,7 +775,7 @@ def get_layers(self): not_measured = set(self.__input_nodes) for cmd in self.__seq: if cmd[0] == "N": - if not cmd[1] in self.output_nodes: + if cmd[1] not in self.output_nodes: not_measured = not_measured | {cmd[1]} depth = 0 l_k = dict() @@ -814,6 +814,7 @@ def connected_edges(self, node, edges): connected: set of tuple set of connected edges """ + connected = set() for edge in edges: if edge[0] == node: @@ -1125,10 +1126,10 @@ def connected_nodes(self, node, prepared=None): if not ind == "end": # end -> 'node' is isolated while self.__seq[ind][0] == "E": if self.__seq[ind][1][0] == node: - if not self.__seq[ind][1][1] in prepared: + if self.__seq[ind][1][1] not in prepared: node_list.append(self.__seq[ind][1][1]) elif self.__seq[ind][1][1] == node: - if not self.__seq[ind][1][0] in prepared: + if self.__seq[ind][1][0] not in prepared: node_list.append(self.__seq[ind][1][0]) ind += 1 return node_list @@ -1217,7 +1218,7 @@ def _reorder_pattern(self, meas_commands): # add isolated nodes for cmd in self.__seq: if cmd[0] == "N": - if not cmd[1] in prepared: + if cmd[1] not in prepared: new.append(["N", cmd[1]]) for cmd in self.__seq: if cmd[0] == "E": diff --git a/graphix/random_objects.py b/graphix/random_objects.py index b583c2849..e88353f4e 100644 --- a/graphix/random_objects.py +++ b/graphix/random_objects.py @@ -1,5 +1,6 @@ +from __future__ import annotations + import numpy as np -import numpy.typing as npt import scipy.linalg from scipy.stats import unitary_group @@ -8,48 +9,48 @@ from graphix.sim.density_matrix import DensityMatrix -def rand_herm(l: int): +def rand_herm(sz: int): """ - generate random hermitian matrix of size l*l + generate random hermitian matrix of size sz*sz """ - tmp = np.random.rand(l, l) + 1j * np.random.rand(l, l) + tmp = np.random.rand(sz, sz) + 1j * np.random.rand(sz, sz) return tmp + tmp.conj().T -def rand_unit(l: int): +def rand_unit(sz: int): """ - generate haar random unitary matrix of size l*l + generate haar random unitary matrix of size sz*sz """ - if l == 1: + if sz == 1: return np.array([np.exp(1j * np.random.rand(1) * 2 * np.pi)]) else: - return unitary_group.rvs(l) + return unitary_group.rvs(sz) UNITS = np.array([1, 1j]) -def rand_dm(dim: int, rank: int = None, dm_dtype=True) -> DensityMatrix: - """Returns a "density matrix" as a DensityMatrix object ie a positive-semidefinite (hence Hermitian) matrix with unit trace - Note, not a proper DM since its dim can be something else than a power of 2. - The rank is random between 1 (pure) and dim if not specified - Thanks to Ulysse Chabaud. +def rand_dm(dim: int, rank: int | None = None, dm_dtype=True) -> DensityMatrix | np.ndarray: + """Utility to generate random density matrices (positive semi-definite matrices with unit trace). + Returns either a :class:`graphix.sim.density_matrix.DensityMatrix` or a :class:`np.ndarray` depending on the parameter `dm_dtype`. :param dim: Linear dimension of the (square) matrix :type dim: int - :param rank: If rank not specified then random between 1 and matrix dimension - If rank is one : then pure state else mixed state. Defaults to None + :param rank: Rank of the density matrix (1 = pure state). If not specified then sent to dim (maximal rank). + Defaults to None :type rank: int, optional - :param dm_dtype: If True returns a :class:`graphix.sim.density_matrix.DensityMatrix` or a numpy.ndarray if False. Defaults to True. + :param dm_dtype: If `True` returns a :class:`graphix.sim.density_matrix.DensityMatrix` object. If `False`returns a :class:`np.ndarray` :type dm_dtype: bool, optional - :return: Random density matrix as a :class:`graphix.sim.density_matrix.DensityMatrix` object or a numpy.ndarray. - :rtype: :class:`graphix.sim.density_matrix.DensityMatrix` or numpy.ndarray. - + :return: the density matrix in the specified format. + :rtype: DensityMatrix | np.ndarray + .. note:: + Thanks to Ulysse Chabaud. + .. warning:: + Note that setting `dm_dtype=False` allows to generate "density matrices" inconsistent with qubits i.e. with dimensions not being powers of 2. """ - # if not provided, use a random value. if rank is None: - rank = np.random.randint(1, dim + 1) + rank = dim evals = np.random.rand(rank) @@ -68,7 +69,7 @@ def rand_dm(dim: int, rank: int = None, dm_dtype=True) -> DensityMatrix: return dm -def rand_gauss_cpx_mat(dim: int, sig: float = 1 / np.sqrt(2)) -> npt.NDArray: +def rand_gauss_cpx_mat(dim: int, sig: float = 1 / np.sqrt(2)) -> np.ndarray: """ Returns a square array of standard normal complex random variates. Code from QuTiP: https://qutip.org/docs/4.0.2/modules/qutip/random_objects.html diff --git a/graphix/sim/density_matrix.py b/graphix/sim/density_matrix.py index e139c88ec..fa4d21ac2 100644 --- a/graphix/sim/density_matrix.py +++ b/graphix/sim/density_matrix.py @@ -3,57 +3,91 @@ Simulate MBQC with density matrix representation. """ +from __future__ import annotations + +import collections +import numbers +import typing from copy import deepcopy import numpy as np import graphix.sim.base_backend +import graphix.states +import graphix.types from graphix.channels import KrausChannel from graphix.clifford import CLIFFORD -from graphix.linalg_validations import check_hermitian, check_square, check_unit_trace +from graphix.linalg_validations import check_psd, check_square, check_unit_trace from graphix.ops import Ops -from graphix.sim.statevec import CNOT_TENSOR, CZ_TENSOR, SWAP_TENSOR, meas_op +from graphix.sim.statevec import CNOT_TENSOR, CZ_TENSOR, SWAP_TENSOR, Statevec class DensityMatrix: """DensityMatrix object.""" - def __init__(self, data=None, plus_state=True, nqubit=1): - """ - Parameters - ---------- - data : DensityMatrix, list, tuple, np.ndarray or None - Density matrix of shape (2**nqubits, 2**nqubits). - nqubit : int - Number of qubits. Default is 1. If both `data` and `nqubit` are specified, `nqubit` is ignored. + def __init__( + self, + data: Data = graphix.states.BasicStates.PLUS, + nqubit: graphix.types.PositiveOrNullInt | None = None, + ): + """Initialize density matrix objects. The behaviour builds on theo ne 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: graphix.states.State | "DensityMatrix" | Statevec | collections.abc.Iterable[graphix.states.State] |collections.abc.Iterable[numbers.Number] | collections.abc.Iterable[collections.abc.Iterable[numbers.Number]], optional + :param nqubit: number of qubits to prepare, defaults to None + :type nqubit: int, optional """ - if data is None: - assert nqubit >= 0 - self.Nqubit = nqubit - if plus_state: - self.rho = np.ones((2**nqubit, 2**nqubit)) / 2**nqubit - else: - self.rho = np.zeros((2**nqubit, 2**nqubit)) - self.rho[0, 0] = 1.0 - else: - if isinstance(data, DensityMatrix): - data = data.rho - elif isinstance(data, (list, tuple)): - data = np.asarray(data, dtype=complex) - elif isinstance(data, np.ndarray): - pass - else: - raise TypeError("data must be DensityMatrix, list, tuple, or np.ndarray.") + assert nqubit is None or isinstance(nqubit, numbers.Integral) and nqubit >= 0 - assert check_square(data) - self.Nqubit = len(data).bit_length() - 1 + def check_size_consistency(mat): + 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}." + ) - self.rho = data - assert check_hermitian(self.rho) - assert check_unit_trace(self.rho) + if isinstance(data, DensityMatrix): + check_size_consistency(data) + # safe: https://numpy.org/doc/stable/reference/generated/numpy.ndarray.copy.html + self.rho = data.rho.copy() + self.Nqubit = data.Nqubit + return + if isinstance(data, collections.abc.Iterable): + input_list = list(data) + if len(input_list) != 0: + # needed since Object is iterable but not subscribable! + try: + if isinstance(input_list[0], collections.abc.Iterable) and isinstance( + input_list[0][0], numbers.Number + ): + self.rho = np.array(input_list) + assert check_square(self.rho) + check_size_consistency(self.rho) + assert check_unit_trace(self.rho) + assert check_psd(self.rho) + self.Nqubit = self.rho.shape[0].bit_length() - 1 + return + except TypeError: + pass + statevec = Statevec(data, nqubit) + # NOTE this works since np.outer flattens the inputs! + self.rho = np.outer(statevec.psi, statevec.psi.conj()) + self.Nqubit = len(statevec.dims()) def __repr__(self): - return f"DensityMatrix, data={self.rho}, shape={self.dims()}" + return f"DensityMatrix object, with density matrix {self.rho} and shape {self.dims()}." def evolve_single(self, op, i): """Single-qubit operation. @@ -283,7 +317,7 @@ def apply_channel(self, channel: KrausChannel, qargs): class DensityMatrixBackend(graphix.sim.base_backend.Backend): """MBQC simulator with density matrix method.""" - def __init__(self, pattern, max_qubit_num=12, pr_calc=True): + def __init__(self, pattern, max_qubit_num=12, pr_calc=True, input_state: Data = graphix.states.BasicStates.PLUS): """ Parameters ---------- @@ -295,10 +329,11 @@ def __init__(self, pattern, max_qubit_num=12, pr_calc=True): pr_calc : bool whether or not to compute the probability distribution before choosing the measurement result. if False, measurements yield results 0/1 with 50% probabilities each. + input_state: same syntax as `graphix.statevec.DensityMatrix` constructor. """ - # check that pattern has output nodes configured - # assert len(pattern.output_nodes) > 0 self.pattern = pattern + if pattern._pauli_preprocessed and input_state != graphix.states.BasicStates.PLUS: + raise ValueError("Pauli preprocessing is currently only available when inputs are initialized in |+> state") self.results = deepcopy(pattern.results) self.state = None self.node_index = [] @@ -308,7 +343,10 @@ def __init__(self, pattern, max_qubit_num=12, pr_calc=True): raise ValueError("Pattern.max_space is larger than max_qubit_num. Increase max_qubit_num and try again.") super().__init__(pr_calc) - def add_nodes(self, nodes, qubit_to_add=None): + # initialize input qubits to desired init_state + self.add_nodes(pattern.input_nodes, input_state) + + def add_nodes(self, nodes, input_state: Data = graphix.states.BasicStates.PLUS): """add new qubit to the internal density matrix and asign the corresponding node number to list self.node_index. @@ -322,12 +360,7 @@ def add_nodes(self, nodes, qubit_to_add=None): if not self.state: self.state = DensityMatrix(nqubit=0) n = len(nodes) - if qubit_to_add is None: - dm_to_add = DensityMatrix(nqubit=n) - else: - assert isinstance(qubit_to_add, DensityMatrix) - assert qubit_to_add.nqubit == 1 - dm_to_add = qubit_to_add + dm_to_add = DensityMatrix(nqubit=n, data=input_state) self.state.tensor(dm_to_add) self.node_index.extend(nodes) self.Nqubit += n @@ -404,3 +437,24 @@ def finalize(self): """To be run at the end of pattern simulation.""" self.sort_qubits() self.state.normalize() + + +## Python <3.10: +## TypeError: unsupported operand type(s) for |: 'ABCMeta' and 'type' +## TypeError: 'ABCMeta' object is not subscriptable +# Data = ( +# graphix.states.State +# | DensityMatrix +# | Statevec +# | collections.abc.Iterable[graphix.states.State] +# | collections.abc.Iterable[numbers.Number] +# | collections.abc.Iterable[collections.abc.Iterable[numbers.Number]] +# ) +Data = typing.Union[ + graphix.states.State, + DensityMatrix, + Statevec, + typing.Iterable[graphix.states.State], + typing.Iterable[numbers.Number], + typing.Iterable[typing.Iterable[numbers.Number]], +] diff --git a/graphix/sim/statevec.py b/graphix/sim/statevec.py index 8ebbc5f08..755565cc1 100644 --- a/graphix/sim/statevec.py +++ b/graphix/sim/statevec.py @@ -1,23 +1,37 @@ +from __future__ import annotations + +import collections +import functools +import numbers +import typing from copy import deepcopy import numpy as np +import graphix.pauli import graphix.sim.base_backend -from graphix.clifford import CLIFFORD, CLIFFORD_CONJ, CLIFFORD_MUL +import graphix.states +import graphix.types +from graphix.clifford import CLIFFORD, CLIFFORD_CONJ from graphix.ops import Ops class StatevectorBackend(graphix.sim.base_backend.Backend): """MBQC simulator with statevector method.""" - def __init__(self, pattern, max_qubit_num=20, pr_calc=True): + def __init__( + self, + pattern, + input_state: Data = graphix.states.BasicStates.PLUS, + max_qubit_num=20, + pr_calc=True, + ): """ Parameters ----------- pattern : :class:`graphix.pattern.Pattern` object MBQC pattern to be simulated. - backend : str, 'statevector' - optional argument for simulation. + input_state: same syntax as `graphix.statevec.Statevec` constructor. max_qubit_num : int optional argument specifying the maximum number of qubits to be stored in the statevector at a time. @@ -28,6 +42,10 @@ def __init__(self, pattern, max_qubit_num=20, pr_calc=True): # check that pattern has output nodes configured # assert len(pattern.output_nodes) > 0 self.pattern = pattern + if pattern._pauli_preprocessed and input_state != graphix.states.BasicStates.PLUS: + raise NotImplementedError( + "Pauli preprocessing is currently only available when inputs are initialized in |+> state (see https://github.com/TeamGraphix/graphix/issues/168 )." + ) self.results = deepcopy(pattern.results) self.state = None self.node_index = [] @@ -36,9 +54,12 @@ def __init__(self, pattern, max_qubit_num=20, pr_calc=True): self.to_trace_loc = [] self.max_qubit_num = max_qubit_num if pattern.max_space() > max_qubit_num: - raise ValueError("Pattern.max_space is larger than max_qubit_num. Increase max_qubit_num and try again") + raise ValueError("Pattern.max_space is larger than max_qubit_num. Increase max_qubit_num and try again.") super().__init__(pr_calc) + # initialize input qubits to desired init_state + self.add_nodes(pattern.input_nodes, input_state) + def qubit_dim(self): """Returns the qubit number in the internal statevector @@ -48,7 +69,7 @@ def qubit_dim(self): """ return len(self.state.dims()) - def add_nodes(self, nodes): + def add_nodes(self, nodes, input_state=graphix.states.BasicStates.PLUS): """add new qubit to internal statevector and assign the corresponding node number to list self.node_index. @@ -60,7 +81,7 @@ def add_nodes(self, nodes): if not self.state: self.state = Statevec(nqubit=0) n = len(nodes) - sv_to_add = Statevec(nqubit=n) + sv_to_add = Statevec(nqubit=n, data=input_state) self.state.tensor(sv_to_add) self.node_index.extend(nodes) self.Nqubit += n @@ -180,26 +201,89 @@ def meas_op(angle, vop=0, plane="XY", choice=0): class Statevec: - """Simple statevector simulator""" + """Statevector object""" + + def __init__( + self, + data: Data = graphix.states.BasicStates.PLUS, + nqubit: graphix.types.PositiveOrNullInt | None = None, + ): + """Initialize statevector objects. The behaviour is as follows. `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. + + + :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, optional + :param nqubit: number of qubits to prepare, defaults to None + :type nqubit: int, optional + """ - def __init__(self, nqubit=1, plus_states=True): - """Initialize statevector + assert nqubit is None or isinstance(nqubit, numbers.Integral) and nqubit >= 0 - Parameters - ---------- - nqubit : int, optional: - number of qubits. Defaults to 1. - plus_states : bool, optional - whether or not to start all qubits in + state or 0 state. Defaults to + - """ - if plus_states: - self.psi = np.ones((2,) * nqubit) / 2 ** (nqubit / 2) + 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 + + if isinstance(data, graphix.states.State): + if nqubit is None: + nqubit = 1 + input_list = [data] * nqubit + elif isinstance(data, collections.abc.Iterable): + input_list = list(data) else: - self.psi = np.zeros((2,) * nqubit) - self.psi[(0,) * nqubit] = 1 + 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) + + else: + if isinstance(input_list[0], graphix.states.State): + graphix.types.check_list_elements(input_list, graphix.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.") + list_of_sv = [s.get_statevector() for s in input_list] + tmp_psi = functools.reduce(np.kron, list_of_sv) + # reshape + self.psi = tmp_psi.reshape((2,) * nqubit) + elif isinstance(input_list[0], numbers.Number): + graphix.types.check_list_elements(input_list, numbers.Number) + 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) + if 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 __repr__(self): - return f"Statevec, data={self.psi}, shape={self.dims()}" + return f"Statevec object with statevector {self.psi} and length {self.dims()}." def evolve_single(self, op, i): """Single-qubit operation @@ -258,6 +342,8 @@ def ptrace(self, qargs): rho = np.tensordot(psi, psi.conj(), axes=(qargs, qargs)) # density matrix rho = np.reshape(rho, (2**nqubit_after, 2**nqubit_after)) evals, evecs = np.linalg.eig(rho) # back to statevector + # NOTE works since only one 1 in the eigenvalues corresponding to the state + # TODO use np.eigh since rho is Hermitian? self.psi = np.reshape(evecs[:, np.argmax(evals)], (2,) * nqubit_after) def remove_qubit(self, qarg): @@ -330,6 +416,7 @@ def tensor(self, other): """ psi_self = self.psi.flatten() psi_other = other.psi.flatten() + total_num = len(self.dims()) + len(other.dims()) self.psi = np.kron(psi_self, psi_other).reshape((2,) * total_num) @@ -412,3 +499,20 @@ def expectation_value(self, op, qargs): def _get_statevec_norm(psi): """returns norm of the state""" return np.sqrt(np.sum(psi.flatten().conj() * psi.flatten())) + + +## Python <3.10: +## TypeError: unsupported operand type(s) for |: 'ABCMeta' and 'type' +## TypeError: 'ABCMeta' object is not subscriptable +# Data = ( +# graphix.states.State +# | Statevec +# | collections.abc.Iterable[graphix.states.State] +# | collections.abc.Iterable[numbers.Number] +# ) +Data = typing.Union[ + graphix.states.State, + Statevec, + typing.Iterable[graphix.states.State], + typing.Iterable[numbers.Number], +] diff --git a/graphix/sim/tensornet.py b/graphix/sim/tensornet.py index 30d942dca..fc9928cb8 100644 --- a/graphix/sim/tensornet.py +++ b/graphix/sim/tensornet.py @@ -6,7 +6,8 @@ from quimb.tensor import Tensor, TensorNetwork from graphix.clifford import CLIFFORD, CLIFFORD_CONJ, CLIFFORD_MUL -from graphix.ops import Ops, States +from graphix.ops import Ops +from graphix.states import BasicStates class TensorNetworkBackend: @@ -15,7 +16,7 @@ class TensorNetworkBackend: Executes the measurement pattern using TN expression of graph states. """ - def __init__(self, pattern, graph_prep="auto", **kwargs): + def __init__(self, pattern, graph_prep="auto", input_state=BasicStates.PLUS, **kwargs): """ Parameters @@ -32,8 +33,13 @@ def __init__(self, pattern, graph_prep="auto", **kwargs): In this strategy, All N and E commands executed sequentially. 'auto'(default) : Automatically select a preparation strategy based on the max degree of a graph + input_state : preparation for input states (only BasicStates.PLUS is supported for tensor networks yet), **kwargs : Additional keyword args to be passed to quimb.tensor.TensorNetwork. """ + if input_state != BasicStates.PLUS: + raise NotImplementedError( + "TensorNetworkBackend currently only supports |+> input state (see https://github.com/TeamGraphix/graphix/issues/167 )." + ) self.pattern = pattern self.output_nodes = pattern.output_nodes self.results = deepcopy(pattern.results) @@ -66,6 +72,9 @@ def __init__(self, pattern, graph_prep="auto", **kwargs): self._decomposed_cz = _get_decomposed_cz() self._isolated_nodes = pattern.get_isolated_nodes() + # initialize input qubits to desired init_state + self.add_nodes(pattern.input_nodes) + def add_nodes(self, nodes): """Add nodes to the network @@ -259,17 +268,17 @@ def add_qubit(self, index, state="plus"): ind = gen_str() tag = str(index) if state == "plus": - vec = States.plus + vec = BasicStates.PLUS.get_statevector() elif state == "minus": - vec = States.minus + vec = BasicStates.MINUS.get_statevector() elif state == "zero": - vec = States.zero + vec = BasicStates.ZERO.get_statevector() elif state == "one": - vec = States.one + vec = BasicStates.ONE.get_statevector() elif state == "iplus": - vec = States.iplus + vec = BasicStates.PLUS_I.get_statevector() elif state == "iminus": - vec = States.iminus + vec = BasicStates.MINUS_I.get_statevector() else: assert state.shape == (2,), "state must be 2-element np.ndarray" assert np.isclose(np.linalg.norm(state), 1), "state must be normalized" @@ -357,17 +366,17 @@ def measure_single(self, index, basis="Z", bypass_probability_calculation=True, raise Warning("Measurement outcome is chosen but the basis state was given.") proj_vec = basis elif basis == "Z" and result == 0: - proj_vec = States.zero + proj_vec = BasicStates.ZERO.get_statevector() elif basis == "Z" and result == 1: - proj_vec = States.one + proj_vec = BasicStates.ONE.get_statevector() elif basis == "X" and result == 0: - proj_vec = States.plus + proj_vec = BasicStates.PLUS.get_statevector() elif basis == "X" and result == 1: - proj_vec = States.minus + proj_vec = BasicStates.MINUS.get_statevector() elif basis == "Y" and result == 0: - proj_vec = States.iplus + proj_vec = BasicStates.PLUS_I.get_statevector() elif basis == "Y" and result == 1: - proj_vec = States.iminus + proj_vec = BasicStates.MINUS_I.get_statevector() else: raise ValueError("Invalid measurement basis.") else: @@ -414,13 +423,17 @@ def set_graph_state(self, nodes, edges): if node not in ind_dict.keys(): ind = gen_str() self._dangling[str(node)] = ind - self.add_tensor(Tensor(States.plus, [ind], [str(node), "Open"])) + self.add_tensor(Tensor(BasicStates.PLUS.get_statevector(), [ind], [str(node), "Open"])) continue dim_tensor = len(vec_dict[node]) tensor = np.array( [ - outer_product([States.vec[0 + 2 * vec_dict[node][i]] for i in range(dim_tensor)]), - outer_product([States.vec[1 + 2 * vec_dict[node][i]] for i in range(dim_tensor)]), + outer_product( + [BasicStates.VEC[0 + 2 * vec_dict[node][i]].get_statevector() for i in range(dim_tensor)] + ), + outer_product( + [BasicStates.VEC[1 + 2 * vec_dict[node][i]].get_statevector() for i in range(dim_tensor)] + ), ] ) * 2 ** (dim_tensor / 4 - 1.0 / 2) self.add_tensor(Tensor(tensor, ind_dict[node], [str(node), "Open"])) @@ -452,10 +465,10 @@ def get_basis_coefficient(self, basis, normalize=True, indices=None, **kwagrs): node = str(indices[i]) exp = len(indices) - i - 1 if (basis // 2**exp) == 1: - state_out = States.one # project onto |1> + state_out = BasicStates.ONE.get_statevector() # project onto |1> basis -= 2**exp else: - state_out = States.zero # project onto |0> + state_out = BasicStates.ZERO.get_statevector() # project onto |0> tensor = Tensor(state_out, [tn._dangling[node]], [node, f"qubit {i}", "Close"]) # retag old_ind = tn._dangling[node] @@ -707,13 +720,13 @@ def proj_basis(angle, vop, plane, choice): projected state """ if plane == "XY": - vec = States.vec[0 + choice] + vec = BasicStates.VEC[0 + choice].get_statevector() rotU = Ops.Rz(angle) elif plane == "YZ": - vec = States.vec[4 + choice] + vec = BasicStates.VEC[4 + choice].get_statevector() rotU = Ops.Rx(angle) elif plane == "XZ": - vec = States.vec[0 + choice] + vec = States.VEC[0 + choice].get_statevector() rotU = Ops.Ry(-angle) vec = np.matmul(rotU, vec) vec = np.matmul(CLIFFORD[CLIFFORD_CONJ[vop]], vec) diff --git a/graphix/simulator.py b/graphix/simulator.py index 2681d2077..65d3fdbdd 100644 --- a/graphix/simulator.py +++ b/graphix/simulator.py @@ -44,14 +44,12 @@ def __init__(self, pattern, backend="statevector", noise_model=None, **kwargs): elif backend == "densitymatrix": if noise_model is None: self.noise_model = None - # no noise: no need to compute probabilities self.backend = DensityMatrixBackend(pattern, **kwargs) warnings.warn( "Simulating using densitymatrix backend with no noise. To add noise to the simulation, give an object of `graphix.noise_models.Noisemodel` to `noise_model` keyword argument." ) if noise_model is not None: self.set_noise_model(noise_model) - # if noise: have to compute the probabilities self.backend = DensityMatrixBackend(pattern, pr_calc=True, **kwargs) elif backend in {"tensornetwork", "mps"} and noise_model is None: self.noise_model = None @@ -84,8 +82,6 @@ def run(self): the output quantum state, in the representation depending on the backend used. """ - - self.backend.add_nodes(self.pattern.input_nodes) if self.noise_model is None: for cmd in self.pattern: if cmd[0] == "N": diff --git a/graphix/states.py b/graphix/states.py new file mode 100644 index 000000000..0b0b8b232 --- /dev/null +++ b/graphix/states.py @@ -0,0 +1,75 @@ +""" +quantum states and operators +""" + +import abc + +import numpy as np +import numpy.typing as npt +import pydantic + +import graphix.pauli + + +# generic class State for all States +class State(abc.ABC): + """Abstract base class for single qubit states objects. + Only requirement for concrete classes is to have + a get_statevector() method that returns the statevector + representation of the state + """ + + @abc.abstractmethod + def get_statevector(self) -> npt.NDArray: + pass + + def get_densitymatrix(self) -> npt.NDArray: + # return DM in 2**n x 2**n dim (2x2 here) + return np.outer(self.get_statevector(), self.get_statevector().conj()) + + +class PlanarState(pydantic.BaseModel, State): + """Light object used to instantiate backends. + doesn't cover all possible states but this is + covered in :class:`graphix.sim.statevec.Statevec` + and :class:`graphix.sim.densitymatrix.DensityMatrix` + constructors. + + :param plane: One of the three planes (XY, XZ, YZ) + :type plane: :class:`graphix.pauli.Plane` + :param angle: angle IN RADIANS + :type angle: complex + :return: State + :rtype: :class:`graphix.states.State` object + """ + + plane: graphix.pauli.Plane + angle: float + + def __repr__(self) -> str: + return f"PlanarState object defined in plane {self.plane} with angle {self.angle}." + + def get_statevector(self) -> npt.NDArray: + if self.plane == graphix.pauli.Plane.XY: + return np.array([1, np.exp(1j * self.angle)]) / np.sqrt(2) + + if self.plane == graphix.pauli.Plane.YZ: + return np.array([np.cos(self.angle / 2), 1j * np.sin(self.angle / 2)]) + + if self.plane == graphix.pauli.Plane.XZ: + return np.array([np.cos(self.angle / 2), np.sin(self.angle / 2)]) + # other case never happens since exhaustive + assert False + + +# States namespace for input initialization. +class BasicStates: + ZERO = PlanarState(plane=graphix.pauli.Plane.XZ, angle=0) + ONE = PlanarState(plane=graphix.pauli.Plane.XZ, angle=np.pi) + PLUS = PlanarState(plane=graphix.pauli.Plane.XY, angle=0) + MINUS = PlanarState(plane=graphix.pauli.Plane.XY, angle=np.pi) + PLUS_I = PlanarState(plane=graphix.pauli.Plane.XY, angle=np.pi / 2) + MINUS_I = PlanarState(plane=graphix.pauli.Plane.XY, angle=-np.pi / 2) + # remove that in the end + # need in TN backend + VEC = [PLUS, MINUS, ZERO, ONE, PLUS_I, MINUS_I] diff --git a/graphix/transpiler.py b/graphix/transpiler.py index e64c6ed04..d7c31858e 100644 --- a/graphix/transpiler.py +++ b/graphix/transpiler.py @@ -8,15 +8,15 @@ import dataclasses from copy import deepcopy -from typing import Optional, Sequence +from typing import Sequence import numpy as np import graphix.pauli import graphix.sim.base_backend +import graphix.sim.statevec from graphix.ops import Ops from graphix.pattern import Pattern -from graphix.sim.statevec import Statevec @dataclasses.dataclass @@ -41,7 +41,7 @@ class SimulateResult: classical_measures : tuple[int,...], classical measures """ - statevec: Statevec + statevec: graphix.sim.statevec.Statevec classical_measures: tuple[int, ...] @@ -1358,8 +1358,8 @@ def _sort_outputs(self, pattern: Pattern, output_nodes: Sequence[int]): elif cmd[1] in old_out: cmd[1] = output_nodes[old_out.index(cmd[1])] - def simulate_statevector(self, input_state: Optional[Statevec] = None) -> SimulateResult: - """Run statevector simultion of the gate sequence, using graphix.Statevec + def simulate_statevector(self, input_state: graphix.sim.statevec.Data | None = None) -> SimulateResult: + """Run statevector simulation of the gate sequence, using graphix.Statevec Parameters ---------- @@ -1372,9 +1372,9 @@ def simulate_statevector(self, input_state: Optional[Statevec] = None) -> Simula """ if input_state is None: - state = Statevec(nqubit=self.width) + state = graphix.sim.statevec.Statevec(nqubit=self.width) else: - state = input_state + state = graphix.sim.statevec.Statevec(nqubit=self.width, data=input_state) classical_measures = [] diff --git a/graphix/types.py b/graphix/types.py new file mode 100644 index 000000000..c761768d9 --- /dev/null +++ b/graphix/types.py @@ -0,0 +1,10 @@ +import annotated_types +import typing_extensions + +PositiveOrNullInt = typing_extensions.Annotated[int, annotated_types.Ge(0)] # includes 0 + + +def check_list_elements(l, ty): + for index, item in enumerate(l): + if not isinstance(item, ty): + raise TypeError(f"data[{index}] has type {type(item)} whereas {ty} is expected") diff --git a/pyproject.toml b/pyproject.toml index 23c7c814a..b72e82b44 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -50,6 +50,9 @@ py_version = 38 line_length = 120 wrap_length = 120 +[tool.ruff] +line-length = 120 + [tool.pytest.ini_options] # Silence cotengra warning filterwarnings = ["ignore:Couldn't import `kahypar`"] diff --git a/tests/conftest.py b/tests/conftest.py index 4806bdb33..055509c1e 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,7 +1,11 @@ import pytest from numpy.random import PCG64, Generator +import graphix.transpiler +import tests.random_circuit + SEED = 25 +DEPTH = 1 @pytest.fixture() @@ -12,3 +16,25 @@ def fx_rng() -> Generator: @pytest.fixture() def fx_bg() -> PCG64: return PCG64(SEED) + + +@pytest.fixture +def hadamardpattern() -> graphix.pattern.Pattern: + circ = graphix.transpiler.Circuit(1) + circ.h(0) + return circ.transpile().pattern + + +@pytest.fixture +def nqb(fx_rng: Generator) -> int: + return fx_rng.integers(2, 5) + + +@pytest.fixture +def rand_circ(nqb, fx_rng: Generator) -> graphix.transpiler.Circuit: + return tests.random_circuit.get_rand_circuit(nqb, DEPTH, fx_rng) + + +@pytest.fixture +def randpattern(rand_circ) -> graphix.pattern.Pattern: + return rand_circ.transpile().pattern diff --git a/tests/test_density_matrix.py b/tests/test_density_matrix.py index bf9d8fdfc..1ed06a112 100644 --- a/tests/test_density_matrix.py +++ b/tests/test_density_matrix.py @@ -1,5 +1,6 @@ from __future__ import annotations +import functools import random from copy import deepcopy from typing import TYPE_CHECKING @@ -8,7 +9,9 @@ import numpy.typing as npt import pytest +import graphix.pauli import graphix.random_objects as randobj +import graphix.states from graphix import Circuit from graphix.channels import KrausChannel, dephasing_channel, depolarising_channel from graphix.ops import Ops @@ -35,9 +38,9 @@ class TestDensityMatrix: def test_init_without_data_fail(self) -> None: with pytest.raises(AssertionError): DensityMatrix(nqubit=-2) - with pytest.raises(TypeError): + with pytest.raises(AssertionError): DensityMatrix(nqubit="hello") - with pytest.raises(TypeError): + with pytest.raises(AssertionError): DensityMatrix(nqubit=[]) def test_init_with_invalid_data_fail(self, fx_rng: Generator) -> None: @@ -45,14 +48,12 @@ def test_init_with_invalid_data_fail(self, fx_rng: Generator) -> None: DensityMatrix("hello") with pytest.raises(TypeError): DensityMatrix(1) - # deprecated data shape (these test might be unnecessary) - with pytest.raises(ValueError): + with pytest.raises(TypeError): DensityMatrix([1, 2, [3]]) # check with hermitian dm but not unit trace with pytest.raises(ValueError): DensityMatrix(data=randobj.rand_herm(2 ** fx_rng.integers(2, 5))) - # check with non hermitian dm but unit trace tmp = _randdm_raw(fx_rng.integers(2, 5), fx_rng) with pytest.raises(ValueError): @@ -60,16 +61,13 @@ def test_init_with_invalid_data_fail(self, fx_rng: Generator) -> None: # check with non hermitian dm and not unit trace with pytest.raises(ValueError): DensityMatrix(data=_randdm_raw(fx_rng.integers(2, 5), fx_rng)) - # check not square matrix with pytest.raises(ValueError): # l = 2 ** fx_rng.integers(2, 5) # fx_rng.integers(2, 20) DensityMatrix(data=fx_rng.uniform(size=(3, 2))) - # check higher dimensional matrix - with pytest.raises(ValueError): - DensityMatrix(data=(2, 2, 3)) - + with pytest.raises(TypeError): + DensityMatrix(data=fx_rng.uniform(size=(2, 2, 3))) # check square and hermitian but with incorrect dimension (non-qubit type) data = randobj.rand_herm(5) data /= np.trace(data) @@ -85,7 +83,7 @@ def test_init_without_data_success(self, n: int) -> None: assert dm.rho.shape == (2**n, 2**n) assert np.allclose(dm.rho, expected_density_matrix) - dm = DensityMatrix(plus_state=False, nqubit=n) + dm = DensityMatrix(data=graphix.states.BasicStates.ZERO, nqubit=n) expected_density_matrix = np.zeros((2**n, 2**n)) expected_density_matrix[0, 0] = 1 assert dm.Nqubit == n @@ -96,13 +94,91 @@ def test_init_without_data_success(self, n: int) -> None: def test_init_with_data_success(self) -> None: # don't use rand_dm here since want to check for n in range(3): - data = randobj.rand_herm(2**n) + _data = randobj.rand_herm(2**n) + + def test_init_with_state_sucess(self, fx_rng: Generator) -> None: + # both "numerical" statevec and Statevec object + # relies on Statevec constructor validation + + nqb = fx_rng.integers(2, 5) + print(f"nqb is {nqb}") + rand_angles = fx_rng.random(nqb) * 2 * np.pi + rand_planes = fx_rng.choice([i for i in graphix.pauli.Plane], nqb) + states = [graphix.states.PlanarState(plane=i, angle=j) for i, j in zip(rand_planes, rand_angles)] + vec = Statevec(data=states) + # flattens input! + expected_dm = np.outer(vec.psi, vec.psi.conj()) + + # input with a State object + dm = DensityMatrix(data=states) + assert dm.dims() == (2**nqb, 2**nqb) + assert np.allclose(dm.rho, expected_dm) + + def test_init_with_state_fail(self, fx_rng: Generator) -> None: + nqb = 2 + rand_angles = fx_rng.random(nqb) * 2 * np.pi + rand_planes = fx_rng.choice(np.array([i for i in graphix.pauli.Plane]), nqb) + states = [graphix.states.PlanarState(plane=i, angle=j) for i, j in zip(rand_planes, rand_angles)] + + with pytest.raises(ValueError): + _dm = DensityMatrix(nqubit=1, data=states) + + with pytest.raises(ValueError): + _dm = DensityMatrix(nqubit=3, data=states) + + def test_init_with_statevec_sucess(self, fx_rng: Generator) -> None: + # both "numerical" statevec and Statevec object + # relies on Statevec constructor validation + + nqb = fx_rng.integers(2, 5) + rand_angles = fx_rng.random(nqb) * 2 * np.pi + rand_planes = fx_rng.choice(np.array([i for i in graphix.pauli.Plane]), nqb) + states = [graphix.states.PlanarState(plane=i, angle=j) for i, j in zip(rand_planes, rand_angles)] + vec = Statevec(data=states) + # flattens input! + expected_dm = np.outer(vec.psi, vec.psi.conj()) + + # input with a Statevec object + dm = DensityMatrix(data=vec) + assert dm.dims() == (2**nqb, 2**nqb) + assert np.allclose(dm.rho, expected_dm) + + sv_list = [state.get_statevector() for state in states] + sv = functools.reduce(np.kron, sv_list) - data /= np.trace(data) - dm = DensityMatrix(data=data) - assert dm.Nqubit == n - assert dm.rho.shape == (2**n, 2**n) - assert np.allclose(dm.rho, data) + # input with a statevector DATA (not Statevec object) + dm2 = DensityMatrix(data=sv) + + print("dims", dm.dims()) + assert dm2.dims() == (2**nqb, 2**nqb) + assert np.allclose(dm2.rho, dm.rho) + assert np.allclose(dm2.rho, expected_dm) + + def test_init_with_densitymatrix_sucess(self, fx_rng: Generator) -> None: + # both "numerical" densitymatrix and DensityMatrix object + + nqb = fx_rng.integers(2, 5) + rand_angles = fx_rng.random(nqb) * 2 * np.pi + rand_planes = fx_rng.choice(np.array([i for i in graphix.pauli.Plane]), nqb) + print("planes", rand_planes) + states = [graphix.states.PlanarState(plane=i, angle=j) for i, j in zip(rand_planes, rand_angles)] + vec = Statevec(data=states) + expected_dm = np.outer(vec.psi, vec.psi.conj()) + + # input with a huge density matrix + dm_list = [state.get_densitymatrix() for state in states] + num_dm = functools.reduce(np.kron, dm_list) + + dm = DensityMatrix(data=num_dm) + + assert dm.dims() == (2**nqb, 2**nqb) + assert np.allclose(dm.rho, expected_dm) + + # check copying + dm2 = DensityMatrix(dm) + assert dm2.dims() == (2**nqb, 2**nqb) + assert np.allclose(dm2.rho, expected_dm) + assert np.allclose(dm2.rho, dm.rho) def test_evolve_single_fail(self) -> None: dm = DensityMatrix(nqubit=2) @@ -441,7 +517,7 @@ def test_evolve_fail(self, fx_rng: Generator) -> None: # TODO: the test for normalization is done at initialization with data. # Now check that all operations conserve the norm. def test_normalize(self, fx_rng: Generator) -> None: - data = randobj.rand_herm(2 ** fx_rng.integers(2, 4)) + data = randobj.rand_dm(2 ** fx_rng.integers(2, 4), dm_dtype=False) dm = DensityMatrix(data / data.trace()) dm.normalize() @@ -496,9 +572,7 @@ def test_apply_dephasing_channel(self, fx_rng: Generator) -> None: # check on single qubit first # # create random density matrix # data = randobj.rand_herm(2 ** fx_rng.integers(2, 4)) - data = randobj.rand_herm(2) - data /= np.trace(data) - dm = DensityMatrix(data=data) + dm = randobj.rand_dm(2) # copy of initial dm rho_test = deepcopy(dm.rho) @@ -580,9 +654,7 @@ def test_apply_depolarising_channel(self, fx_rng: Generator) -> None: # check on single qubit first # # create random density matrix # data = randobj.rand_herm(2 ** fx_rng.integers(2, 4)) - data = randobj.rand_herm(2) - data /= np.trace(data) - dm = DensityMatrix(data=data) + dm = randobj.rand_dm(2) # copy of initial dm rho_test = deepcopy(dm.rho) @@ -793,27 +865,64 @@ def test_apply_channel_fail(self, fx_rng: Generator) -> None: class TestDensityMatrixBackend: """Test for DensityMatrixBackend class.""" - def test_init_fail(self) -> None: + # test initialization only + def test_init_success(self, fx_rng: Generator, hadamardpattern, randpattern, nqb) -> None: + # plus state (default) + backend = DensityMatrixBackend(hadamardpattern) + dm = DensityMatrix(nqubit=1) + assert np.allclose(dm.rho, backend.state.rho) + # assert backend.state.Nqubit == 1 + assert backend.state.dims() == (2, 2) + + # minus state + backend = DensityMatrixBackend(randpattern, input_state=graphix.states.BasicStates.MINUS) + dm = DensityMatrix(nqubit=nqb, data=graphix.states.BasicStates.MINUS) + assert np.allclose(dm.rho, backend.state.rho) + # assert backend.state.Nqubit == 1 + assert backend.state.dims() == (2**nqb, 2**nqb) + + rand_angles = fx_rng.random(nqb) * 2 * np.pi + rand_planes = fx_rng.choice(np.array([i for i in graphix.pauli.Plane]), nqb) + states = [graphix.states.PlanarState(plane=i, angle=j) for i, j in zip(rand_planes, rand_angles)] + + expected_dm = DensityMatrix(data=states).rho + + backend = DensityMatrixBackend(randpattern, input_state=states) + dm = backend.state + + assert dm.dims() == (2**nqb, 2**nqb) + assert np.allclose(dm.rho, expected_dm) + assert backend.Nqubit == nqb + + def test_init_fail(self, fx_rng: Generator, nqb, randpattern) -> None: + rand_angles = fx_rng.random(nqb + 1) * 2 * np.pi + rand_planes = fx_rng.choice(np.array([i for i in graphix.pauli.Plane]), nqb + 1) + states = [graphix.states.PlanarState(plane=i, angle=j) for i, j in zip(rand_planes, rand_angles)] + + # test init from State Iterable with incorrect size + with pytest.raises(ValueError): + _backend = DensityMatrixBackend(randpattern, input_state=states) + + # don't provide required pattern argument with pytest.raises(TypeError): DensityMatrixBackend() - def test_init_success(self) -> None: + def test_init_success_2(self) -> None: circ = Circuit(1) circ.rx(0, np.pi / 2) pattern = circ.transpile().pattern backend = DensityMatrixBackend(pattern) assert backend.pattern == pattern assert backend.results == pattern.results - assert backend.state is None - assert backend.node_index == [] - assert backend.Nqubit == 0 + assert backend.node_index == [0] + assert backend.Nqubit == 1 assert backend.max_qubit_num == 12 def test_add_nodes(self) -> None: circ = Circuit(1) pattern = circ.transpile().pattern backend = DensityMatrixBackend(pattern) - backend.add_nodes([0, 1]) + backend.add_nodes([1]) expected_matrix = np.array([0.25] * 16).reshape(4, 4) assert np.allclose(backend.state.rho, expected_matrix) @@ -821,7 +930,7 @@ def test_entangle_nodes(self) -> None: circ = Circuit(1) pattern = circ.transpile().pattern backend = DensityMatrixBackend(pattern) - backend.add_nodes([0, 1]) + backend.add_nodes([1]) backend.entangle_nodes((0, 1)) expected_matrix = np.array([[1, 1, 1, -1], [1, 1, 1, -1], [1, 1, 1, -1], [-1, -1, -1, 1]]) / 4 assert np.allclose(backend.state.rho, expected_matrix) @@ -835,7 +944,7 @@ def test_measure(self) -> None: pattern = circ.transpile().pattern backend = DensityMatrixBackend(pattern) - backend.add_nodes([0, 1, 2]) + backend.add_nodes([1, 2]) backend.entangle_nodes((0, 1)) backend.entangle_nodes((1, 2)) backend.measure(backend.pattern[-4]) @@ -851,7 +960,7 @@ def test_measure_pr_calc(self) -> None: pattern = circ.transpile().pattern backend = DensityMatrixBackend(pattern, pr_calc=True) - backend.add_nodes([0, 1, 2]) + backend.add_nodes([1, 2]) backend.entangle_nodes((0, 1)) backend.entangle_nodes((1, 2)) backend.measure(backend.pattern[-4]) @@ -868,7 +977,8 @@ def test_correct_byproduct(self) -> None: pattern = circ.transpile().pattern backend = DensityMatrixBackend(pattern) - backend.add_nodes([0, 1, 2]) + # node 0 initialized in Backend + backend.add_nodes([1, 2]) backend.entangle_nodes((0, 1)) backend.entangle_nodes((1, 2)) backend.measure(backend.pattern[-4]) @@ -879,7 +989,8 @@ def test_correct_byproduct(self) -> None: rho = backend.state.rho backend = StatevectorBackend(pattern) - backend.add_nodes([0, 1, 2]) + # node 0 initialized in Backend + backend.add_nodes([1, 2]) backend.entangle_nodes((0, 1)) backend.entangle_nodes((1, 2)) backend.measure(backend.pattern[-4]) diff --git a/tests/test_kraus.py b/tests/test_kraus.py index 95bccf15e..6b763d642 100644 --- a/tests/test_kraus.py +++ b/tests/test_kraus.py @@ -129,7 +129,6 @@ def test_init_with_data_fail(self, fx_rng: Generator) -> None: randobj.rand_channel_kraus(dim=2**2, rank=20) def test_dephasing_channel(self, fx_rng: Generator) -> None: - prob = fx_rng.uniform() data = [ {"coef": np.sqrt(1 - prob), "operator": np.array([[1.0, 0.0], [0.0, 1.0]])}, @@ -146,7 +145,6 @@ def test_dephasing_channel(self, fx_rng: Generator) -> None: assert np.allclose(dephase_channel.kraus_ops[i]["operator"], data[i]["operator"]) def test_depolarising_channel(self, fx_rng: Generator) -> None: - prob = fx_rng.uniform() data = [ {"coef": np.sqrt(1 - prob), "operator": np.eye(2)}, @@ -167,7 +165,6 @@ def test_depolarising_channel(self, fx_rng: Generator) -> None: assert np.allclose(depol_channel.kraus_ops[i]["operator"], data[i]["operator"]) def test_2_qubit_depolarising_channel(self, fx_rng: Generator) -> None: - prob = fx_rng.uniform() data = [ {"coef": np.sqrt(1 - prob), "operator": np.kron(np.eye(2), np.eye(2))}, @@ -200,7 +197,6 @@ def test_2_qubit_depolarising_channel(self, fx_rng: Generator) -> None: assert np.allclose(depol_channel_2_qubit.kraus_ops[i]["operator"], data[i]["operator"]) def test_2_qubit_depolarising_tensor_channel(self, fx_rng: Generator) -> None: - prob = fx_rng.uniform() data = [ {"coef": 1 - prob, "operator": np.kron(np.eye(2), np.eye(2))}, diff --git a/tests/test_pattern.py b/tests/test_pattern.py index 2308a00a1..bf44b3695 100644 --- a/tests/test_pattern.py +++ b/tests/test_pattern.py @@ -1,14 +1,18 @@ from __future__ import annotations -import itertools import sys from typing import TYPE_CHECKING, Literal import numpy as np import pytest +import graphix.ops +import graphix.sim.base_backend +import graphix.states import tests.random_circuit as rc from graphix.pattern import CommandNode, Pattern +from graphix.sim.density_matrix import DensityMatrix +from graphix.sim.statevec import Statevec from graphix.simulator import PatternSimulator from graphix.transpiler import Circuit @@ -18,6 +22,15 @@ from numpy.random import PCG64, Generator +def compare_backend_result_with_statevec(backend: str, backend_state, statevec: Statevec) -> float: + if backend == "statevector": + return np.abs(np.dot(backend_state.flatten().conjugate(), statevec.flatten())) + elif backend == "densitymatrix": + return np.abs(np.dot(backend_state.rho.flatten().conjugate(), DensityMatrix(statevec).rho.flatten())) + else: + raise NotImplementedError(backend) + + class TestPattern: # this fails without behaviour modification def test_manual_generation(self) -> None: @@ -83,7 +96,7 @@ def simulate_and_measure(): nb_shots = 1000 nb_ones = sum(1 for _ in range(nb_shots) if simulate_and_measure()) - assert abs(nb_ones - nb_shots / 2) < nb_shots / 20 + assert abs(nb_ones - nb_shots / 2) < nb_shots / 10 def test_minimize_space_graph_maxspace_with_flow(self, fx_rng: Generator) -> None: max_qubits = 20 @@ -122,7 +135,11 @@ def test_shift_signals(self, fx_bg: PCG64, jumps: int) -> None: assert np.abs(np.dot(state_mbqc.flatten().conjugate(), state.flatten())) == pytest.approx(1) @pytest.mark.parametrize("jumps", range(1, 11)) - def test_pauli_measurment(self, fx_bg: PCG64, jumps: int, use_rustworkx: bool = True) -> None: + @pytest.mark.parametrize("backend", ["statevector", "densitymatrix"]) + # TODO: tensor network backend is excluded because "parallel preparation strategy does not support not-standardized pattern". + def test_pauli_measurement_random_circuit( + self, fx_bg: PCG64, jumps: int, backend: graphix.sim.base_backend.Backend, use_rustworkx: bool = True + ) -> None: rng = Generator(fx_bg.jumped(jumps)) nqubits = 3 depth = 3 @@ -133,11 +150,13 @@ def test_pauli_measurment(self, fx_bg: PCG64, jumps: int, use_rustworkx: bool = pattern.perform_pauli_measurements(use_rustworkx=use_rustworkx) pattern.minimize_space() state = circuit.simulate_statevector().statevec - state_mbqc = pattern.simulate_pattern() - assert np.abs(np.dot(state_mbqc.flatten().conjugate(), state.flatten())) == pytest.approx(1) + state_mbqc = pattern.simulate_pattern(backend) + assert compare_backend_result_with_statevec(backend, state_mbqc, state) == pytest.approx(1) @pytest.mark.parametrize("jumps", range(1, 11)) - def test_pauli_measurment_leave_input(self, fx_bg: PCG64, jumps: int, use_rustworkx: bool = True) -> None: + def test_pauli_measurement_leave_input_random_circuit( + self, fx_bg: PCG64, jumps: int, use_rustworkx: bool = True + ) -> None: rng = Generator(fx_bg.jumped(jumps)) nqubits = 3 depth = 3 @@ -152,7 +171,7 @@ def test_pauli_measurment_leave_input(self, fx_bg: PCG64, jumps: int, use_rustwo assert np.abs(np.dot(state_mbqc.flatten().conjugate(), state.flatten())) == pytest.approx(1) @pytest.mark.parametrize("jumps", range(1, 11)) - def test_pauli_measurment_opt_gate(self, fx_bg: PCG64, jumps: int, use_rustworkx: bool = True) -> None: + def test_pauli_measurement_opt_gate(self, fx_bg: PCG64, jumps: int, use_rustworkx: bool = True) -> None: rng = Generator(fx_bg.jumped(jumps)) nqubits = 3 depth = 3 @@ -167,7 +186,7 @@ def test_pauli_measurment_opt_gate(self, fx_bg: PCG64, jumps: int, use_rustworkx assert np.abs(np.dot(state_mbqc.flatten().conjugate(), state.flatten())) == pytest.approx(1) @pytest.mark.parametrize("jumps", range(1, 11)) - def test_pauli_measurment_opt_gate_transpiler(self, fx_bg: PCG64, jumps: int, use_rustworkx: bool = True) -> None: + def test_pauli_measurement_opt_gate_transpiler(self, fx_bg: PCG64, jumps: int, use_rustworkx: bool = True) -> None: rng = Generator(fx_bg.jumped(jumps)) nqubits = 3 depth = 3 @@ -182,7 +201,7 @@ def test_pauli_measurment_opt_gate_transpiler(self, fx_bg: PCG64, jumps: int, us assert np.abs(np.dot(state_mbqc.flatten().conjugate(), state.flatten())) == pytest.approx(1) @pytest.mark.parametrize("jumps", range(1, 11)) - def test_pauli_measurment_opt_gate_transpiler_without_signalshift( + def test_pauli_measurement_opt_gate_transpiler_without_signalshift( self, fx_bg: PCG64, jumps: int, @@ -515,12 +534,24 @@ def test_pauli_measurement_end_with_measure(self) -> None: p.add(["M", 1, "XY", 0, [], []]) p.perform_pauli_measurements() + @pytest.mark.parametrize("backend", ["statevector", "densitymatrix"]) + def test_arbitrary_inputs(self, fx_rng: Generator, nqb: int, rand_circ: Circuit, backend: str) -> None: + rand_angles = fx_rng.random(nqb) * 2 * np.pi + rand_planes = fx_rng.choice(np.array([i for i in graphix.pauli.Plane]), nqb) + states = [graphix.states.PlanarState(plane=i, angle=j) for i, j in zip(rand_planes, rand_angles)] + randpattern = rand_circ.transpile().pattern + out = randpattern.simulate_pattern(backend=backend, input_state=states) + out_circ = rand_circ.simulate_statevector(input_state=states).statevec + assert compare_backend_result_with_statevec(backend, out, out_circ) == pytest.approx(1) + + def test_arbitrary_inputs_tn(self, fx_rng: Generator, nqb: int, rand_circ: Circuit) -> None: + rand_angles = fx_rng.random(nqb) * 2 * np.pi + rand_planes = fx_rng.choice(np.array([i for i in graphix.pauli.Plane]), nqb) + states = [graphix.states.PlanarState(plane=i, angle=j) for i, j in zip(rand_planes, rand_angles)] + randpattern = rand_circ.transpile().pattern + with pytest.raises(NotImplementedError): + randpattern.simulate_pattern(backend="tensornetwork", graph_prep="sequential", input_state=states) + def assert_equal_edge(edge: Sequence[int], ref: Sequence[int]) -> bool: - ans = True - for ei, ri in zip(edge, ref): - ans &= ei == ri - ansr = True - for ei, ri in zip(edge, reversed(ref)): - ansr &= ei == ri - return ans or ansr + return any(all(ei == ri for ei, ri in zip(edge, other)) for other in (ref, reversed(ref))) diff --git a/tests/test_random_utilities.py b/tests/test_random_utilities.py index ebde16524..fab85f94d 100644 --- a/tests/test_random_utilities.py +++ b/tests/test_random_utilities.py @@ -32,7 +32,6 @@ def test_rand_unit(self, fx_rng: Generator) -> None: assert np.allclose(tmp.conj().T @ tmp, np.eye(d), atol=1e-15) def test_random_channel_success(self, fx_rng: Generator) -> None: - nqb = int(fx_rng.integers(1, 5)) dim = 2**nqb # fx_rng.integers(2, 8) @@ -57,7 +56,6 @@ def test_random_channel_success(self, fx_rng: Generator) -> None: assert channel.is_normalized def test_random_channel_fail(self) -> None: - # incorrect rank type with pytest.raises(TypeError): _ = randobj.rand_channel_kraus(dim=2**2, rank=3.0) @@ -67,7 +65,6 @@ def test_random_channel_fail(self) -> None: _ = randobj.rand_channel_kraus(dim=2**2, rank=0) def test_rand_gauss_cpx(self, fx_rng: Generator) -> None: - nsample = int(1e4) dim = fx_rng.integers(2, 20) @@ -78,7 +75,6 @@ def test_rand_gauss_cpx(self, fx_rng: Generator) -> None: assert next(iter(dimset)) == (dim, dim) def test_check_psd_success(self, fx_rng: Generator) -> None: - # Generate a random mixed state from state vectors with same probability # We know this is PSD @@ -98,14 +94,13 @@ def test_check_psd_success(self, fx_rng: Generator) -> None: assert check_psd(dm) def test_check_psd_fail(self, fx_rng: Generator) -> None: - # not hermitian # don't use dim = 2, too easy to have a PSD matrix. # NOTE useless test since eigvalsh treats the matrix as hermitian and takes only the L or U part - l = fx_rng.integers(5, 20) + lst = fx_rng.integers(5, 20) - mat = fx_rng.uniform(size=(l, l)) + 1j * fx_rng.uniform(size=(l, l)) + mat = fx_rng.uniform(size=(lst, lst)) + 1j * fx_rng.uniform(size=(lst, lst)) # eigvalsh doesn't raise a LinAlgError since just use upper or lower part of the matrix. # instead Value error @@ -113,7 +108,7 @@ def test_check_psd_fail(self, fx_rng: Generator) -> None: check_psd(mat) # hermitian but not positive eigenvalues - mat = randobj.rand_herm(l) + mat = randobj.rand_herm(lst) with pytest.raises(ValueError): check_psd(mat) @@ -134,7 +129,6 @@ def test_rand_dm_fail(self, fx_rng: Generator) -> None: _ = randobj.rand_dm(2 ** fx_rng.integers(2, 5) + 1) def test_rand_dm_rank(self, fx_rng: Generator) -> None: - rk = 3 dm = randobj.rand_dm(2 ** fx_rng.integers(2, 5), rank=rk) @@ -162,7 +156,6 @@ def test_pauli_tensor_ops(self, fx_rng: Generator) -> None: assert np.all(dims == (2**nqb, 2**nqb)) def test_pauli_tensor_ops_fail(self, fx_rng: Generator) -> None: - with pytest.raises(TypeError): _ = Ops.build_tensor_Pauli_ops(fx_rng.integers(2, 6) + 0.5) @@ -170,7 +163,6 @@ def test_pauli_tensor_ops_fail(self, fx_rng: Generator) -> None: _ = Ops.build_tensor_Pauli_ops(0) def test_random_pauli_channel_success(self, fx_rng: Generator) -> None: - nqb = int(fx_rng.integers(2, 6)) rk = int(fx_rng.integers(1, 2**nqb + 1)) pauli_channel = randobj.rand_Pauli_channel_kraus(dim=2**nqb, rank=rk) # default is full rank diff --git a/tests/test_statevec.py b/tests/test_statevec.py new file mode 100644 index 000000000..92575b098 --- /dev/null +++ b/tests/test_statevec.py @@ -0,0 +1,139 @@ +import functools + +import numpy as np +import pytest + +import graphix.pauli +from graphix.sim.statevec import Statevec +from graphix.states import BasicStates, PlanarState + + +class TestStatevec: + """Test for Statevec class. Particularly new constructor.""" + + # test injitializing 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: np.random.Generator) -> None: + nqb = 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.get_statevector() for _ in range(nqb)] + sv = functools.reduce(np.kron, 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 * np.pi + rand_plane = fx_rng.choice(np.array([i for i in graphix.pauli.Plane])) + state = PlanarState(plane=rand_plane, angle=rand_angle) + vec = Statevec(nqubit=nqb, data=state) + sv_list = [state.get_statevector() for _ in range(nqb)] + sv = functools.reduce(np.kron, 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 * np.pi + rand_planes = fx_rng.choice(np.array([i for i in graphix.pauli.Plane]), nqb) + states = [PlanarState(plane=i, angle=j) for i, j in zip(rand_planes, rand_angles)] + vec = Statevec(nqubit=nqb, data=states) + sv_list = [state.get_statevector() for state in states] + sv = functools.reduce(np.kron, sv_list) + assert np.allclose(vec.psi, sv.reshape((2,) * nqb)) + assert len(vec.dims()) == nqb + + def test_data_success(self, fx_rng: np.random.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 + + # fail: incorrect len + def test_data_dim_fail(self, fx_rng: np.random.Generator) -> None: + length = 5 + 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) + + # fail: with less qubit than number of qubits inferred from a correct state vect + def test_data_dim_fail_mismatch(self, fx_rng: np.random.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)) + with pytest.raises(ValueError): + _vec = Statevec(nqubit=2, data=rand_vec) + + # fail: not normalized + def test_data_norm_fail(self, fx_rng: np.random.Generator) -> None: + nqb = fx_rng.integers(2, 5) + length = 2**nqb + rand_vec = fx_rng.random(length) + 1j * fx_rng.random(length) + with pytest.raises(ValueError): + _vec = Statevec(data=rand_vec) + + def test_defaults_to_one(self) -> None: + vec = Statevec() + assert len(vec.dims()) == 1 + + # try copying Statevec input + def test_copy_success(self, fx_rng: np.random.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) + + assert np.allclose(vec.psi, test_vec.psi) + assert len(vec.dims()) == len(test_vec.dims()) + + # try calling with incorrect number of qubits compared to inferred one + def test_copy_fail(self, fx_rng: np.random.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) + + with pytest.raises(ValueError): + _vec = Statevec(nqubit=length - 1, data=test_vec) diff --git a/tests/test_statevec_backend.py b/tests/test_statevec_backend.py index aba2efd12..f49fbc44d 100644 --- a/tests/test_statevec_backend.py +++ b/tests/test_statevec_backend.py @@ -1,13 +1,17 @@ from __future__ import annotations from copy import deepcopy +from typing import TYPE_CHECKING import numpy as np -import numpy.typing as npt import pytest -from graphix.ops import States -from graphix.sim.statevec import Statevec, meas_op +import graphix.pauli +from graphix.sim.statevec import Statevec, StatevectorBackend, meas_op +from graphix.states import BasicStates, PlanarState + +if TYPE_CHECKING: + from numpy.random import Generator class TestStatevec: @@ -28,12 +32,15 @@ def test_remove_one_qubit(self) -> None: assert np.abs(sv.psi.flatten().dot(sv2.psi.flatten().conj())) == pytest.approx(1) - @pytest.mark.parametrize("state", [States.plus, States.zero, States.one, States.iplus, States.iminus]) - def test_measurement_into_each_xyz_basis(self, state: npt.NDArray) -> None: + @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: BasicStates) -> None: n = 3 k = 0 # for measurement into |-> returns [[0, 0], ..., [0, 0]] (whose norm is zero) - m_op = np.outer(state, state.T.conjugate()) + statevector = state.get_statevector() + m_op = np.outer(statevector, statevector.T.conjugate()) sv = Statevec(nqubit=n) sv.evolve(m_op, [k]) sv.remove_qubit(k) @@ -44,8 +51,47 @@ def test_measurement_into_each_xyz_basis(self, state: npt.NDArray) -> None: def test_measurement_into_minus_state(self) -> None: n = 3 k = 0 - m_op = np.outer(States.minus, States.minus.T.conjugate()) + m_op = np.outer(BasicStates.MINUS.get_statevector(), BasicStates.MINUS.get_statevector().T.conjugate()) sv = Statevec(nqubit=n) sv.evolve(m_op, [k]) with pytest.raises(AssertionError): sv.remove_qubit(k) + + +class TestStatevecNew: + # more tests not really needed since redundant with Statevec constructor tests + + # test initialization only + def test_init_success(self, hadamardpattern, fx_rng: Generator) -> None: + # plus state (default) + backend = StatevectorBackend(hadamardpattern) + vec = Statevec(nqubit=1) + assert np.allclose(vec.psi, backend.state.psi) + assert len(backend.state.dims()) == 1 + + # minus state + backend = StatevectorBackend(hadamardpattern, input_state=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 * np.pi + rand_plane = fx_rng.choice(np.array([i for i in graphix.pauli.Plane])) + state = PlanarState(plane=rand_plane, angle=rand_angle) + backend = StatevectorBackend(hadamardpattern, input_state=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, fx_rng: Generator) -> None: + rand_angle = fx_rng.random(2) * 2 * np.pi + rand_plane = fx_rng.choice(np.array([i for i in graphix.pauli.Plane]), 2) + + state = PlanarState(plane=rand_plane[0], angle=rand_angle[0]) + state2 = PlanarState(plane=rand_plane[1], angle=rand_angle[1]) + with pytest.raises(ValueError): + StatevectorBackend(hadamardpattern, input_state=[state, state2]) diff --git a/tests/test_tnsim.py b/tests/test_tnsim.py index 5c16e7be9..d956ae3a8 100644 --- a/tests/test_tnsim.py +++ b/tests/test_tnsim.py @@ -10,8 +10,9 @@ import tests.random_circuit as rc from graphix.clifford import CLIFFORD -from graphix.ops import Ops, States +from graphix.ops import Ops from graphix.sim.tensornet import MBQCTensorNet, gen_str +from graphix.states import BasicStates from graphix.transpiler import Circuit @@ -25,7 +26,7 @@ def random_op(sites: int, dtype: type, rng: Generator) -> npt.NDArray: CZ = Ops.cz -plus = States.plus +plus = BasicStates.PLUS.get_statevector() class TestTN: diff --git a/tests/test_transpiler.py b/tests/test_transpiler.py index ff2c0a878..d6647cfe2 100644 --- a/tests/test_transpiler.py +++ b/tests/test_transpiler.py @@ -165,6 +165,6 @@ def simulate_and_measure() -> int: assert circuit_simulate.classical_measures[0] == (circuit_simulate.statevec.psi[0][1].imag > 0) return circuit_simulate.classical_measures[0] - nb_shots = 1000 + nb_shots = 10000 count = sum(1 for _ in range(nb_shots) if simulate_and_measure()) assert abs(count - nb_shots / 2) < nb_shots / 20