From 77bf8328d3e9570757aeee82defdb020f85f0ce6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Henrique=20Silv=C3=A9rio?= Date: Wed, 10 May 2023 10:39:11 +0200 Subject: [PATCH] Defining the Result class (#506) * Preliminary implementation of the Result class * First `plot_histogram()` implementation * Changing nomenclature in Results class * Update typing * Renaming `Result.sample()` to `Result.get_samples()` * Adding `Result.get_state()` * Refactoring test_simresults.py * Finish UTs --- pulser-core/pulser/result.py | 157 ++++++++++++++ .../pulser_simulation/qutip_result.py | 196 ++++++++++++++++++ .../pulser_simulation/simresults.py | 172 ++++----------- .../pulser_simulation/simulation.py | 19 +- tests/test_result.py | 129 ++++++++++++ tests/test_simresults.py | 188 +++++++++-------- tests/test_simulation.py | 2 +- 7 files changed, 644 insertions(+), 219 deletions(-) create mode 100644 pulser-core/pulser/result.py create mode 100644 pulser-simulation/pulser_simulation/qutip_result.py create mode 100644 tests/test_result.py diff --git a/pulser-core/pulser/result.py b/pulser-core/pulser/result.py new file mode 100644 index 000000000..35db67750 --- /dev/null +++ b/pulser-core/pulser/result.py @@ -0,0 +1,157 @@ +# Copyright 2023 Pulser Development Team +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +"""Classes to store measurement results.""" +from __future__ import annotations + +from abc import ABC, abstractmethod +from collections import Counter +from dataclasses import dataclass +from typing import Any + +import matplotlib.pyplot as plt +import numpy as np + +from pulser.register import QubitId + + +@dataclass +class Result(ABC): + """Base class for storing the result of a sequence run.""" + + atom_order: tuple[QubitId, ...] + meas_basis: str + + @property + def sampling_dist(self) -> dict[str, float]: + """Sampling distribution of the measured bitstring. + + Args: + atom_order: The order of the atoms in the bitstrings that + represent the measured states. + meas_basis: The measurement basis. + """ + n = self._size + return { + np.binary_repr(ind, width=n): prob + for ind, prob in enumerate(self._weights()) + if prob != 0 + } + + @property + @abstractmethod + def sampling_errors(self) -> dict[str, float]: + """The sampling error associated to each bitstring's sampling rate. + + Uses the standard error of the mean as a quantifier for sampling error. + """ + pass + + @property + def _size(self) -> int: + return len(self.atom_order) + + @abstractmethod + def _weights(self) -> np.ndarray: + """The sampling rate for every state in an ordered array.""" + pass + + def get_samples(self, n_samples: int) -> Counter[str]: + """Takes multiple samples from the sampling distribution. + + Args: + n_samples: Number of samples to return. + + Returns: + Samples of bitstrings corresponding to measured quantum states. + """ + dist = np.random.multinomial(n_samples, self._weights()) + return Counter( + { + np.binary_repr(i, self._size): dist[i] + for i in np.nonzero(dist)[0] + } + ) + + def get_state(self) -> Any: + """Gets the quantum state associated with the result. + + Can only be defined for emulation results that don't resort to + sampling a quantum state (which is the case for certain types of + noise). + """ + raise NotImplementedError( + f"`{self.__class__.__name__}.get_state()` is not implemented." + ) + + def plot_histogram( + self, + min_rate: float = 0.001, + max_n_bitstrings: int | None = None, + show: bool = True, + ) -> None: + """Plots the result in an histogram. + + Args: + min_rate: The minimum sampling rate a bitstring must have to be + displayed. + max_n_bitstrings: An optional limit on the number of bitrstrings + displayed. + show: Whether or not to call `plt.show()` before returning. + """ + # TODO: Add error bars + probs = np.array( + Counter(self.sampling_dist).most_common(max_n_bitstrings), + dtype=object, + ) + probs = probs[probs[:, 1] >= min_rate] + plt.bar(probs[:, 0], probs[:, 1]) + plt.xticks(rotation="vertical") + plt.ylabel("Probabilites") + if show: + plt.show() + + +@dataclass +class SampledResult(Result): + """Represents the result of a run from a series of samples. + + Args: + atom_order: The order of the atoms in the bitstrings that + represent the measured states. + meas_basis: The measurement basis. + bitstring_counts: The number of times each bitstring was + measured. + """ + + bitstring_counts: dict[str, int] + + def __post_init__(self) -> None: + self.n_samples = sum(self.bitstring_counts.values()) + + @property + def sampling_errors(self) -> dict[str, float]: + """The sampling error associated to each bitstring's sampling rate. + + Uses the standard error of the mean as a quantifier for sampling error. + """ + return { + bitstr: np.sqrt(p * (1 - p) / self.n_samples) + for bitstr, p in self.sampling_dist.items() + } + + def _weights(self) -> np.ndarray: + weights = np.zeros(2**self._size) + for bitstr, counts in self.bitstring_counts.items(): + weights[int(bitstr, base=2)] = counts / self.n_samples + return weights / sum(weights) diff --git a/pulser-simulation/pulser_simulation/qutip_result.py b/pulser-simulation/pulser_simulation/qutip_result.py new file mode 100644 index 000000000..30c4528b9 --- /dev/null +++ b/pulser-simulation/pulser_simulation/qutip_result.py @@ -0,0 +1,196 @@ +# Copyright 2023 Pulser Development Team +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +"""Defines a special Result subclass for simulation runs returning states.""" +from __future__ import annotations + +from dataclasses import dataclass +from typing import Union, cast + +import numpy as np +import qutip + +from pulser.result import Result + + +@dataclass +class QutipResult(Result): + """Represents the result of a run as a Qutip QObj. + + Args: + atom_order: The order of the atoms in the bitstrings that + represent the measured states. + meas_basis: The measurement basis. + state: The Qobj representing the state. Can be a statevector + or a density matrix. + matching_meas_basis: Whether the measurement basis is the + same as the state's basis. + """ + + state: qutip.Qobj + matching_meas_basis: bool + + @property + def sampling_errors(self) -> dict[str, float]: + """The sampling error associated to each bitstring's sampling rate. + + Uses the standard error of the mean as a quantifier for sampling error. + """ + return {bitstr: 0.0 for bitstr in self.sampling_dist} + + @property + def _dim(self) -> int: + full_state_size = np.prod(self.state.shape) + if not self.state.isket: + full_state_size = np.sqrt(full_state_size) + return cast( + int, np.rint(full_state_size ** (1 / self._size)).astype(int) + ) + + @property + def _basis_name(self) -> str: + if self._dim > 2: + return "all" + if self.meas_basis == "XY": + return "XY" + if not self.matching_meas_basis: + return ( + "digital" + if self.meas_basis == "ground-rydberg" + else "ground-rydberg" + ) + return self.meas_basis + + def _weights(self) -> np.ndarray: + n = self._size + if not self.state.isket: + probs = np.abs(self.state.diag()) + else: + probs = (np.abs(self.state.full()) ** 2).flatten() + + if self._dim == 2: + if self.matching_meas_basis: + # State vector ordered with r first for 'ground_rydberg' + # e.g. n=2: [rr, rg, gr, gg] -> [11, 10, 01, 00] + # Invert the order -> [00, 01, 10, 11] correspondence + # The same applies in XY mode, which is ordered with u first + weights = ( + probs if self.meas_basis == "digital" else probs[::-1] + ) + else: + # Only 000...000 is measured + weights = np.zeros(probs.size) + weights[0] = 1.0 + + elif self._dim == 3: + if self.meas_basis == "ground-rydberg": + one_state = 0 # 1 = |r> + ex_one = slice(1, 3) + elif self.meas_basis == "digital": + one_state = 2 # 1 = |h> + ex_one = slice(0, 2) + else: + raise RuntimeError( + f"Unknown measurement basis '{self.meas_basis}' " + "for a three-level system.'" + ) + probs = probs.reshape([3] * n) + weights = np.zeros(2**n) + for dec_val in range(2**n): + ind: list[Union[int, slice]] = [] + for v in np.binary_repr(dec_val, width=n): + if v == "0": + ind.append(ex_one) + else: + ind.append(one_state) + # Eg: 'digital' basis : |1> = index2, |0> = index0, 1 = 0:2 + # p_11010 = sum(probs[2, 2, 0:2, 2, 0:2]) + # We sum all probabilites that correspond to measuring + # 11010, namely hhghg, hhrhg, hhghr, hhrhr + weights[dec_val] = np.sum(probs[tuple(ind)]) + else: + raise NotImplementedError( + "Cannot sample system with single-atom state vectors of " + "dimension > 3." + ) + # Takes care of numerical artefacts in case sum(weights) != 1 + return cast(np.ndarray, weights / sum(weights)) + + def get_state( + self, + reduce_to_basis: str | None = None, + ignore_global_phase: bool = True, + tol: float = 1e-6, + normalize: bool = True, + ) -> qutip.Qobj: + """Gets the state with some optional post-processing. + + Args: + reduce_to_basis: Reduces the full state vector + to the given basis ("ground-rydberg" or "digital"), if the + population of the states to be ignored is negligible. Doesn't + apply to XY mode. + ignore_global_phase: If True and if the final state is a vector, + changes the final state's global phase such that the largest + term (in absolute value) is real. + tol: Maximum allowed population of each eliminated state. + normalize: Whether to normalize the reduced state. + + Returns: + The resulting state. + + Raises: + TypeError: If trying to reduce to a basis that would eliminate + states with significant occupation probabilites. + """ + state = self.state.copy() + is_density_matrix = state.isoper + if ignore_global_phase and not is_density_matrix: + full = state.full() + global_ph = float(np.angle(full[np.argmax(np.abs(full))])) + state *= np.exp(-1j * global_ph) + if self._dim != 3: + if reduce_to_basis not in [None, self._basis_name]: + raise TypeError( + f"Can't reduce a system in {self._basis_name}" + + f" to the {reduce_to_basis} basis." + ) + elif reduce_to_basis is not None: + if is_density_matrix: # pragma: no cover + # Not tested as noise in digital or all basis not implemented + raise NotImplementedError( + "Reduce to basis not implemented for density matrix" + " states." + ) + if reduce_to_basis == "ground-rydberg": + ex_state = "2" + elif reduce_to_basis == "digital": + ex_state = "0" + else: + raise ValueError( + "'reduce_to_basis' must be 'ground-rydberg' " + + f"or 'digital', not '{reduce_to_basis}'." + ) + ex_inds = [ + i + for i in range(3**self._size) + if ex_state in np.base_repr(i, base=3).zfill(self._size) + ] + ex_probs = np.abs(state.extract_states(ex_inds).full()) ** 2 + if not np.all(np.isclose(ex_probs, 0, atol=tol)): + raise TypeError( + "Can't reduce to chosen basis because the population of a " + "state to eliminate is above the allowed tolerance." + ) + state = state.eliminate_states(ex_inds, normalize=normalize) + return state.tidyup() diff --git a/pulser-simulation/pulser_simulation/simresults.py b/pulser-simulation/pulser_simulation/simresults.py index d35c5e9c4..7961486ee 100644 --- a/pulser-simulation/pulser_simulation/simresults.py +++ b/pulser-simulation/pulser_simulation/simresults.py @@ -16,10 +16,11 @@ from __future__ import annotations import collections.abc +import typing from abc import ABC, abstractmethod from collections import Counter from functools import lru_cache -from typing import Mapping, Optional, Tuple, Union, cast +from typing import Mapping, Optional, Tuple, TypeVar, Union, cast, overload import matplotlib.pyplot as plt import numpy as np @@ -27,8 +28,13 @@ from numpy.typing import ArrayLike from qutip.piqs import isdiagonal +from pulser.result import Result, SampledResult +from pulser_simulation.qutip_result import QutipResult -class SimulationResults(ABC): +ResultType = TypeVar("ResultType", bound=Result) + + +class SimulationResults(ABC, typing.Sequence[ResultType]): """Results of a simulation run of a pulse sequence. Parent class for NoisyResults and CoherentResults. @@ -36,6 +42,9 @@ class SimulationResults(ABC): from them. """ + # Use the pseudo-density matrix when calculating expectation values + _use_pseudo_dens: bool = False + def __init__( self, size: int, basis_name: str, sim_times: np.ndarray ) -> None: @@ -57,9 +66,7 @@ def __init__( ) self._basis_name = basis_name self._sim_times = sim_times - self._results: Union[list[Counter], list[qutip.Qobj]] - # Use the pseudo-density matrix when calculating expectation values - self._use_pseudo_dens: bool = False + self._results: list[ResultType] @property @abstractmethod @@ -77,19 +84,14 @@ def get_final_state(self) -> qutip.Qobj: """Returns the final state of the system.""" pass - @abstractmethod - def _calc_weights(self, t_index: int) -> np.ndarray: - """Computes the bitstring probabilities for sampled states.""" - pass - def expect( self, obs_list: collections.abc.Sequence[Union[qutip.Qobj, ArrayLike]] ) -> list[Union[float, complex, ArrayLike]]: """Returns the expectation values of operators in obs_list. Args: - obs_list: Input observable - list. ArrayLike objects will be converted to qutip.Qobj. + obs_list: Input observable list. ArrayLike objects will + be converted to qutip.Qobj. Returns: Expectation values of obs_list. @@ -120,8 +122,7 @@ def expect( if not isdiagonal(obs): raise ValueError(f"Observable {obs!r} is non-diagonal.") states = [ - self._calc_pseudo_density(ind) - for ind in range(len(self._results)) + self._calc_pseudo_density(ind) for ind in range(len(self)) ] else: states = self.states @@ -144,13 +145,7 @@ def sample_state( measured quantum states at time t. """ t_index = self._get_index_from_time(t, t_tol) - dist = np.random.multinomial(n_samples, self._calc_weights(t_index)) - return Counter( - { - np.binary_repr(i, self._size): dist[i] - for i in np.nonzero(dist)[0] - } - ) + return self[t_index].get_samples(n_samples) def sample_final_state(self, N_samples: int = 1000) -> Counter: """Returns the result of multiple measurements of the final state. @@ -213,7 +208,7 @@ def _proj_from_bitstring(bitstring: str) -> qutip.Qobj: ) return proj - w = self._calc_weights(t_index) + w = self[t_index]._weights() return sum( w[i] * _proj_from_bitstring(np.binary_repr(i, width=self._size)) for i in np.nonzero(w)[0] @@ -231,6 +226,24 @@ def _meas_projector(self, state_n: int) -> qutip.Qobj: # 0 = |g or d> = |1>; 1 = |r or u> = |0> return qutip.basis(2, 1 - state_n).proj() + @overload + def __getitem__(self, key: int) -> ResultType: + pass + + @overload + def __getitem__(self, key: slice) -> list[ResultType]: + pass + + def __getitem__(self, key: int | slice) -> ResultType | list[ResultType]: + return self._results[key] + + def __len__(self) -> int: + return len(self._results) + + def __iter__(self) -> collections.abc.Iterator[ResultType]: + for res in self._results: + yield res + class NoisyResults(SimulationResults): """Results of a noisy simulation run of a pulse sequence. @@ -242,9 +255,11 @@ class NoisyResults(SimulationResults): information from them. """ + _use_pseudo_dens: bool = True + def __init__( self, - run_output: list[Counter], + run_output: list[SampledResult], size: int, basis_name: str, sim_times: np.ndarray, @@ -277,8 +292,7 @@ def __init__( basis_name_ = "digital" if basis_name == "all" else basis_name super().__init__(size, basis_name_, sim_times) self.n_measures = n_measures - self._results = run_output - self._use_pseudo_dens = True + self._results: list[SampledResult] = run_output @property def states(self) -> list[qutip.Qobj]: @@ -288,7 +302,7 @@ def states(self) -> list[qutip.Qobj]: @property def results(self) -> list[Counter]: """Probability distribution of the bitstrings.""" - return self._results + return [Counter(res.sampling_dist) for res in self] def get_state(self, t: float, t_tol: float = 1.0e-3) -> qutip.Qobj: """Gets the state at time t as a diagonal density matrix. @@ -320,12 +334,6 @@ def get_final_state(self) -> qutip.Qobj: """ return self.get_state(self._sim_times[-1]) - def _calc_weights(self, t_index: int) -> np.ndarray: - weights = np.zeros(2**self._size) - for bin_rep, prob in self._results[t_index].items(): - weights[int(bin_rep, base=2)] = prob - return weights - def plot( self, op: qutip.Qobj, @@ -373,7 +381,7 @@ class CoherentResults(SimulationResults): def __init__( self, - run_output: list[qutip.Qobj], + run_output: list[QutipResult], size: int, basis_name: str, sim_times: np.ndarray, @@ -422,7 +430,7 @@ def __init__( @property def states(self) -> list[qutip.Qobj]: """List of ``qutip.Qobj`` for each state in the simulation.""" - return list(self._results) + return [res.state for res in self] def get_state( self, @@ -459,47 +467,9 @@ def get_state( states with significant occupation probabilites. """ t_index = self._get_index_from_time(t, t_tol) - state = cast(qutip.Qobj, self._results[t_index].copy()) - is_density_matrix = state.isoper - if ignore_global_phase and not is_density_matrix: - full = state.full() - global_ph = float(np.angle(full[np.argmax(np.abs(full))])) - state *= np.exp(-1j * global_ph) - if self._dim != 3: - if reduce_to_basis not in [None, self._basis_name]: - raise TypeError( - f"Can't reduce a system in {self._basis_name}" - + f" to the {reduce_to_basis} basis." - ) - elif reduce_to_basis is not None: - if is_density_matrix: # pragma: no cover - # Not tested as noise in digital or all basis not implemented - raise NotImplementedError( - "Reduce to basis not implemented for density matrix" - " states." - ) - if reduce_to_basis == "ground-rydberg": - ex_state = "2" - elif reduce_to_basis == "digital": - ex_state = "0" - else: - raise ValueError( - "'reduce_to_basis' must be 'ground-rydberg' " - + f"or 'digital', not '{reduce_to_basis}'." - ) - ex_inds = [ - i - for i in range(3**self._size) - if ex_state in np.base_repr(i, base=3).zfill(self._size) - ] - ex_probs = np.abs(state.extract_states(ex_inds).full()) ** 2 - if not np.all(np.isclose(ex_probs, 0, atol=tol)): - raise TypeError( - "Can't reduce to chosen basis because the population of a " - "state to eliminate is above the allowed tolerance." - ) - state = state.eliminate_states(ex_inds, normalize=normalize) - return state.tidyup() + return self[t_index].get_state( + reduce_to_basis, ignore_global_phase, tol, normalize + ) def get_final_state( self, @@ -538,58 +508,6 @@ def get_final_state( normalize, ) - def _calc_weights(self, t_index: int) -> np.ndarray: - n = self._size - state_t = cast(qutip.Qobj, self._results[t_index]).unit() - if state_t.type != "ket": - probs = np.abs(state_t.diag()) - else: - probs = (np.abs(state_t.full()) ** 2).flatten() - - if self._dim == 2: - if self._meas_basis == self._basis_name: - # State vector ordered with r first for 'ground_rydberg' - # e.g. n=2: [rr, rg, gr, gg] -> [11, 10, 01, 00] - # Invert the order -> [00, 01, 10, 11] correspondence - # The same applies in XY mode, which is ordered with u first - weights = ( - probs if self._meas_basis == "digital" else probs[::-1] - ) - else: - # Only 000...000 is measured - weights = np.zeros(probs.size) - weights[0] = 1.0 - - elif self._dim == 3: - if self._meas_basis == "ground-rydberg": - one_state = 0 # 1 = |r> - ex_one = slice(1, 3) - elif self._meas_basis == "digital": - one_state = 2 # 1 = |h> - ex_one = slice(0, 2) - probs = probs.reshape([3] * n) - weights = np.zeros(2**n) - for dec_val in range(2**n): - ind: list[Union[int, slice]] = [] - for v in np.binary_repr(dec_val, width=n): - if v == "0": - ind.append(ex_one) - else: - ind.append(one_state) - # Eg: 'digital' basis : |1> = index2, |0> = index0, 1 = 0:2 - # p_11010 = sum(probs[2, 2, 0:2, 2, 0:2]) - # We sum all probabilites that correspond to measuring - # 11010, namely hhghg, hhrhg, hhghr, hhrhr - weights[dec_val] = np.sum(probs[tuple(ind)]) - else: - raise NotImplementedError( - "Cannot sample system with single-atom state vectors of " - "dimension > 3." - ) - # Takes care of numerical artefacts in case sum(weights) != 1 - weights /= sum(weights) - return cast(np.ndarray, weights) - def _meas_projector(self, state_n: int) -> qutip.Qobj: if self._meas_errors: err_param = ( diff --git a/pulser-simulation/pulser_simulation/simulation.py b/pulser-simulation/pulser_simulation/simulation.py index 9659a5813..bef615b80 100644 --- a/pulser-simulation/pulser_simulation/simulation.py +++ b/pulser-simulation/pulser_simulation/simulation.py @@ -30,8 +30,10 @@ import pulser.sampler as sampler from pulser import Pulse, Sequence from pulser.register import QubitId +from pulser.result import SampledResult from pulser.sampler.samples import _TargetSlot from pulser.sequence._seq_drawer import draw_sequence +from pulser_simulation.qutip_result import QutipResult from pulser_simulation.simconfig import SimConfig from pulser_simulation.simresults import ( CoherentResults, @@ -970,8 +972,17 @@ def _run_solver() -> CoherentResults: progress_bar=p_bar, options=solv_ops, ) + results = [ + QutipResult( + tuple(self._qdict), + self._meas_basis, + state, + self._meas_basis == self.basis_name, + ) + for state in result.states + ] return CoherentResults( - result.states, + results, self._size, self.basis_name, self._eval_times_array, @@ -1036,12 +1047,12 @@ def _run_solver() -> CoherentResults: ] ) n_measures = self.config.runs * self.config.samples_per_run - total_run_prob = [ - Counter({k: v / n_measures for k, v in total_count[t].items()}) + results = [ + SampledResult(tuple(self._qdict), self._meas_basis, total_count[t]) for t in time_indices ] return NoisyResults( - total_run_prob, + results, self._size, self.basis_name, self._eval_times_array, diff --git a/tests/test_result.py b/tests/test_result.py new file mode 100644 index 000000000..9e41ebb96 --- /dev/null +++ b/tests/test_result.py @@ -0,0 +1,129 @@ +# Copyright 2023 Pulser Development Team +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +from __future__ import annotations + +import re +from collections import Counter + +import numpy as np +import pytest +import qutip + +from pulser.result import SampledResult +from pulser_simulation.qutip_result import QutipResult + + +def test_sampled_result(patch_plt_show): + samples = Counter({"000": 50, "111": 50}) + result = SampledResult( + atom_order=("a", "b", "c"), + meas_basis="ground-rydberg", + bitstring_counts=samples, + ) + assert result.n_samples == 100 + assert result.sampling_dist == {"000": 0.5, "111": 0.5} + sampling_err = np.sqrt(0.5**2 / 100) + assert result.sampling_errors == { + "000": sampling_err, + "111": sampling_err, + } + n_samples = 100 + np.random.seed(3052023) + new_samples = result.get_samples(100) + new_samples.subtract(samples) + assert all( + abs(counts_diff) < sampling_err * n_samples + for counts_diff in new_samples.values() + ) + + with pytest.raises( + NotImplementedError, + match=re.escape("`SampledResult.get_state()` is not implemented"), + ): + result.get_state() + + result.plot_histogram() + + +def test_qutip_result(): + qutrit_state = qutip.tensor(qutip.basis(3, 0), qutip.basis(3, 1)) + result = QutipResult( + atom_order=("q0", "q1"), + meas_basis="ground-rydberg", + state=qutrit_state, + matching_meas_basis=True, + ) + assert result.sampling_dist == {"10": 1.0} + assert result.sampling_errors == {"10": 0.0} + assert result._basis_name == "all" + + assert result.get_state() == qutrit_state + qubit_state = qutip.tensor(qutip.basis(2, 0), qutip.basis(2, 1)) + np.testing.assert_array_equal( + result.get_state(reduce_to_basis="ground-rydberg").full(), + qubit_state.full(), + ) + + result.meas_basis = "digital" + assert result.sampling_dist == {"00": 1.0} + + result.meas_basis = "XY" + with pytest.raises(RuntimeError, match="Unknown measurement basis 'XY'"): + result.sampling_dist + + new_result = QutipResult( + atom_order=("q0", "q1"), + meas_basis="digital", + state=qubit_state, + matching_meas_basis=True, + ) + assert new_result.sampling_dist == {"01": 1.0} + + new_result.meas_basis = "ground-rydberg" + assert new_result.sampling_dist == {"10": 1.0} + + new_result.matching_meas_basis = False + assert new_result.sampling_dist == {"00": 1.0} + # The state basis should be infered to be "digital" + with pytest.raises( + TypeError, + match="Can't reduce a system in digital to the ground-rydberg basis", + ): + new_result.get_state(reduce_to_basis="ground-rydberg") + + oversized_state = qutip.Qobj(np.eye(16) / 16) + result.state = oversized_state + assert result._dim == 4 + with pytest.raises( + NotImplementedError, + match="Cannot sample system with single-atom state vectors of" + " dimension > 3", + ): + result.sampling_dist + + density_matrix = qutip.Qobj(np.eye(4) / 4) + result = QutipResult( + atom_order=("a", "b"), + meas_basis="ground-rydberg", + state=density_matrix, + matching_meas_basis=True, + ) + assert result.state.isoper + assert result._dim == 2 + assert result.sampling_dist == { + "00": 0.25, + "01": 0.25, + "10": 0.25, + "11": 0.25, + } diff --git a/tests/test_simresults.py b/tests/test_simresults.py index 28017125c..20adcf15d 100644 --- a/tests/test_simresults.py +++ b/tests/test_simresults.py @@ -12,7 +12,7 @@ # See the License for the specific language governing permissions and # limitations under the License. from collections import Counter -from copy import deepcopy +from typing import cast import numpy as np import pytest @@ -25,79 +25,97 @@ from pulser_simulation import SimConfig, Simulation from pulser_simulation.simresults import CoherentResults, NoisyResults -np.random.seed(123) -q_dict = { - "A": np.array([0.0, 0.0]), - "B": np.array([0.0, 10.0]), -} -reg = Register(q_dict) -duration = 1000 -pi_pulse = Pulse.ConstantDetuning(BlackmanWaveform(duration, np.pi), 0.0, 0) +@pytest.fixture +def reg(): + q_dict = { + "A": np.array([0.0, 0.0]), + "B": np.array([0.0, 10.0]), + } + return Register(q_dict) -seq = Sequence(reg, Chadoq2) -# Declare Channels -seq.declare_channel("ryd", "rydberg_global") -seq.add(pi_pulse, "ryd") -seq_no_meas = deepcopy(seq) -seq_no_meas_noisy = deepcopy(seq) -seq.measure("ground-rydberg") +@pytest.fixture +def pi_pulse(): + return Pulse.ConstantDetuning(BlackmanWaveform(1000, np.pi), 0.0, 0) -cfg_noisy = SimConfig(noise=("SPAM", "doppler", "amplitude"), amp_sigma=1e-3) -cfg_noisych = SimConfig(noise="dephasing", dephasing_prob=0.01) -sim = Simulation(seq) -sim_noisy = Simulation(seq, config=cfg_noisy) -sim_noisych = Simulation(seq, config=cfg_noisych) +@pytest.fixture +def seq_no_meas(reg, pi_pulse): + seq = Sequence(reg, Chadoq2) + seq.declare_channel("ryd", "rydberg_global") + seq.add(pi_pulse, "ryd") + return seq -results = sim.run() -results_noisy = sim_noisy.run() -results_noisych = sim_noisych.run() -state = qutip.tensor([qutip.basis(2, 0), qutip.basis(2, 0)]) -ground = qutip.tensor([qutip.basis(2, 1), qutip.basis(2, 1)]) +@pytest.fixture +def sim(seq_no_meas): + seq_no_meas.measure("ground-rydberg") + np.random.seed(123) + return Simulation(seq_no_meas) + + +@pytest.fixture +def results_noisy(sim): + sim.add_config( + SimConfig(noise=("SPAM", "doppler", "amplitude"), amp_sigma=1e-3) + ) + return sim.run() -def test_initialization(): +@pytest.fixture +def results(sim): + return sim.run() + + +def test_initialization(results): + rr_state = qutip.tensor([qutip.basis(2, 0), qutip.basis(2, 0)]) with pytest.raises(ValueError, match="`basis_name` must be"): - CoherentResults(state, 2, "bad_basis", None, [0]) + CoherentResults(rr_state, 2, "bad_basis", None, [0]) with pytest.raises( ValueError, match="`meas_basis` must be 'ground-rydberg' or 'digital'." ): - CoherentResults(state, 1, "all", None, "XY") + CoherentResults(rr_state, 1, "all", None, "XY") with pytest.raises( ValueError, match="`meas_basis` and `basis_name` must have the same value.", ): CoherentResults( - state, 1, "ground-rydberg", [0], "wrong_measurement_basis" + rr_state, 1, "ground-rydberg", [0], "wrong_measurement_basis" ) with pytest.raises(ValueError, match="`basis_name` must be"): - NoisyResults(state, 2, "bad_basis", [0], 123) + NoisyResults(rr_state, 2, "bad_basis", [0], 123) with pytest.raises( ValueError, match="only values of 'epsilon' and 'epsilon_prime'" ): CoherentResults( - state, + rr_state, 1, "ground-rydberg", [0], "ground-rydberg", - cfg_noisy.spam_dict, + {"eta": 0.1, "epsilon": 0.0, "epsilon_prime": 0.4}, ) assert results._dim == 2 assert results._size == 2 assert results._basis_name == "ground-rydberg" assert results._meas_basis == "ground-rydberg" - assert results.states[0] == ground + assert results.states[0] == qutip.tensor( + [qutip.basis(2, 1), qutip.basis(2, 1)] + ) @pytest.mark.parametrize("noisychannel", [True, False]) -def test_get_final_state(noisychannel): - _results = results_noisych if noisychannel else results - assert _results.get_final_state().isoper or not noisychannel +def test_get_final_state( + noisychannel, sim: Simulation, results, reg, pi_pulse +): + if noisychannel: + sim.add_config(SimConfig(noise="dephasing", dephasing_prob=0.01)) + _results = sim.run() + assert isinstance(_results, CoherentResults) + final_state = _results.get_final_state() + assert final_state.isoper if noisychannel else final_state.isket with pytest.raises(TypeError, match="Can't reduce"): _results.get_final_state(reduce_to_basis="digital") assert ( @@ -131,6 +149,7 @@ def test_get_final_state(noisychannel): sim_ = Simulation(seq_) results_ = sim_.run() + results_ = cast(CoherentResults, results_) with pytest.raises(ValueError, match="'reduce_to_basis' must be"): results_.get_final_state(reduce_to_basis="all") @@ -155,7 +174,7 @@ def test_get_final_state(noisychannel): ) -def test_get_final_state_noisy(): +def test_get_final_state_noisy(reg, pi_pulse): np.random.seed(123) seq_ = Sequence(reg, Chadoq2) seq_.declare_channel("ram", "raman_local", initial_target="A") @@ -177,7 +196,7 @@ def test_get_final_state_noisy(): ) -def test_get_state_float_time(): +def test_get_state_float_time(results): with pytest.raises(IndexError, match="is absent from"): results.get_state(-1.0) with pytest.raises(IndexError, match="is absent from"): @@ -199,7 +218,7 @@ def test_get_state_float_time(): ).all() -def test_expect(): +def test_expect(results, pi_pulse, reg): with pytest.raises(TypeError, match="must be a list"): results.expect("bad_observable") with pytest.raises(TypeError, match="Incompatible type"): @@ -215,7 +234,7 @@ def test_expect(): op = [qutip.basis(2, 0).proj()] exp = results_single.expect(op)[0] assert np.isclose(exp[-1], 1) - assert len(exp) == duration + 1 # +1 for the final instant + assert len(exp) == pi_pulse.duration + 1 # +1 for the final instant np.testing.assert_almost_equal( results_single._calc_pseudo_density(-1).full(), np.array([[1, 0], [0, 0]]), @@ -252,59 +271,63 @@ def test_expect(): assert np.isclose(exp3dim[0][-1], 1.89690200e-14) -def test_expect_noisy(): +def test_expect_noisy(results_noisy): np.random.seed(123) bad_op = qutip.tensor([qutip.qeye(2), qutip.sigmap()]) with pytest.raises(ValueError, match="non-diagonal"): results_noisy.expect([bad_op]) op = qutip.tensor([qutip.qeye(2), qutip.basis(2, 0).proj()]) - assert np.isclose(results_noisy.expect([op])[0][-1], 0.7333333333333334) + assert np.isclose(results_noisy.expect([op])[0][-1], 0.7466666666666666) -def test_plot(): +def test_plot(results_noisy, results): op = qutip.tensor([qutip.qeye(2), qutip.basis(2, 0).proj()]) results_noisy.plot(op) results_noisy.plot(op, error_bars=False) results.plot(op) -def test_sample_final_state(): - np.random.seed(123) +def test_sim_without_measurement(seq_no_meas): + assert not seq_no_meas.is_measured() sim_no_meas = Simulation(seq_no_meas, config=SimConfig(runs=1)) results_no_meas = sim_no_meas.run() assert results_no_meas.sample_final_state() == Counter( - {"00": 88, "01": 156, "10": 188, "11": 568} + {"00": 80, "01": 164, "10": 164, "11": 592} ) - with pytest.raises(NotImplementedError, match="dimension > 3"): - results_large_dim = deepcopy(results) - results_large_dim._dim = 7 - results_large_dim.sample_final_state() + +def test_sample_final_state(results): sampling = results.sample_final_state(1234) assert len(sampling) == 4 # Check that all states were observed. - results._meas_basis = "digital" + # Switch the measurement basis in the result + results[-1].matching_meas_basis = False sampling0 = results.sample_final_state(N_samples=911) assert sampling0 == {"00": 911} + + +def test_sample_final_state_three_level(seq_no_meas, pi_pulse): seq_no_meas.declare_channel("raman", "raman_local", "B") seq_no_meas.add(pi_pulse, "raman") res_3level = Simulation(seq_no_meas).run() # Raman pi pulse on one atom will not affect other, # even with global pi on rydberg assert len(res_3level.sample_final_state()) == 2 - res_3level._meas_basis = "ground-rydberg" - sampling_three_levelB = res_3level.sample_final_state() + + seq_no_meas.measure("ground-rydberg") + res_3level_gb = Simulation(seq_no_meas).run() + sampling_three_levelB = res_3level_gb.sample_final_state() # Rydberg will affect both: assert len(sampling_three_levelB) == 4 -def test_sample_final_state_noisy(): +def test_sample_final_state_noisy(seq_no_meas, results_noisy): np.random.seed(123) assert results_noisy.sample_final_state(N_samples=1234) == Counter( - {"11": 725, "10": 265, "01": 192, "00": 52} + {"11": 772, "10": 190, "01": 161, "00": 111} ) res_3level = Simulation( - seq_no_meas_noisy, config=SimConfig(noise=("SPAM", "doppler"), runs=10) + seq_no_meas, config=SimConfig(noise=("SPAM", "doppler"), runs=10) ) final_state = res_3level.run().states[-1] assert np.isclose( @@ -320,51 +343,42 @@ def test_sample_final_state_noisy(): ).all() -def test_results_xy(): - q_dict = { - "A": np.array([0.0, 0.0]), - "B": np.array([0.0, 10.0]), - } - reg = Register(q_dict) - duration = 1000 - pi_pulse = Pulse.ConstantDetuning( - BlackmanWaveform(duration, np.pi), 0.0, 0 - ) - seq = Sequence(reg, MockDevice) +def test_results_xy(reg, pi_pulse): + seq_ = Sequence(reg, MockDevice) # Declare Channels - seq.declare_channel("ch0", "mw_global") - seq.add(pi_pulse, "ch0") - seq.measure("XY") - - sim = Simulation(seq) - results = sim.run() + seq_.declare_channel("ch0", "mw_global") + seq_.add(pi_pulse, "ch0") + seq_.measure("XY") - ground = qutip.tensor([qutip.basis(2, 1), qutip.basis(2, 1)]) + sim_ = Simulation(seq_) + results_ = sim_.run() - assert results._dim == 2 - assert results._size == 2 - assert results._basis_name == "XY" - assert results._meas_basis == "XY" - assert results.states[0] == ground + assert results_._dim == 2 + assert results_._size == 2 + assert results_._basis_name == "XY" + assert results_._meas_basis == "XY" + assert results_.states[0] == qutip.tensor( + [qutip.basis(2, 1), qutip.basis(2, 1)] + ) with pytest.raises(TypeError, match="Can't reduce a system in"): - results.get_final_state(reduce_to_basis="all") + results_.get_final_state(reduce_to_basis="all") with pytest.raises(TypeError, match="Can't reduce a system in"): - results.get_final_state(reduce_to_basis="ground-rydberg") + results_.get_final_state(reduce_to_basis="ground-rydberg") with pytest.raises(TypeError, match="Can't reduce a system in"): - results.get_final_state(reduce_to_basis="digital") + results_.get_final_state(reduce_to_basis="digital") - state = results.get_final_state(reduce_to_basis="XY") + state = results_.get_final_state(reduce_to_basis="XY") assert np.all( np.isclose( - np.abs(state.full()), np.abs(results.states[-1].full()), atol=1e-5 + np.abs(state.full()), np.abs(results_.states[-1].full()), atol=1e-5 ) ) # Check that measurement projectors are correct - assert results._meas_projector(0) == qutip.basis(2, 1).proj() - assert results._meas_projector(1) == qutip.basis(2, 0).proj() + assert results_._meas_projector(0) == qutip.basis(2, 1).proj() + assert results_._meas_projector(1) == qutip.basis(2, 0).proj() diff --git a/tests/test_simulation.py b/tests/test_simulation.py index 8c54479c2..616feb754 100644 --- a/tests/test_simulation.py +++ b/tests/test_simulation.py @@ -643,7 +643,7 @@ def test_depolarizing(): seq, sampling_rate=0.01, config=SimConfig(noise="depolarizing") ) assert sim.run().sample_final_state() == Counter({"0": 587, "1": 413}) - trace_2 = sim.run()._results[-1] ** 2 + trace_2 = sim.run().states[-1] ** 2 assert np.trace(trace_2) < 1 and not np.isclose(np.trace(trace_2), 1) assert len(sim._collapse_ops) != 0 with pytest.warns(UserWarning, match="first-order"):