From b03567e4c07e60600ae69e185f0b0e2f1414a198 Mon Sep 17 00:00:00 2001 From: a_corni Date: Thu, 26 Oct 2023 10:07:43 +0200 Subject: [PATCH 01/12] Adding Hamiltonian class --- .../pulser_simulation/hamiltonian.py | 670 ++++++++++++++++++ .../pulser_simulation/simulation.py | 600 +--------------- 2 files changed, 695 insertions(+), 575 deletions(-) create mode 100644 pulser-simulation/pulser_simulation/hamiltonian.py diff --git a/pulser-simulation/pulser_simulation/hamiltonian.py b/pulser-simulation/pulser_simulation/hamiltonian.py new file mode 100644 index 000000000..1d3994b99 --- /dev/null +++ b/pulser-simulation/pulser_simulation/hamiltonian.py @@ -0,0 +1,670 @@ +# 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 the QutipEmulator, used to simulate a Sequence or its samples.""" + +from __future__ import annotations + +import itertools +import warnings +from collections import defaultdict +from collections.abc import Mapping +from dataclasses import asdict +from typing import Any, Union, cast + +import numpy as np +import qutip + +from pulser.devices._device_datacls import BaseDevice +from pulser.register.base_register import QubitId +from pulser.sampler.samples import SequenceSamples, _PulseTargetSlot +from pulser_simulation.simconfig import SimConfig + + +def adapt_to_sampling_rate( + full_array: np.ndarray, sampling_rate: float, duration: int +) -> np.ndarray: + """Adapt list to correspond to sampling rate.""" + indices = np.linspace( + 0, + len(full_array) - 1, + int(sampling_rate * duration), + dtype=int, + ) + return cast(np.ndarray, full_array[indices]) + + +class Hamiltonian: + r"""Hamiltonian generated from a sampled sequence and noise. + + Args: + sampled_obj: A sampled pulse sequence extended such that all its + channels have same duration. + device: The device specifications. + sampling_rate: The fraction of samples that we wish to extract from + the samples to simulate. Has to be a value between 0.05 and 1.0. + interaction: The type of interaction between the atoms ("XY", + "ground-rydberg", "raman"). + qdict: A dictionary associating positions to qubit ids. + config: Noise configuration to use. + """ + def __init__( + self, + samples_obj: SequenceSamples, + device: BaseDevice, + sampling_rate: float, + interaction: str, + qdict: dict[QubitId, np.ndarray], + config: SimConfig, + ): + """Instantiates a Hamiltonian object.""" + self.samples_obj: SequenceSamples = samples_obj + self.duration = self.samples_obj.max_duration + self._device = device + self._sampling_rate = sampling_rate + self.sampling_times = adapt_to_sampling_rate( + # Include extra time step for final instruction from samples: + np.arange(self.duration, dtype=np.double) / 1000, + self._sampling_rate, + self.duration, + ) + self._interaction = interaction + self._qdict: dict[QubitId, np.ndarray] = qdict + self._config: SimConfig + + # Type hints for attributes defined outside of __init__ + self.basis_name: str + self.op_matrix: dict[str, qutip.Qobj] + self.basis: dict[str, qutip.Qobj] + self.dim: int + self._eval_times_array: np.ndarray + self._bad_atoms: dict[Union[str, int], bool] = {} + self._doppler_detune: dict[Union[str, int], float] = {} + self.samples: dict + + # Stores the qutip operators used in building the Hamiltonian + self.operators: dict[str, defaultdict[str, dict]] = { + addr: defaultdict(dict) for addr in ["Global", "Local"] + } + self._collapse_ops: list[qutip.Qobj] = [] + + # Initializing qubit infos + self._size = len(self._qdict) + self._qid_index = {qid: i for i, qid in enumerate(self._qdict)} + + # Builds the hamiltonian + self.set_config(config) + + @property + def config(self) -> SimConfig: + """The current configuration, as a SimConfig instance.""" + return self._config + + def build_operator(self, operations: Union[list, tuple]) -> qutip.Qobj: + """Creates an operator with non-trivial actions on some qubits. + + Takes as argument a list of tuples ``[(operator_1, qubits_1), + (operator_2, qubits_2)...]``. Returns the operator given by the tensor + product of {``operator_i`` applied on ``qubits_i``} and Id on the rest. + ``(operator, 'global')`` returns the sum for all ``j`` of operator + applied at ``qubit_j`` and identity elsewhere. + + Example for 4 qubits: ``[(Z, [1, 2]), (Y, [3])]`` returns `ZZYI` + and ``[(X, 'global')]`` returns `XIII + IXII + IIXI + IIIX` + + Args: + operations: List of tuples `(operator, qubits)`. + `operator` can be a ``qutip.Quobj`` or a string key for + ``self.op_matrix``. `qubits` is the list on which operator + will be applied. The qubits can be passed as their + index or their label in the register. + + Returns: + The final operator. + """ + op_list = [self.op_matrix["I"] for j in range(self._size)] + + if not isinstance(operations, list): + operations = [operations] + + for operator, qubits in operations: + if qubits == "global": + return sum( + self.build_operator([(operator, [q_id])]) + for q_id in self._qdict + ) + else: + qubits_set = set(qubits) + if len(qubits_set) < len(qubits): + raise ValueError("Duplicate atom ids in argument list.") + if not qubits_set.issubset(self._qdict.keys()): + raise ValueError( + "Invalid qubit names: " + f"{qubits_set - self._qdict.keys()}" + ) + if isinstance(operator, str): + try: + operator = self.op_matrix[operator] + except KeyError: + raise ValueError(f"{operator} is not a valid operator") + for qubit in qubits: + k = self._qid_index[qubit] + op_list[k] = operator + return qutip.tensor(op_list) + + def _build_collapse_operators(self, prev_config: SimConfig) -> None: + kraus_ops = [] + self._collapse_ops = [] + if "dephasing" in self.config.noise: + if self.basis_name == "digital" or self.basis_name == "all": + # Go back to previous config + self.set_config(prev_config) + raise NotImplementedError( + "Cannot include dephasing noise in digital- or all-basis." + ) + # Probability of phase (Z) flip: + # First order in prob + prob = self.config.dephasing_prob / 2 + n = self._size + if prob > 0.1 and n > 1: + warnings.warn( + "The dephasing model is a first-order approximation in the" + f" dephasing probability. p = {2*prob} is too large for " + "realistic results.", + stacklevel=2, + ) + k = np.sqrt(prob * (1 - prob) ** (n - 1)) + + self._collapse_ops += [ + np.sqrt((1 - prob) ** n) + * qutip.tensor([self.op_matrix["I"] for _ in range(n)]) + ] + kraus_ops.append(k * qutip.sigmaz()) + + if "depolarizing" in self.config.noise: + if self.basis_name == "digital" or self.basis_name == "all": + # Go back to previous config + self.set_config(prev_config) + raise NotImplementedError( + "Cannot include depolarizing " + + "noise in digital- or all-basis." + ) + # Probability of error occurrence + + prob = self.config.depolarizing_prob / 4 + n = self._size + if prob > 0.1 and n > 1: + warnings.warn( + "The depolarizing model is a first-order approximation" + f" in the depolarizing probability. p = {4*prob}" + " is too large for realistic results.", + stacklevel=2, + ) + + k = np.sqrt((prob) * (1 - 3 * prob) ** (n - 1)) + self._collapse_ops += [ + np.sqrt((1 - 3 * prob) ** n) + * qutip.tensor([self.op_matrix["I"] for _ in range(n)]) + ] + kraus_ops.append(k * qutip.sigmax()) + kraus_ops.append(k * qutip.sigmay()) + kraus_ops.append(k * qutip.sigmaz()) + + if "eff_noise" in self.config.noise: + if self.basis_name == "digital" or self.basis_name == "all": + # Go back to previous config + self.set_config(prev_config) + raise NotImplementedError( + "Cannot include general " + + "noise in digital- or all-basis." + ) + # Probability distribution of error occurences + n = self._size + m = len(self.config.eff_noise_opers) + if n > 1: + for i in range(1, m): + prob_i = self.config.eff_noise_probs[i] + if prob_i > 0.1: + warnings.warn( + "The effective noise model is a first-order" + " approximation in the noise probability." + f"p={prob_i} is large for realistic results.", + stacklevel=2, + ) + break + # Deriving Kraus operators + prob_id = self.config.eff_noise_probs[0] + self._collapse_ops += [ + np.sqrt(prob_id**n) + * qutip.tensor([self.op_matrix["I"] for _ in range(n)]) + ] + for i in range(1, m): + k = np.sqrt( + self.config.eff_noise_probs[i] * prob_id ** (n - 1) + ) + k_op = k * self.config.eff_noise_opers[i] + kraus_ops.append(k_op) + + # Building collapse operators + for operator in kraus_ops: + self._collapse_ops += [ + self.build_operator([(operator, [qid])]) + for qid in self._qid_index + ] + + def set_config(self, cfg: SimConfig) -> None: + """Sets current config to cfg and updates simulation parameters. + + Args: + cfg: New configuration. + """ + if not isinstance(cfg, SimConfig): + raise ValueError(f"Object {cfg} is not a valid `SimConfig`.") + not_supported = ( + set(cfg.noise) - cfg.supported_noises[self._interaction] + ) + if not_supported: + raise NotImplementedError( + f"Interaction mode '{self._interaction}' does not support " + f"simulation of noise types: {', '.join(not_supported)}." + ) + prev_config = self.config if hasattr(self, "_config") else SimConfig() + self._config = cfg + if not ("SPAM" in self.config.noise and self.config.eta > 0): + self._bad_atoms = {qid: False for qid in self._qid_index} + if "doppler" not in self.config.noise: + self._doppler_detune = {qid: 0.0 for qid in self._qid_index} + self.generate_hamiltonian() + # Update collapse operators + self._build_collapse_operators(prev_config) + + def add_config(self, config: SimConfig) -> None: + """Updates the current configuration with parameters of another one. + + Mostly useful when dealing with multiple noise types in different + configurations and wanting to merge these configurations together. + Adds simulation parameters to noises that weren't available in the + former SimConfig. Noises specified in both SimConfigs will keep + former noise parameters. + + Args: + config: SimConfig to retrieve parameters from. + """ + if not isinstance(config, SimConfig): + raise ValueError(f"Object {config} is not a valid `SimConfig`") + + not_supported = ( + set(config.noise) - config.supported_noises[self._interaction] + ) + if not_supported: + raise NotImplementedError( + f"Interaction mode '{self._interaction}' does not support " + f"simulation of noise types: {', '.join(not_supported)}." + ) + + old_noise_set = set(self.config.noise) + new_noise_set = old_noise_set.union(config.noise) + diff_noise_set = new_noise_set - old_noise_set + # Create temporary param_dict to add noise parameters: + param_dict: dict[str, Any] = asdict(self._config) + # Begin populating with added noise parameters: + param_dict["noise"] = tuple(new_noise_set) + if "SPAM" in diff_noise_set: + param_dict["eta"] = config.eta + param_dict["epsilon"] = config.epsilon + param_dict["epsilon_prime"] = config.epsilon_prime + if "doppler" in diff_noise_set: + param_dict["temperature"] = config.temperature + if "amplitude" in diff_noise_set: + param_dict["laser_waist"] = config.laser_waist + if "dephasing" in diff_noise_set: + param_dict["dephasing_prob"] = config.dephasing_prob + if "depolarizing" in diff_noise_set: + param_dict["depolarizing_prob"] = config.depolarizing_prob + if "eff_noise" in diff_noise_set: + param_dict["eff_noise_opers"] = config.eff_noise_opers + param_dict["eff_noise_probs"] = config.eff_noise_probs + param_dict["temperature"] *= 1.0e6 + # update runs: + param_dict["runs"] = config.runs + param_dict["samples_per_run"] = config.samples_per_run + + # set config with the new parameters: + self.set_config(SimConfig(**param_dict)) + + def _update_noise(self) -> None: + """Updates noise random parameters. + + Used at the start of each run. If SPAM isn't in chosen noises, all + atoms are set to be correctly prepared. + """ + if "SPAM" in self.config.noise and self.config.eta > 0: + dist = ( + np.random.uniform(size=len(self._qid_index)) + < self.config.spam_dict["eta"] + ) + self._bad_atoms = dict(zip(self._qid_index, dist)) + if "doppler" in self.config.noise: + detune = np.random.normal( + 0, self.config.doppler_sigma, size=len(self._qid_index) + ) + self._doppler_detune = dict(zip(self._qid_index, detune)) + + def _extract_samples(self) -> None: + """Populates samples dictionary with every pulse in the sequence.""" + local_noises = True + if set(self.config.noise).issubset( + {"dephasing", "SPAM", "depolarizing", "eff_noise"} + ): + local_noises = "SPAM" in self.config.noise and self.config.eta > 0 + samples = self.samples_obj.to_nested_dict(all_local=local_noises) + + def add_noise( + slot: _PulseTargetSlot, + samples_dict: Mapping[QubitId, dict[str, np.ndarray]], + is_global_pulse: bool, + ) -> None: + """Builds hamiltonian coefficients. + + Taking into account, if necessary, noise effects, which are local + and depend on the qubit's id qid. + """ + noise_amp_base = max( + 0, np.random.normal(1.0, self.config.amp_sigma) + ) + for qid in slot.targets: + if "doppler" in self.config.noise: + noise_det = self._doppler_detune[qid] + samples_dict[qid]["det"][slot.ti : slot.tf] += noise_det + # Gaussian beam loss in amplitude for global pulses only + # Noise is drawn at random for each pulse + if "amplitude" in self.config.noise and is_global_pulse: + position = self._qdict[qid] + r = np.linalg.norm(position) + w0 = self.config.laser_waist + noise_amp = noise_amp_base * np.exp(-((r / w0) ** 2)) + samples_dict[qid]["amp"][slot.ti : slot.tf] *= noise_amp + + if local_noises: + for ch, ch_samples in self.samples_obj.channel_samples.items(): + addr = self.samples_obj._ch_objs[ch].addressing + basis = self.samples_obj._ch_objs[ch].basis + samples_dict = samples["Local"][basis] + for slot in ch_samples.slots: + add_noise(slot, samples_dict, addr == "Global") + # Delete samples for badly prepared atoms + for basis in samples["Local"]: + for qid in samples["Local"][basis]: + if self._bad_atoms[qid]: + for qty in ("amp", "det", "phase"): + samples["Local"][basis][qid][qty] = 0.0 + self.samples = samples + + def _build_basis_and_op_matrices(self) -> None: + """Determine dimension, basis and projector operators.""" + if self._interaction == "XY": + self.basis_name = "XY" + self.dim = 2 + basis = ["u", "d"] + projectors = ["uu", "du", "ud", "dd"] + else: + if "digital" not in self.samples_obj.used_bases: + self.basis_name = "ground-rydberg" + self.dim = 2 + basis = ["r", "g"] + projectors = ["gr", "rr", "gg"] + elif "ground-rydberg" not in self.samples_obj.used_bases: + self.basis_name = "digital" + self.dim = 2 + basis = ["g", "h"] + projectors = ["hg", "hh", "gg"] + else: + self.basis_name = "all" # All three states + self.dim = 3 + basis = ["r", "g", "h"] + projectors = ["gr", "hg", "rr", "gg", "hh"] + + self.basis = {b: qutip.basis(self.dim, i) for i, b in enumerate(basis)} + self.op_matrix = {"I": qutip.qeye(self.dim)} + + for proj in projectors: + self.op_matrix["sigma_" + proj] = ( + self.basis[proj[0]] * self.basis[proj[1]].dag() + ) + + def generate_hamiltonian(self, update: bool = True) -> None: + """Constructs the hamiltonian from the Sequence. + + Also builds qutip.Qobjs related to the Sequence if not built already, + and refreshes potential noise parameters by drawing new at random. + + Args: + update: Whether to update the noise parameters. + """ + if update: + self._update_noise() + self._extract_samples() + if not hasattr(self, "basis_name"): + self._build_basis_and_op_matrices() + + def make_vdw_term(q1: QubitId, q2: QubitId) -> qutip.Qobj: + """Construct the Van der Waals interaction Term. + + For each pair of qubits, calculate the distance between them, + then assign the local operator "sigma_rr" at each pair. + The units are given so that the coefficient includes a + 1/hbar factor. + """ + dist = np.linalg.norm(self._qdict[q1] - self._qdict[q2]) + U = 0.5 * self._device.interaction_coeff / dist**6 + return U * self.build_operator([("sigma_rr", [q1, q2])]) + + def make_xy_term(q1: QubitId, q2: QubitId) -> qutip.Qobj: + """Construct the XY interaction Term. + + For each pair of qubits, calculate the distance between them, + then assign the local operator "sigma_ud * sigma_du" at each pair. + The units are given so that the coefficient + includes a 1/hbar factor. + """ + dist = np.linalg.norm(self._qdict[q1] - self._qdict[q2]) + coords_dim = len(self._qdict[q1]) + mag_field = cast(np.ndarray, self.samples_obj._magnetic_field)[ + :coords_dim + ] + mag_norm = np.linalg.norm(mag_field) + if mag_norm < 1e-8: + cosine = 0.0 + else: + cosine = np.dot( + (self._qdict[q1] - self._qdict[q2]), + mag_field, + ) / (dist * mag_norm) + U = ( + 0.5 + * cast(float, self._device.interaction_coeff_xy) + * (1 - 3 * cosine**2) + / dist**3 + ) + return U * self.build_operator( + [("sigma_ud", [q1]), ("sigma_du", [q2])] + ) + + def make_interaction_term(masked: bool = False) -> qutip.Qobj: + if masked: + # Calculate the total number of good, unmasked qubits + effective_size = self._size - sum(self._bad_atoms.values()) + for q in self.samples_obj._slm_mask.targets: + if not self._bad_atoms[q]: + effective_size -= 1 + if effective_size < 2: + return 0 * self.build_operator([("I", "global")]) + + # make interaction term + dipole_interaction = cast(qutip.Qobj, 0) + for q1, q2 in itertools.combinations(self._qdict.keys(), r=2): + if ( + self._bad_atoms[q1] + or self._bad_atoms[q2] + or ( + masked + and self._interaction == "XY" + and ( + q1 in self.samples_obj._slm_mask.targets + or q2 in self.samples_obj._slm_mask.targets + ) + ) + ): + continue + + if self._interaction == "XY": + dipole_interaction += make_xy_term(q1, q2) + else: + dipole_interaction += make_vdw_term(q1, q2) + return dipole_interaction + + def build_coeffs_ops(basis: str, addr: str) -> list[list]: + """Build coefficients and operators for the hamiltonian QobjEvo.""" + samples = self.samples[addr][basis] + operators = self.operators[addr][basis] + # Choose operator names according to addressing: + if basis == "ground-rydberg": + op_ids = ["sigma_gr", "sigma_rr"] + elif basis == "digital": + op_ids = ["sigma_hg", "sigma_gg"] + elif basis == "XY": + op_ids = ["sigma_du", "sigma_uu"] + + terms = [] + if addr == "Global": + coeffs = [ + 0.5 * samples["amp"] * np.exp(-1j * samples["phase"]), + -0.5 * samples["det"], + ] + for op_id, coeff in zip(op_ids, coeffs): + if np.any(coeff != 0): + # Build once global operators as they are needed + if op_id not in operators: + operators[op_id] = self.build_operator( + [(op_id, "global")] + ) + terms.append( + [ + operators[op_id], + adapt_to_sampling_rate( + coeff, self._sampling_rate, self.duration + ), + ] + ) + elif addr == "Local": + for q_id, samples_q in samples.items(): + if q_id not in operators: + operators[q_id] = {} + coeffs = [ + 0.5 + * samples_q["amp"] + * np.exp(-1j * samples_q["phase"]), + -0.5 * samples_q["det"], + ] + for coeff, op_id in zip(coeffs, op_ids): + if np.any(coeff != 0): + if op_id not in operators[q_id]: + operators[q_id][op_id] = self.build_operator( + [(op_id, [q_id])] + ) + terms.append( + [ + operators[q_id][op_id], + adapt_to_sampling_rate( + coeff, + self._sampling_rate, + self.duration, + ), + ] + ) + self.operators[addr][basis] = operators + return terms + + qobj_list = [] + # Time independent term: + effective_size = self._size - sum(self._bad_atoms.values()) + if self.basis_name != "digital" and effective_size > 1: + # Build time-dependent or time-independent interaction term based + # on whether an SLM mask was defined or not + if ( + self.samples_obj._slm_mask.end > 0 + and self._interaction == "XY" + ): + # Build an array of binary coefficients for the interaction + # term of unmasked qubits + coeff = np.ones(self.duration - 1) + coeff[0 : self.samples_obj._slm_mask.end] = 0 + # Build the interaction term for unmasked qubits + qobj_list = [ + [ + make_interaction_term(), + adapt_to_sampling_rate( + coeff, self._sampling_rate, self.duration + ), + ] + ] + # Build the interaction term for masked qubits + qobj_list += [ + [ + make_interaction_term(masked=True), + adapt_to_sampling_rate( + np.logical_not(coeff).astype(int), + self._sampling_rate, + self.duration, + ), + ] + ] + else: + qobj_list = [make_interaction_term()] + + # Time dependent terms: + for addr in self.samples: + for basis in self.samples[addr]: + if self.samples[addr][basis]: + qobj_list += cast(list, build_coeffs_ops(basis, addr)) + + if not qobj_list: # If qobj_list ends up empty + qobj_list = [0 * self.build_operator([("I", "global")])] + + ham = qutip.QobjEvo(qobj_list, tlist=self.sampling_times) + ham = ham + ham.dag() + ham.compress() + self._hamiltonian = ham + + def __call__(self, time: float) -> qutip.Qobj: + r"""Get the Hamiltonian created from the sequence at a fixed time. + + Note: + The whole Hamiltonian is divided by :math:`\hbar`, so its + units are rad/µs. + + Args: + time: The specific time at which we want to extract the + Hamiltonian (in ns). + + Returns: + A new Qobj for the Hamiltonian with coefficients + extracted from the effective sequence (determined by + `self.sampling_rate`) at the specified time. + """ + if time < 0: + raise ValueError( + f"Provided time (`time` = {time}) must be " + "greater than or equal to 0." + ) + return self._hamiltonian(time / 1000) # Creates new Qutip.Qobj diff --git a/pulser-simulation/pulser_simulation/simulation.py b/pulser-simulation/pulser_simulation/simulation.py index 1d6c91c7d..59ccbf95b 100644 --- a/pulser-simulation/pulser_simulation/simulation.py +++ b/pulser-simulation/pulser_simulation/simulation.py @@ -15,11 +15,10 @@ from __future__ import annotations -import itertools import warnings -from collections import Counter, defaultdict +from collections import Counter from collections.abc import Mapping -from dataclasses import asdict, replace +from dataclasses import replace from typing import Any, Optional, Union, cast import matplotlib.pyplot as plt @@ -30,10 +29,11 @@ import pulser.sampler as sampler from pulser import Sequence from pulser.devices._device_datacls import BaseDevice -from pulser.register.base_register import BaseRegister, QubitId +from pulser.register.base_register import BaseRegister from pulser.result import SampledResult -from pulser.sampler.samples import SequenceSamples, _PulseTargetSlot +from pulser.sampler.samples import SequenceSamples from pulser.sequence._seq_drawer import draw_samples, draw_sequence +from pulser_simulation.hamiltonian import Hamiltonian, adapt_to_sampling_rate from pulser_simulation.qutip_result import QutipResult from pulser_simulation.simconfig import SimConfig from pulser_simulation.simresults import ( @@ -140,20 +140,8 @@ def __init__( # Initializing qubit infos self._qdict = self._register.qubits - self._size = len(self._qdict) - self._qid_index = {qid: i for i, qid in enumerate(self._qdict)} - - # Type hints for attributes defined outside of __init__ - self.basis_name: str - self._config: SimConfig - self.op_matrix: dict[str, qutip.Qobj] - self.basis: dict[str, qutip.Qobj] - self.dim: int - self._eval_times_array: np.ndarray - self._bad_atoms: dict[Union[str, int], bool] = {} - self._doppler_detune: dict[Union[str, int], float] = {} - - # Initializing sampling and evalutaion times + + # Initializing sampling and evaluation times if not (0 < sampling_rate <= 1.0): raise ValueError( "The sampling rate (`sampling_rate` = " @@ -165,21 +153,23 @@ def __init__( "`sampling_rate` is too small, less than 4 data points." ) self._sampling_rate = sampling_rate - self.sampling_times = self._adapt_to_sampling_rate( + self.sampling_times = adapt_to_sampling_rate( # Include extra time step for final instruction from samples: - np.arange(self._tot_duration + 1, dtype=np.double) - / 1000 + np.arange(self._tot_duration + 1, dtype=np.double) / 1000, + self._sampling_rate, + self._tot_duration + 1, ) self.set_evaluation_times(evaluation_times) - # Stores the qutip operators used in building the Hamiltonian - self.operators: dict[str, defaultdict[str, dict]] = { - addr: defaultdict(dict) for addr in ["Global", "Local"] - } - self._collapse_ops: list[qutip.Qobj] = [] - # Sets the config as well as builds the hamiltonian - self.set_config(config) if config else self.set_config(SimConfig()) + self.hamiltonian = Hamiltonian( + self.samples_obj, + self._device, + self._sampling_rate, + self._interaction, + self._qdict, + config if config else SimConfig(), + ) if self.samples_obj._measurement: self._meas_basis = self.samples_obj._measurement else: @@ -189,192 +179,10 @@ def __init__( self._meas_basis = self.basis_name self.set_initial_state("all-ground") - @property - def config(self) -> SimConfig: - """The current configuration, as a SimConfig instance.""" - return self._config - - def set_config(self, cfg: SimConfig) -> None: - """Sets current config to cfg and updates simulation parameters. - - Args: - cfg: New configuration. - """ - if not isinstance(cfg, SimConfig): - raise ValueError(f"Object {cfg} is not a valid `SimConfig`.") - not_supported = ( - set(cfg.noise) - cfg.supported_noises[self._interaction] - ) - if not_supported: - raise NotImplementedError( - f"Interaction mode '{self._interaction}' does not support " - f"simulation of noise types: {', '.join(not_supported)}." - ) - prev_config = self.config if hasattr(self, "_config") else SimConfig() - self._config = cfg - if not ("SPAM" in self.config.noise and self.config.eta > 0): - self._bad_atoms = {qid: False for qid in self._qid_index} - if "doppler" not in self.config.noise: - self._doppler_detune = {qid: 0.0 for qid in self._qid_index} - # Noise, samples and Hamiltonian update routine - self._construct_hamiltonian() - - kraus_ops = [] - self._collapse_ops = [] - if "dephasing" in self.config.noise: - if self.basis_name == "digital" or self.basis_name == "all": - # Go back to previous config - self.set_config(prev_config) - raise NotImplementedError( - "Cannot include dephasing noise in digital- or all-basis." - ) - # Probability of phase (Z) flip: - # First order in prob - prob = self.config.dephasing_prob / 2 - n = self._size - if prob > 0.1 and n > 1: - warnings.warn( - "The dephasing model is a first-order approximation in the" - f" dephasing probability. p = {2*prob} is too large for " - "realistic results.", - stacklevel=2, - ) - k = np.sqrt(prob * (1 - prob) ** (n - 1)) - - self._collapse_ops += [ - np.sqrt((1 - prob) ** n) - * qutip.tensor([self.op_matrix["I"] for _ in range(n)]) - ] - kraus_ops.append(k * qutip.sigmaz()) - - if "depolarizing" in self.config.noise: - if self.basis_name == "digital" or self.basis_name == "all": - # Go back to previous config - self.set_config(prev_config) - raise NotImplementedError( - "Cannot include depolarizing " - + "noise in digital- or all-basis." - ) - # Probability of error occurrence - - prob = self.config.depolarizing_prob / 4 - n = self._size - if prob > 0.1 and n > 1: - warnings.warn( - "The depolarizing model is a first-order approximation" - f" in the depolarizing probability. p = {4*prob}" - " is too large for realistic results.", - stacklevel=2, - ) - - k = np.sqrt((prob) * (1 - 3 * prob) ** (n - 1)) - self._collapse_ops += [ - np.sqrt((1 - 3 * prob) ** n) - * qutip.tensor([self.op_matrix["I"] for _ in range(n)]) - ] - kraus_ops.append(k * qutip.sigmax()) - kraus_ops.append(k * qutip.sigmay()) - kraus_ops.append(k * qutip.sigmaz()) - - if "eff_noise" in self.config.noise: - if self.basis_name == "digital" or self.basis_name == "all": - # Go back to previous config - self.set_config(prev_config) - raise NotImplementedError( - "Cannot include general " - + "noise in digital- or all-basis." - ) - # Probability distribution of error occurences - n = self._size - m = len(self.config.eff_noise_opers) - if n > 1: - for i in range(1, m): - prob_i = self.config.eff_noise_probs[i] - if prob_i > 0.1: - warnings.warn( - "The effective noise model is a first-order" - " approximation in the noise probability." - f"p={prob_i} is large for realistic results.", - stacklevel=2, - ) - break - # Deriving Kraus operators - prob_id = self.config.eff_noise_probs[0] - self._collapse_ops += [ - np.sqrt(prob_id**n) - * qutip.tensor([self.op_matrix["I"] for _ in range(n)]) - ] - for i in range(1, m): - k = np.sqrt( - self.config.eff_noise_probs[i] * prob_id ** (n - 1) - ) - k_op = k * self.config.eff_noise_opers[i] - kraus_ops.append(k_op) - - # Building collapse operators - for operator in kraus_ops: - self._collapse_ops += [ - self.build_operator([(operator, [qid])]) - for qid in self._qid_index - ] - - def add_config(self, config: SimConfig) -> None: - """Updates the current configuration with parameters of another one. - - Mostly useful when dealing with multiple noise types in different - configurations and wanting to merge these configurations together. - Adds simulation parameters to noises that weren't available in the - former SimConfig. Noises specified in both SimConfigs will keep - former noise parameters. - - Args: - config: SimConfig to retrieve parameters from. - """ - if not isinstance(config, SimConfig): - raise ValueError(f"Object {config} is not a valid `SimConfig`") - - not_supported = ( - set(config.noise) - config.supported_noises[self._interaction] - ) - if not_supported: - raise NotImplementedError( - f"Interaction mode '{self._interaction}' does not support " - f"simulation of noise types: {', '.join(not_supported)}." - ) - - old_noise_set = set(self.config.noise) - new_noise_set = old_noise_set.union(config.noise) - diff_noise_set = new_noise_set - old_noise_set - # Create temporary param_dict to add noise parameters: - param_dict: dict[str, Any] = asdict(self._config) - # Begin populating with added noise parameters: - param_dict["noise"] = tuple(new_noise_set) - if "SPAM" in diff_noise_set: - param_dict["eta"] = config.eta - param_dict["epsilon"] = config.epsilon - param_dict["epsilon_prime"] = config.epsilon_prime - if "doppler" in diff_noise_set: - param_dict["temperature"] = config.temperature - if "amplitude" in diff_noise_set: - param_dict["laser_waist"] = config.laser_waist - if "dephasing" in diff_noise_set: - param_dict["dephasing_prob"] = config.dephasing_prob - if "depolarizing" in diff_noise_set: - param_dict["depolarizing_prob"] = config.depolarizing_prob - if "eff_noise" in diff_noise_set: - param_dict["eff_noise_opers"] = config.eff_noise_opers - param_dict["eff_noise_probs"] = config.eff_noise_probs - param_dict["temperature"] *= 1.0e6 - # update runs: - param_dict["runs"] = config.runs - param_dict["samples_per_run"] = config.samples_per_run - - # set config with the new parameters: - self.set_config(SimConfig(**param_dict)) def show_config(self, solver_options: bool = False) -> None: """Shows current configuration.""" - print(self._config.__str__(solver_options)) + print(self.config.__str__(solver_options)) def reset_config(self) -> None: """Resets configuration to default.""" @@ -490,362 +298,6 @@ def set_evaluation_times( ) self._eval_times_instruction = value - def _extract_samples(self) -> None: - """Populates samples dictionary with every pulse in the sequence.""" - local_noises = True - if set(self.config.noise).issubset( - {"dephasing", "SPAM", "depolarizing", "eff_noise"} - ): - local_noises = "SPAM" in self.config.noise and self.config.eta > 0 - samples = self.samples_obj.to_nested_dict(all_local=local_noises) - - def add_noise( - slot: _PulseTargetSlot, - samples_dict: Mapping[QubitId, dict[str, np.ndarray]], - is_global_pulse: bool, - ) -> None: - """Builds hamiltonian coefficients. - - Taking into account, if necessary, noise effects, which are local - and depend on the qubit's id qid. - """ - noise_amp_base = max( - 0, np.random.normal(1.0, self.config.amp_sigma) - ) - for qid in slot.targets: - if "doppler" in self.config.noise: - noise_det = self._doppler_detune[qid] - samples_dict[qid]["det"][slot.ti : slot.tf] += noise_det - # Gaussian beam loss in amplitude for global pulses only - # Noise is drawn at random for each pulse - if "amplitude" in self.config.noise and is_global_pulse: - position = self._qdict[qid] - r = np.linalg.norm(position) - w0 = self.config.laser_waist - noise_amp = noise_amp_base * np.exp(-((r / w0) ** 2)) - samples_dict[qid]["amp"][slot.ti : slot.tf] *= noise_amp - - if local_noises: - for ch, ch_samples in self.samples_obj.channel_samples.items(): - addr = self.samples_obj._ch_objs[ch].addressing - basis = self.samples_obj._ch_objs[ch].basis - samples_dict = samples["Local"][basis] - for slot in ch_samples.slots: - add_noise(slot, samples_dict, addr == "Global") - # Delete samples for badly prepared atoms - for basis in samples["Local"]: - for qid in samples["Local"][basis]: - if self._bad_atoms[qid]: - for qty in ("amp", "det", "phase"): - samples["Local"][basis][qid][qty] = 0.0 - self.samples = samples - - def build_operator(self, operations: Union[list, tuple]) -> qutip.Qobj: - """Creates an operator with non-trivial actions on some qubits. - - Takes as argument a list of tuples ``[(operator_1, qubits_1), - (operator_2, qubits_2)...]``. Returns the operator given by the tensor - product of {``operator_i`` applied on ``qubits_i``} and Id on the rest. - ``(operator, 'global')`` returns the sum for all ``j`` of operator - applied at ``qubit_j`` and identity elsewhere. - - Example for 4 qubits: ``[(Z, [1, 2]), (Y, [3])]`` returns `ZZYI` - and ``[(X, 'global')]`` returns `XIII + IXII + IIXI + IIIX` - - Args: - operations: List of tuples `(operator, qubits)`. - `operator` can be a ``qutip.Quobj`` or a string key for - ``self.op_matrix``. `qubits` is the list on which operator - will be applied. The qubits can be passed as their - index or their label in the register. - - Returns: - The final operator. - """ - op_list = [self.op_matrix["I"] for j in range(self._size)] - - if not isinstance(operations, list): - operations = [operations] - - for operator, qubits in operations: - if qubits == "global": - return sum( - self.build_operator([(operator, [q_id])]) - for q_id in self._qdict - ) - else: - qubits_set = set(qubits) - if len(qubits_set) < len(qubits): - raise ValueError("Duplicate atom ids in argument list.") - if not qubits_set.issubset(self._qdict.keys()): - raise ValueError( - "Invalid qubit names: " - f"{qubits_set - self._qdict.keys()}" - ) - if isinstance(operator, str): - try: - operator = self.op_matrix[operator] - except KeyError: - raise ValueError(f"{operator} is not a valid operator") - for qubit in qubits: - k = self._qid_index[qubit] - op_list[k] = operator - return qutip.tensor(op_list) - - def _adapt_to_sampling_rate(self, full_array: np.ndarray) -> np.ndarray: - """Adapt list to correspond to sampling rate.""" - indices = np.linspace( - 0, - len(full_array) - 1, - int(self._sampling_rate * (self._tot_duration + 1)), - dtype=int, - ) - return cast(np.ndarray, full_array[indices]) - - def _update_noise(self) -> None: - """Updates noise random parameters. - - Used at the start of each run. If SPAM isn't in chosen noises, all - atoms are set to be correctly prepared. - """ - if "SPAM" in self.config.noise and self.config.eta > 0: - dist = ( - np.random.uniform(size=len(self._qid_index)) - < self.config.spam_dict["eta"] - ) - self._bad_atoms = dict(zip(self._qid_index, dist)) - if "doppler" in self.config.noise: - detune = np.random.normal( - 0, self.config.doppler_sigma, size=len(self._qid_index) - ) - self._doppler_detune = dict(zip(self._qid_index, detune)) - - def _build_basis_and_op_matrices(self) -> None: - """Determine dimension, basis and projector operators.""" - if self._interaction == "XY": - self.basis_name = "XY" - self.dim = 2 - basis = ["u", "d"] - projectors = ["uu", "du", "ud", "dd"] - else: - if "digital" not in self.samples_obj.used_bases: - self.basis_name = "ground-rydberg" - self.dim = 2 - basis = ["r", "g"] - projectors = ["gr", "rr", "gg"] - elif "ground-rydberg" not in self.samples_obj.used_bases: - self.basis_name = "digital" - self.dim = 2 - basis = ["g", "h"] - projectors = ["hg", "hh", "gg"] - else: - self.basis_name = "all" # All three states - self.dim = 3 - basis = ["r", "g", "h"] - projectors = ["gr", "hg", "rr", "gg", "hh"] - - self.basis = {b: qutip.basis(self.dim, i) for i, b in enumerate(basis)} - self.op_matrix = {"I": qutip.qeye(self.dim)} - - for proj in projectors: - self.op_matrix["sigma_" + proj] = ( - self.basis[proj[0]] * self.basis[proj[1]].dag() - ) - - def _construct_hamiltonian(self, update: bool = True) -> None: - """Constructs the hamiltonian from the Sequence. - - Also builds qutip.Qobjs related to the Sequence if not built already, - and refreshes potential noise parameters by drawing new at random. - - Args: - update: Whether to update the noise parameters. - """ - if update: - self._update_noise() - self._extract_samples() - if not hasattr(self, "basis_name"): - self._build_basis_and_op_matrices() - - def make_vdw_term(q1: QubitId, q2: QubitId) -> qutip.Qobj: - """Construct the Van der Waals interaction Term. - - For each pair of qubits, calculate the distance between them, - then assign the local operator "sigma_rr" at each pair. - The units are given so that the coefficient includes a - 1/hbar factor. - """ - dist = np.linalg.norm(self._qdict[q1] - self._qdict[q2]) - U = 0.5 * self._device.interaction_coeff / dist**6 - return U * self.build_operator([("sigma_rr", [q1, q2])]) - - def make_xy_term(q1: QubitId, q2: QubitId) -> qutip.Qobj: - """Construct the XY interaction Term. - - For each pair of qubits, calculate the distance between them, - then assign the local operator "sigma_ud * sigma_du" at each pair. - The units are given so that the coefficient - includes a 1/hbar factor. - """ - dist = np.linalg.norm(self._qdict[q1] - self._qdict[q2]) - coords_dim = len(self._qdict[q1]) - mag_field = cast(np.ndarray, self.samples_obj._magnetic_field)[ - :coords_dim - ] - mag_norm = np.linalg.norm(mag_field) - if mag_norm < 1e-8: - cosine = 0.0 - else: - cosine = np.dot( - (self._qdict[q1] - self._qdict[q2]), - mag_field, - ) / (dist * mag_norm) - U = ( - 0.5 - * cast(float, self._device.interaction_coeff_xy) - * (1 - 3 * cosine**2) - / dist**3 - ) - return U * self.build_operator( - [("sigma_ud", [q1]), ("sigma_du", [q2])] - ) - - def make_interaction_term(masked: bool = False) -> qutip.Qobj: - if masked: - # Calculate the total number of good, unmasked qubits - effective_size = self._size - sum(self._bad_atoms.values()) - for q in self.samples_obj._slm_mask.targets: - if not self._bad_atoms[q]: - effective_size -= 1 - if effective_size < 2: - return 0 * self.build_operator([("I", "global")]) - - # make interaction term - dipole_interaction = cast(qutip.Qobj, 0) - for q1, q2 in itertools.combinations(self._qdict.keys(), r=2): - if ( - self._bad_atoms[q1] - or self._bad_atoms[q2] - or ( - masked - and self._interaction == "XY" - and ( - q1 in self.samples_obj._slm_mask.targets - or q2 in self.samples_obj._slm_mask.targets - ) - ) - ): - continue - - if self._interaction == "XY": - dipole_interaction += make_xy_term(q1, q2) - else: - dipole_interaction += make_vdw_term(q1, q2) - return dipole_interaction - - def build_coeffs_ops(basis: str, addr: str) -> list[list]: - """Build coefficients and operators for the hamiltonian QobjEvo.""" - samples = self.samples[addr][basis] - operators = self.operators[addr][basis] - # Choose operator names according to addressing: - if basis == "ground-rydberg": - op_ids = ["sigma_gr", "sigma_rr"] - elif basis == "digital": - op_ids = ["sigma_hg", "sigma_gg"] - elif basis == "XY": - op_ids = ["sigma_du", "sigma_uu"] - - terms = [] - if addr == "Global": - coeffs = [ - 0.5 * samples["amp"] * np.exp(-1j * samples["phase"]), - -0.5 * samples["det"], - ] - for op_id, coeff in zip(op_ids, coeffs): - if np.any(coeff != 0): - # Build once global operators as they are needed - if op_id not in operators: - operators[op_id] = self.build_operator( - [(op_id, "global")] - ) - terms.append( - [ - operators[op_id], - self._adapt_to_sampling_rate(coeff), - ] - ) - elif addr == "Local": - for q_id, samples_q in samples.items(): - if q_id not in operators: - operators[q_id] = {} - coeffs = [ - 0.5 - * samples_q["amp"] - * np.exp(-1j * samples_q["phase"]), - -0.5 * samples_q["det"], - ] - for coeff, op_id in zip(coeffs, op_ids): - if np.any(coeff != 0): - if op_id not in operators[q_id]: - operators[q_id][op_id] = self.build_operator( - [(op_id, [q_id])] - ) - terms.append( - [ - operators[q_id][op_id], - self._adapt_to_sampling_rate(coeff), - ] - ) - self.operators[addr][basis] = operators - return terms - - qobj_list = [] - # Time independent term: - effective_size = self._size - sum(self._bad_atoms.values()) - if self.basis_name != "digital" and effective_size > 1: - # Build time-dependent or time-independent interaction term based - # on whether an SLM mask was defined or not - if ( - self.samples_obj._slm_mask.end > 0 - and self._interaction == "XY" - ): - # Build an array of binary coefficients for the interaction - # term of unmasked qubits - coeff = np.ones(self._tot_duration) - coeff[0 : self.samples_obj._slm_mask.end] = 0 - # Build the interaction term for unmasked qubits - qobj_list = [ - [ - make_interaction_term(), - self._adapt_to_sampling_rate(coeff), - ] - ] - # Build the interaction term for masked qubits - qobj_list += [ - [ - make_interaction_term(masked=True), - self._adapt_to_sampling_rate( - np.logical_not(coeff).astype(int) - ), - ] - ] - else: - qobj_list = [make_interaction_term()] - - # Time dependent terms: - for addr in self.samples: - for basis in self.samples[addr]: - if self.samples[addr][basis]: - qobj_list += cast(list, build_coeffs_ops(basis, addr)) - - if not qobj_list: # If qobj_list ends up empty - qobj_list = [0 * self.build_operator([("I", "global")])] - - ham = qutip.QobjEvo(qobj_list, tlist=self.sampling_times) - ham = ham + ham.dag() - ham.compress() - self._hamiltonian = ham - def get_hamiltonian(self, time: float) -> qutip.Qobj: r"""Get the Hamiltonian created from the sequence at a fixed time. @@ -868,12 +320,7 @@ def get_hamiltonian(self, time: float) -> qutip.Qobj: "less than or equal to the sequence duration " f"({self._tot_duration})." ) - if time < 0: - raise ValueError( - f"Provided time (`time` = {time}) must be " - "greater than or equal to 0." - ) - return self._hamiltonian(time / 1000) # Creates new Qutip.Qobj + return self.hamiltonian(time) # Run Simulation Evolution using Qutip def run( @@ -1021,7 +468,7 @@ def _run_solver() -> CoherentResults: else: reps = 1 # At each run, new random noise: new Hamiltonian - self._construct_hamiltonian(update=update_ham) + self.hamiltonian.generate_hamiltonian(update=update_ham) # Get CoherentResults instance from sequence with added noise: cleanres_noisyseq = _run_solver() # Extract statistics at eval time: @@ -1156,6 +603,9 @@ def from_sequence( evaluation_times, ) + def __getattr__(self, name: str) -> Any: + return getattr(self.hamiltonian, name) + class Simulation: r"""Simulation of a pulse sequence using QuTiP. From bcc67a5a201819adf7374b259c40e6cecb721b03 Mon Sep 17 00:00:00 2001 From: a_corni Date: Thu, 26 Oct 2023 12:07:11 +0200 Subject: [PATCH 02/12] fixing type --- pulser-simulation/pulser_simulation/hamiltonian.py | 3 ++- pulser-simulation/pulser_simulation/simulation.py | 1 - 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/pulser-simulation/pulser_simulation/hamiltonian.py b/pulser-simulation/pulser_simulation/hamiltonian.py index 1d3994b99..22838b29e 100644 --- a/pulser-simulation/pulser_simulation/hamiltonian.py +++ b/pulser-simulation/pulser_simulation/hamiltonian.py @@ -53,11 +53,12 @@ class Hamiltonian: device: The device specifications. sampling_rate: The fraction of samples that we wish to extract from the samples to simulate. Has to be a value between 0.05 and 1.0. - interaction: The type of interaction between the atoms ("XY", + interaction: The type of interaction between the atoms ("XY", "ground-rydberg", "raman"). qdict: A dictionary associating positions to qubit ids. config: Noise configuration to use. """ + def __init__( self, samples_obj: SequenceSamples, diff --git a/pulser-simulation/pulser_simulation/simulation.py b/pulser-simulation/pulser_simulation/simulation.py index 59ccbf95b..80a12450a 100644 --- a/pulser-simulation/pulser_simulation/simulation.py +++ b/pulser-simulation/pulser_simulation/simulation.py @@ -179,7 +179,6 @@ def __init__( self._meas_basis = self.basis_name self.set_initial_state("all-ground") - def show_config(self, solver_options: bool = False) -> None: """Shows current configuration.""" print(self.config.__str__(solver_options)) From 930223613c25cf50f3bdb0ee8f5835db5d0952d8 Mon Sep 17 00:00:00 2001 From: a_corni Date: Thu, 26 Oct 2023 12:42:12 +0200 Subject: [PATCH 03/12] Modifying docstring --- pulser-simulation/pulser_simulation/hamiltonian.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/pulser-simulation/pulser_simulation/hamiltonian.py b/pulser-simulation/pulser_simulation/hamiltonian.py index 22838b29e..15c58c65b 100644 --- a/pulser-simulation/pulser_simulation/hamiltonian.py +++ b/pulser-simulation/pulser_simulation/hamiltonian.py @@ -444,10 +444,11 @@ def _build_basis_and_op_matrices(self) -> None: ) def generate_hamiltonian(self, update: bool = True) -> None: - """Constructs the hamiltonian from the Sequence. + """Generates the hamiltonian from the sampled sequence and noise. - Also builds qutip.Qobjs related to the Sequence if not built already, - and refreshes potential noise parameters by drawing new at random. + Also builds qutip.Qobjs related to the SequenceSamples if not built + already, and refreshes potential noise parameters by drawing new at + random. Args: update: Whether to update the noise parameters. From 8e43982d51a4a9709ad7466573a1d1dc853003ae Mon Sep 17 00:00:00 2001 From: a_corni Date: Thu, 26 Oct 2023 15:51:51 +0200 Subject: [PATCH 04/12] Simplifying refactoring --- .../pulser_simulation/hamiltonian.py | 653 ++++++++++++++++++ .../pulser_simulation/simulation.py | 597 +--------------- 2 files changed, 674 insertions(+), 576 deletions(-) create mode 100644 pulser-simulation/pulser_simulation/hamiltonian.py diff --git a/pulser-simulation/pulser_simulation/hamiltonian.py b/pulser-simulation/pulser_simulation/hamiltonian.py new file mode 100644 index 000000000..a15c8957f --- /dev/null +++ b/pulser-simulation/pulser_simulation/hamiltonian.py @@ -0,0 +1,653 @@ +# Copyright 2020 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 the Hamiltonian class.""" + +from __future__ import annotations + +import itertools +import warnings +from collections import defaultdict +from collections.abc import Mapping +from dataclasses import asdict +from typing import Any, Union, cast + +import numpy as np +import qutip + +from pulser.devices._device_datacls import BaseDevice +from pulser.register.base_register import QubitId +from pulser.sampler.samples import SequenceSamples, _PulseTargetSlot +from pulser_simulation.simconfig import SimConfig + + +def adapt_to_sampling_rate( + full_array: np.ndarray, sampling_rate: float, duration: int +) -> np.ndarray: + """Adapt list to correspond to sampling rate.""" + indices = np.linspace( + 0, + len(full_array) - 1, + int(sampling_rate * duration), + dtype=int, + ) + return cast(np.ndarray, full_array[indices]) + + +class Hamiltonian: + r"""Generates Hamiltonian from a sampled sequence and noise. + + Args: + samples_obj: A sampled sequence whose ChannelSamples have same + duration. + qdict: A dictionary associating coordinates to qubit ids. + device: The device specifications. + sampling_rate: The fraction of samples that we wish to extract from + the samples to simulate. Has to be a value between 0.05 and 1.0. + config: Configuration to be used for this simulation. + """ + + def __init__( + self, + samples_obj: SequenceSamples, + qdict: dict[QubitId, np.ndarray], + device: BaseDevice, + sampling_rate: float, + config: SimConfig, + ) -> None: + """Instantiates a Simulation object.""" + self.samples_obj = samples_obj + self._qdict = qdict + self._device = device + self._sampling_rate = sampling_rate + + # Type hints for attributes defined outside of __init__ + self.basis_name: str + self._config: SimConfig + self.op_matrix: dict[str, qutip.Qobj] + self.basis: dict[str, qutip.Qobj] + self.dim: int + self._eval_times_array: np.ndarray + self._bad_atoms: dict[Union[str, int], bool] = {} + self._doppler_detune: dict[Union[str, int], float] = {} + + # Define interaction + self._interaction = "XY" if self.samples_obj._in_xy else "ising" + + # Initializing qubit infos + self._size = len(self._qdict) + self._qid_index = {qid: i for i, qid in enumerate(self._qdict)} + + # Compute sampling times + self.duration = self.samples_obj.max_duration + self.sampling_times = adapt_to_sampling_rate( + # Include extra time step for final instruction from samples: + np.arange(self.duration, dtype=np.double) / 1000, + self._sampling_rate, + self.duration, + ) + + # Stores the qutip operators used in building the Hamiltonian + self.operators: dict[str, defaultdict[str, dict]] = { + addr: defaultdict(dict) for addr in ["Global", "Local"] + } + self._collapse_ops: list[qutip.Qobj] = [] + + self.set_config(config) + + @property + def config(self) -> SimConfig: + """The current configuration, as a SimConfig instance.""" + return self._config + + def set_config(self, cfg: SimConfig) -> None: + """Sets current config to cfg and updates simulation parameters. + + Args: + cfg: New configuration. + """ + if not isinstance(cfg, SimConfig): + raise ValueError(f"Object {cfg} is not a valid `SimConfig`.") + not_supported = ( + set(cfg.noise) - cfg.supported_noises[self._interaction] + ) + if not_supported: + raise NotImplementedError( + f"Interaction mode '{self._interaction}' does not support " + f"simulation of noise types: {', '.join(not_supported)}." + ) + prev_config = self.config if hasattr(self, "_config") else SimConfig() + self._config = cfg + if not ("SPAM" in self.config.noise and self.config.eta > 0): + self._bad_atoms = {qid: False for qid in self._qid_index} + if "doppler" not in self.config.noise: + self._doppler_detune = {qid: 0.0 for qid in self._qid_index} + # Noise, samples and Hamiltonian update routine + self._construct_hamiltonian() + + kraus_ops = [] + self._collapse_ops = [] + if "dephasing" in self.config.noise: + if self.basis_name == "digital" or self.basis_name == "all": + # Go back to previous config + self.set_config(prev_config) + raise NotImplementedError( + "Cannot include dephasing noise in digital- or all-basis." + ) + # Probability of phase (Z) flip: + # First order in prob + prob = self.config.dephasing_prob / 2 + n = self._size + if prob > 0.1 and n > 1: + warnings.warn( + "The dephasing model is a first-order approximation in the" + f" dephasing probability. p = {2*prob} is too large for " + "realistic results.", + stacklevel=2, + ) + k = np.sqrt(prob * (1 - prob) ** (n - 1)) + + self._collapse_ops += [ + np.sqrt((1 - prob) ** n) + * qutip.tensor([self.op_matrix["I"] for _ in range(n)]) + ] + kraus_ops.append(k * qutip.sigmaz()) + + if "depolarizing" in self.config.noise: + if self.basis_name == "digital" or self.basis_name == "all": + # Go back to previous config + self.set_config(prev_config) + raise NotImplementedError( + "Cannot include depolarizing " + + "noise in digital- or all-basis." + ) + # Probability of error occurrence + + prob = self.config.depolarizing_prob / 4 + n = self._size + if prob > 0.1 and n > 1: + warnings.warn( + "The depolarizing model is a first-order approximation" + f" in the depolarizing probability. p = {4*prob}" + " is too large for realistic results.", + stacklevel=2, + ) + + k = np.sqrt((prob) * (1 - 3 * prob) ** (n - 1)) + self._collapse_ops += [ + np.sqrt((1 - 3 * prob) ** n) + * qutip.tensor([self.op_matrix["I"] for _ in range(n)]) + ] + kraus_ops.append(k * qutip.sigmax()) + kraus_ops.append(k * qutip.sigmay()) + kraus_ops.append(k * qutip.sigmaz()) + + if "eff_noise" in self.config.noise: + if self.basis_name == "digital" or self.basis_name == "all": + # Go back to previous config + self.set_config(prev_config) + raise NotImplementedError( + "Cannot include general " + + "noise in digital- or all-basis." + ) + # Probability distribution of error occurences + n = self._size + m = len(self.config.eff_noise_opers) + if n > 1: + for i in range(1, m): + prob_i = self.config.eff_noise_probs[i] + if prob_i > 0.1: + warnings.warn( + "The effective noise model is a first-order" + " approximation in the noise probability." + f"p={prob_i} is large for realistic results.", + stacklevel=2, + ) + break + # Deriving Kraus operators + prob_id = self.config.eff_noise_probs[0] + self._collapse_ops += [ + np.sqrt(prob_id**n) + * qutip.tensor([self.op_matrix["I"] for _ in range(n)]) + ] + for i in range(1, m): + k = np.sqrt( + self.config.eff_noise_probs[i] * prob_id ** (n - 1) + ) + k_op = k * self.config.eff_noise_opers[i] + kraus_ops.append(k_op) + + # Building collapse operators + for operator in kraus_ops: + self._collapse_ops += [ + self.build_operator([(operator, [qid])]) + for qid in self._qid_index + ] + + def add_config(self, config: SimConfig) -> None: + """Updates the current configuration with parameters of another one. + + Mostly useful when dealing with multiple noise types in different + configurations and wanting to merge these configurations together. + Adds simulation parameters to noises that weren't available in the + former SimConfig. Noises specified in both SimConfigs will keep + former noise parameters. + + Args: + config: SimConfig to retrieve parameters from. + """ + if not isinstance(config, SimConfig): + raise ValueError(f"Object {config} is not a valid `SimConfig`") + + not_supported = ( + set(config.noise) - config.supported_noises[self._interaction] + ) + if not_supported: + raise NotImplementedError( + f"Interaction mode '{self._interaction}' does not support " + f"simulation of noise types: {', '.join(not_supported)}." + ) + + old_noise_set = set(self.config.noise) + new_noise_set = old_noise_set.union(config.noise) + diff_noise_set = new_noise_set - old_noise_set + # Create temporary param_dict to add noise parameters: + param_dict: dict[str, Any] = asdict(self._config) + # Begin populating with added noise parameters: + param_dict["noise"] = tuple(new_noise_set) + if "SPAM" in diff_noise_set: + param_dict["eta"] = config.eta + param_dict["epsilon"] = config.epsilon + param_dict["epsilon_prime"] = config.epsilon_prime + if "doppler" in diff_noise_set: + param_dict["temperature"] = config.temperature + if "amplitude" in diff_noise_set: + param_dict["laser_waist"] = config.laser_waist + if "dephasing" in diff_noise_set: + param_dict["dephasing_prob"] = config.dephasing_prob + if "depolarizing" in diff_noise_set: + param_dict["depolarizing_prob"] = config.depolarizing_prob + if "eff_noise" in diff_noise_set: + param_dict["eff_noise_opers"] = config.eff_noise_opers + param_dict["eff_noise_probs"] = config.eff_noise_probs + param_dict["temperature"] *= 1.0e6 + # update runs: + param_dict["runs"] = config.runs + param_dict["samples_per_run"] = config.samples_per_run + + # set config with the new parameters: + self.set_config(SimConfig(**param_dict)) + + def show_config(self, solver_options: bool = False) -> None: + """Shows current configuration.""" + print(self._config.__str__(solver_options)) + + def reset_config(self) -> None: + """Resets configuration to default.""" + self.set_config(SimConfig()) + + def _extract_samples(self) -> None: + """Populates samples dictionary with every pulse in the sequence.""" + local_noises = True + if set(self.config.noise).issubset( + {"dephasing", "SPAM", "depolarizing", "eff_noise"} + ): + local_noises = "SPAM" in self.config.noise and self.config.eta > 0 + samples = self.samples_obj.to_nested_dict(all_local=local_noises) + + def add_noise( + slot: _PulseTargetSlot, + samples_dict: Mapping[QubitId, dict[str, np.ndarray]], + is_global_pulse: bool, + ) -> None: + """Builds hamiltonian coefficients. + + Taking into account, if necessary, noise effects, which are local + and depend on the qubit's id qid. + """ + noise_amp_base = max( + 0, np.random.normal(1.0, self.config.amp_sigma) + ) + for qid in slot.targets: + if "doppler" in self.config.noise: + noise_det = self._doppler_detune[qid] + samples_dict[qid]["det"][slot.ti : slot.tf] += noise_det + # Gaussian beam loss in amplitude for global pulses only + # Noise is drawn at random for each pulse + if "amplitude" in self.config.noise and is_global_pulse: + position = self._qdict[qid] + r = np.linalg.norm(position) + w0 = self.config.laser_waist + noise_amp = noise_amp_base * np.exp(-((r / w0) ** 2)) + samples_dict[qid]["amp"][slot.ti : slot.tf] *= noise_amp + + if local_noises: + for ch, ch_samples in self.samples_obj.channel_samples.items(): + addr = self.samples_obj._ch_objs[ch].addressing + basis = self.samples_obj._ch_objs[ch].basis + samples_dict = samples["Local"][basis] + for slot in ch_samples.slots: + add_noise(slot, samples_dict, addr == "Global") + # Delete samples for badly prepared atoms + for basis in samples["Local"]: + for qid in samples["Local"][basis]: + if self._bad_atoms[qid]: + for qty in ("amp", "det", "phase"): + samples["Local"][basis][qid][qty] = 0.0 + self.samples = samples + + def build_operator(self, operations: Union[list, tuple]) -> qutip.Qobj: + """Creates an operator with non-trivial actions on some qubits. + + Takes as argument a list of tuples ``[(operator_1, qubits_1), + (operator_2, qubits_2)...]``. Returns the operator given by the tensor + product of {``operator_i`` applied on ``qubits_i``} and Id on the rest. + ``(operator, 'global')`` returns the sum for all ``j`` of operator + applied at ``qubit_j`` and identity elsewhere. + + Example for 4 qubits: ``[(Z, [1, 2]), (Y, [3])]`` returns `ZZYI` + and ``[(X, 'global')]`` returns `XIII + IXII + IIXI + IIIX` + + Args: + operations: List of tuples `(operator, qubits)`. + `operator` can be a ``qutip.Quobj`` or a string key for + ``self.op_matrix``. `qubits` is the list on which operator + will be applied. The qubits can be passed as their + index or their label in the register. + + Returns: + The final operator. + """ + op_list = [self.op_matrix["I"] for j in range(self._size)] + + if not isinstance(operations, list): + operations = [operations] + + for operator, qubits in operations: + if qubits == "global": + return sum( + self.build_operator([(operator, [q_id])]) + for q_id in self._qdict + ) + else: + qubits_set = set(qubits) + if len(qubits_set) < len(qubits): + raise ValueError("Duplicate atom ids in argument list.") + if not qubits_set.issubset(self._qdict.keys()): + raise ValueError( + "Invalid qubit names: " + f"{qubits_set - self._qdict.keys()}" + ) + if isinstance(operator, str): + try: + operator = self.op_matrix[operator] + except KeyError: + raise ValueError(f"{operator} is not a valid operator") + for qubit in qubits: + k = self._qid_index[qubit] + op_list[k] = operator + return qutip.tensor(op_list) + + def _update_noise(self) -> None: + """Updates noise random parameters. + + Used at the start of each run. If SPAM isn't in chosen noises, all + atoms are set to be correctly prepared. + """ + if "SPAM" in self.config.noise and self.config.eta > 0: + dist = ( + np.random.uniform(size=len(self._qid_index)) + < self.config.spam_dict["eta"] + ) + self._bad_atoms = dict(zip(self._qid_index, dist)) + if "doppler" in self.config.noise: + detune = np.random.normal( + 0, self.config.doppler_sigma, size=len(self._qid_index) + ) + self._doppler_detune = dict(zip(self._qid_index, detune)) + + def _build_basis_and_op_matrices(self) -> None: + """Determine dimension, basis and projector operators.""" + if self._interaction == "XY": + self.basis_name = "XY" + self.dim = 2 + basis = ["u", "d"] + projectors = ["uu", "du", "ud", "dd"] + else: + if "digital" not in self.samples_obj.used_bases: + self.basis_name = "ground-rydberg" + self.dim = 2 + basis = ["r", "g"] + projectors = ["gr", "rr", "gg"] + elif "ground-rydberg" not in self.samples_obj.used_bases: + self.basis_name = "digital" + self.dim = 2 + basis = ["g", "h"] + projectors = ["hg", "hh", "gg"] + else: + self.basis_name = "all" # All three states + self.dim = 3 + basis = ["r", "g", "h"] + projectors = ["gr", "hg", "rr", "gg", "hh"] + + self.basis = {b: qutip.basis(self.dim, i) for i, b in enumerate(basis)} + self.op_matrix = {"I": qutip.qeye(self.dim)} + + for proj in projectors: + self.op_matrix["sigma_" + proj] = ( + self.basis[proj[0]] * self.basis[proj[1]].dag() + ) + + def _construct_hamiltonian(self, update: bool = True) -> None: + """Constructs the hamiltonian from the Sequence. + + Also builds qutip.Qobjs related to the Sequence if not built already, + and refreshes potential noise parameters by drawing new at random. + + Args: + update: Whether to update the noise parameters. + """ + if update: + self._update_noise() + self._extract_samples() + if not hasattr(self, "basis_name"): + self._build_basis_and_op_matrices() + + def make_vdw_term(q1: QubitId, q2: QubitId) -> qutip.Qobj: + """Construct the Van der Waals interaction Term. + + For each pair of qubits, calculate the distance between them, + then assign the local operator "sigma_rr" at each pair. + The units are given so that the coefficient includes a + 1/hbar factor. + """ + dist = np.linalg.norm(self._qdict[q1] - self._qdict[q2]) + U = 0.5 * self._device.interaction_coeff / dist**6 + return U * self.build_operator([("sigma_rr", [q1, q2])]) + + def make_xy_term(q1: QubitId, q2: QubitId) -> qutip.Qobj: + """Construct the XY interaction Term. + + For each pair of qubits, calculate the distance between them, + then assign the local operator "sigma_ud * sigma_du" at each pair. + The units are given so that the coefficient + includes a 1/hbar factor. + """ + dist = np.linalg.norm(self._qdict[q1] - self._qdict[q2]) + coords_dim = len(self._qdict[q1]) + mag_field = cast(np.ndarray, self.samples_obj._magnetic_field)[ + :coords_dim + ] + mag_norm = np.linalg.norm(mag_field) + if mag_norm < 1e-8: + cosine = 0.0 + else: + cosine = np.dot( + (self._qdict[q1] - self._qdict[q2]), + mag_field, + ) / (dist * mag_norm) + U = ( + 0.5 + * cast(float, self._device.interaction_coeff_xy) + * (1 - 3 * cosine**2) + / dist**3 + ) + return U * self.build_operator( + [("sigma_ud", [q1]), ("sigma_du", [q2])] + ) + + def make_interaction_term(masked: bool = False) -> qutip.Qobj: + if masked: + # Calculate the total number of good, unmasked qubits + effective_size = self._size - sum(self._bad_atoms.values()) + for q in self.samples_obj._slm_mask.targets: + if not self._bad_atoms[q]: + effective_size -= 1 + if effective_size < 2: + return 0 * self.build_operator([("I", "global")]) + + # make interaction term + dipole_interaction = cast(qutip.Qobj, 0) + for q1, q2 in itertools.combinations(self._qdict.keys(), r=2): + if ( + self._bad_atoms[q1] + or self._bad_atoms[q2] + or ( + masked + and self._interaction == "XY" + and ( + q1 in self.samples_obj._slm_mask.targets + or q2 in self.samples_obj._slm_mask.targets + ) + ) + ): + continue + + if self._interaction == "XY": + dipole_interaction += make_xy_term(q1, q2) + else: + dipole_interaction += make_vdw_term(q1, q2) + return dipole_interaction + + def build_coeffs_ops(basis: str, addr: str) -> list[list]: + """Build coefficients and operators for the hamiltonian QobjEvo.""" + samples = self.samples[addr][basis] + operators = self.operators[addr][basis] + # Choose operator names according to addressing: + if basis == "ground-rydberg": + op_ids = ["sigma_gr", "sigma_rr"] + elif basis == "digital": + op_ids = ["sigma_hg", "sigma_gg"] + elif basis == "XY": + op_ids = ["sigma_du", "sigma_uu"] + + terms = [] + if addr == "Global": + coeffs = [ + 0.5 * samples["amp"] * np.exp(-1j * samples["phase"]), + -0.5 * samples["det"], + ] + for op_id, coeff in zip(op_ids, coeffs): + if np.any(coeff != 0): + # Build once global operators as they are needed + if op_id not in operators: + operators[op_id] = self.build_operator( + [(op_id, "global")] + ) + terms.append( + [ + operators[op_id], + adapt_to_sampling_rate( + coeff, self._sampling_rate, self.duration + ), + ] + ) + elif addr == "Local": + for q_id, samples_q in samples.items(): + if q_id not in operators: + operators[q_id] = {} + coeffs = [ + 0.5 + * samples_q["amp"] + * np.exp(-1j * samples_q["phase"]), + -0.5 * samples_q["det"], + ] + for coeff, op_id in zip(coeffs, op_ids): + if np.any(coeff != 0): + if op_id not in operators[q_id]: + operators[q_id][op_id] = self.build_operator( + [(op_id, [q_id])] + ) + terms.append( + [ + operators[q_id][op_id], + adapt_to_sampling_rate( + coeff, + self._sampling_rate, + self.duration, + ), + ] + ) + self.operators[addr][basis] = operators + return terms + + qobj_list = [] + # Time independent term: + effective_size = self._size - sum(self._bad_atoms.values()) + if self.basis_name != "digital" and effective_size > 1: + # Build time-dependent or time-independent interaction term based + # on whether an SLM mask was defined or not + if ( + self.samples_obj._slm_mask.end > 0 + and self._interaction == "XY" + ): + # Build an array of binary coefficients for the interaction + # term of unmasked qubits + coeff = np.ones(self.duration - 1) + coeff[0 : self.samples_obj._slm_mask.end] = 0 + # Build the interaction term for unmasked qubits + qobj_list = [ + [ + make_interaction_term(), + adapt_to_sampling_rate( + coeff, self._sampling_rate, self.duration + ), + ] + ] + # Build the interaction term for masked qubits + qobj_list += [ + [ + make_interaction_term(masked=True), + adapt_to_sampling_rate( + np.logical_not(coeff).astype(int), + self._sampling_rate, + self.duration, + ), + ] + ] + else: + qobj_list = [make_interaction_term()] + + # Time dependent terms: + for addr in self.samples: + for basis in self.samples[addr]: + if self.samples[addr][basis]: + qobj_list += cast(list, build_coeffs_ops(basis, addr)) + + if not qobj_list: # If qobj_list ends up empty + qobj_list = [0 * self.build_operator([("I", "global")])] + + ham = qutip.QobjEvo(qobj_list, tlist=self.sampling_times) + ham = ham + ham.dag() + ham.compress() + self._hamiltonian = ham diff --git a/pulser-simulation/pulser_simulation/simulation.py b/pulser-simulation/pulser_simulation/simulation.py index 1d6c91c7d..a32bc6b5f 100644 --- a/pulser-simulation/pulser_simulation/simulation.py +++ b/pulser-simulation/pulser_simulation/simulation.py @@ -15,11 +15,10 @@ from __future__ import annotations -import itertools import warnings -from collections import Counter, defaultdict +from collections import Counter from collections.abc import Mapping -from dataclasses import asdict, replace +from dataclasses import replace from typing import Any, Optional, Union, cast import matplotlib.pyplot as plt @@ -30,10 +29,11 @@ import pulser.sampler as sampler from pulser import Sequence from pulser.devices._device_datacls import BaseDevice -from pulser.register.base_register import BaseRegister, QubitId +from pulser.register.base_register import BaseRegister from pulser.result import SampledResult -from pulser.sampler.samples import SequenceSamples, _PulseTargetSlot +from pulser.sampler.samples import SequenceSamples from pulser.sequence._seq_drawer import draw_samples, draw_sequence +from pulser_simulation.hamiltonian import Hamiltonian, adapt_to_sampling_rate from pulser_simulation.qutip_result import QutipResult from pulser_simulation.simconfig import SimConfig from pulser_simulation.simresults import ( @@ -140,20 +140,8 @@ def __init__( # Initializing qubit infos self._qdict = self._register.qubits - self._size = len(self._qdict) - self._qid_index = {qid: i for i, qid in enumerate(self._qdict)} - - # Type hints for attributes defined outside of __init__ - self.basis_name: str - self._config: SimConfig - self.op_matrix: dict[str, qutip.Qobj] - self.basis: dict[str, qutip.Qobj] - self.dim: int - self._eval_times_array: np.ndarray - self._bad_atoms: dict[Union[str, int], bool] = {} - self._doppler_detune: dict[Union[str, int], float] = {} - - # Initializing sampling and evalutaion times + + # Initializing sampling and evaluation times if not (0 < sampling_rate <= 1.0): raise ValueError( "The sampling rate (`sampling_rate` = " @@ -165,21 +153,22 @@ def __init__( "`sampling_rate` is too small, less than 4 data points." ) self._sampling_rate = sampling_rate - self.sampling_times = self._adapt_to_sampling_rate( + self.sampling_times = adapt_to_sampling_rate( # Include extra time step for final instruction from samples: - np.arange(self._tot_duration + 1, dtype=np.double) - / 1000 + np.arange(self._tot_duration + 1, dtype=np.double) / 1000, + self._sampling_rate, + self._tot_duration + 1, ) self.set_evaluation_times(evaluation_times) - # Stores the qutip operators used in building the Hamiltonian - self.operators: dict[str, defaultdict[str, dict]] = { - addr: defaultdict(dict) for addr in ["Global", "Local"] - } - self._collapse_ops: list[qutip.Qobj] = [] - # Sets the config as well as builds the hamiltonian - self.set_config(config) if config else self.set_config(SimConfig()) + self.hamiltonian = Hamiltonian( + self.samples_obj, + self._qdict, + self._device, + self._sampling_rate, + config if config else SimConfig(), + ) if self.samples_obj._measurement: self._meas_basis = self.samples_obj._measurement else: @@ -189,197 +178,6 @@ def __init__( self._meas_basis = self.basis_name self.set_initial_state("all-ground") - @property - def config(self) -> SimConfig: - """The current configuration, as a SimConfig instance.""" - return self._config - - def set_config(self, cfg: SimConfig) -> None: - """Sets current config to cfg and updates simulation parameters. - - Args: - cfg: New configuration. - """ - if not isinstance(cfg, SimConfig): - raise ValueError(f"Object {cfg} is not a valid `SimConfig`.") - not_supported = ( - set(cfg.noise) - cfg.supported_noises[self._interaction] - ) - if not_supported: - raise NotImplementedError( - f"Interaction mode '{self._interaction}' does not support " - f"simulation of noise types: {', '.join(not_supported)}." - ) - prev_config = self.config if hasattr(self, "_config") else SimConfig() - self._config = cfg - if not ("SPAM" in self.config.noise and self.config.eta > 0): - self._bad_atoms = {qid: False for qid in self._qid_index} - if "doppler" not in self.config.noise: - self._doppler_detune = {qid: 0.0 for qid in self._qid_index} - # Noise, samples and Hamiltonian update routine - self._construct_hamiltonian() - - kraus_ops = [] - self._collapse_ops = [] - if "dephasing" in self.config.noise: - if self.basis_name == "digital" or self.basis_name == "all": - # Go back to previous config - self.set_config(prev_config) - raise NotImplementedError( - "Cannot include dephasing noise in digital- or all-basis." - ) - # Probability of phase (Z) flip: - # First order in prob - prob = self.config.dephasing_prob / 2 - n = self._size - if prob > 0.1 and n > 1: - warnings.warn( - "The dephasing model is a first-order approximation in the" - f" dephasing probability. p = {2*prob} is too large for " - "realistic results.", - stacklevel=2, - ) - k = np.sqrt(prob * (1 - prob) ** (n - 1)) - - self._collapse_ops += [ - np.sqrt((1 - prob) ** n) - * qutip.tensor([self.op_matrix["I"] for _ in range(n)]) - ] - kraus_ops.append(k * qutip.sigmaz()) - - if "depolarizing" in self.config.noise: - if self.basis_name == "digital" or self.basis_name == "all": - # Go back to previous config - self.set_config(prev_config) - raise NotImplementedError( - "Cannot include depolarizing " - + "noise in digital- or all-basis." - ) - # Probability of error occurrence - - prob = self.config.depolarizing_prob / 4 - n = self._size - if prob > 0.1 and n > 1: - warnings.warn( - "The depolarizing model is a first-order approximation" - f" in the depolarizing probability. p = {4*prob}" - " is too large for realistic results.", - stacklevel=2, - ) - - k = np.sqrt((prob) * (1 - 3 * prob) ** (n - 1)) - self._collapse_ops += [ - np.sqrt((1 - 3 * prob) ** n) - * qutip.tensor([self.op_matrix["I"] for _ in range(n)]) - ] - kraus_ops.append(k * qutip.sigmax()) - kraus_ops.append(k * qutip.sigmay()) - kraus_ops.append(k * qutip.sigmaz()) - - if "eff_noise" in self.config.noise: - if self.basis_name == "digital" or self.basis_name == "all": - # Go back to previous config - self.set_config(prev_config) - raise NotImplementedError( - "Cannot include general " - + "noise in digital- or all-basis." - ) - # Probability distribution of error occurences - n = self._size - m = len(self.config.eff_noise_opers) - if n > 1: - for i in range(1, m): - prob_i = self.config.eff_noise_probs[i] - if prob_i > 0.1: - warnings.warn( - "The effective noise model is a first-order" - " approximation in the noise probability." - f"p={prob_i} is large for realistic results.", - stacklevel=2, - ) - break - # Deriving Kraus operators - prob_id = self.config.eff_noise_probs[0] - self._collapse_ops += [ - np.sqrt(prob_id**n) - * qutip.tensor([self.op_matrix["I"] for _ in range(n)]) - ] - for i in range(1, m): - k = np.sqrt( - self.config.eff_noise_probs[i] * prob_id ** (n - 1) - ) - k_op = k * self.config.eff_noise_opers[i] - kraus_ops.append(k_op) - - # Building collapse operators - for operator in kraus_ops: - self._collapse_ops += [ - self.build_operator([(operator, [qid])]) - for qid in self._qid_index - ] - - def add_config(self, config: SimConfig) -> None: - """Updates the current configuration with parameters of another one. - - Mostly useful when dealing with multiple noise types in different - configurations and wanting to merge these configurations together. - Adds simulation parameters to noises that weren't available in the - former SimConfig. Noises specified in both SimConfigs will keep - former noise parameters. - - Args: - config: SimConfig to retrieve parameters from. - """ - if not isinstance(config, SimConfig): - raise ValueError(f"Object {config} is not a valid `SimConfig`") - - not_supported = ( - set(config.noise) - config.supported_noises[self._interaction] - ) - if not_supported: - raise NotImplementedError( - f"Interaction mode '{self._interaction}' does not support " - f"simulation of noise types: {', '.join(not_supported)}." - ) - - old_noise_set = set(self.config.noise) - new_noise_set = old_noise_set.union(config.noise) - diff_noise_set = new_noise_set - old_noise_set - # Create temporary param_dict to add noise parameters: - param_dict: dict[str, Any] = asdict(self._config) - # Begin populating with added noise parameters: - param_dict["noise"] = tuple(new_noise_set) - if "SPAM" in diff_noise_set: - param_dict["eta"] = config.eta - param_dict["epsilon"] = config.epsilon - param_dict["epsilon_prime"] = config.epsilon_prime - if "doppler" in diff_noise_set: - param_dict["temperature"] = config.temperature - if "amplitude" in diff_noise_set: - param_dict["laser_waist"] = config.laser_waist - if "dephasing" in diff_noise_set: - param_dict["dephasing_prob"] = config.dephasing_prob - if "depolarizing" in diff_noise_set: - param_dict["depolarizing_prob"] = config.depolarizing_prob - if "eff_noise" in diff_noise_set: - param_dict["eff_noise_opers"] = config.eff_noise_opers - param_dict["eff_noise_probs"] = config.eff_noise_probs - param_dict["temperature"] *= 1.0e6 - # update runs: - param_dict["runs"] = config.runs - param_dict["samples_per_run"] = config.samples_per_run - - # set config with the new parameters: - self.set_config(SimConfig(**param_dict)) - - def show_config(self, solver_options: bool = False) -> None: - """Shows current configuration.""" - print(self._config.__str__(solver_options)) - - def reset_config(self) -> None: - """Resets configuration to default.""" - self.set_config(SimConfig()) - @property def initial_state(self) -> qutip.Qobj: """The initial state of the simulation. @@ -490,362 +288,6 @@ def set_evaluation_times( ) self._eval_times_instruction = value - def _extract_samples(self) -> None: - """Populates samples dictionary with every pulse in the sequence.""" - local_noises = True - if set(self.config.noise).issubset( - {"dephasing", "SPAM", "depolarizing", "eff_noise"} - ): - local_noises = "SPAM" in self.config.noise and self.config.eta > 0 - samples = self.samples_obj.to_nested_dict(all_local=local_noises) - - def add_noise( - slot: _PulseTargetSlot, - samples_dict: Mapping[QubitId, dict[str, np.ndarray]], - is_global_pulse: bool, - ) -> None: - """Builds hamiltonian coefficients. - - Taking into account, if necessary, noise effects, which are local - and depend on the qubit's id qid. - """ - noise_amp_base = max( - 0, np.random.normal(1.0, self.config.amp_sigma) - ) - for qid in slot.targets: - if "doppler" in self.config.noise: - noise_det = self._doppler_detune[qid] - samples_dict[qid]["det"][slot.ti : slot.tf] += noise_det - # Gaussian beam loss in amplitude for global pulses only - # Noise is drawn at random for each pulse - if "amplitude" in self.config.noise and is_global_pulse: - position = self._qdict[qid] - r = np.linalg.norm(position) - w0 = self.config.laser_waist - noise_amp = noise_amp_base * np.exp(-((r / w0) ** 2)) - samples_dict[qid]["amp"][slot.ti : slot.tf] *= noise_amp - - if local_noises: - for ch, ch_samples in self.samples_obj.channel_samples.items(): - addr = self.samples_obj._ch_objs[ch].addressing - basis = self.samples_obj._ch_objs[ch].basis - samples_dict = samples["Local"][basis] - for slot in ch_samples.slots: - add_noise(slot, samples_dict, addr == "Global") - # Delete samples for badly prepared atoms - for basis in samples["Local"]: - for qid in samples["Local"][basis]: - if self._bad_atoms[qid]: - for qty in ("amp", "det", "phase"): - samples["Local"][basis][qid][qty] = 0.0 - self.samples = samples - - def build_operator(self, operations: Union[list, tuple]) -> qutip.Qobj: - """Creates an operator with non-trivial actions on some qubits. - - Takes as argument a list of tuples ``[(operator_1, qubits_1), - (operator_2, qubits_2)...]``. Returns the operator given by the tensor - product of {``operator_i`` applied on ``qubits_i``} and Id on the rest. - ``(operator, 'global')`` returns the sum for all ``j`` of operator - applied at ``qubit_j`` and identity elsewhere. - - Example for 4 qubits: ``[(Z, [1, 2]), (Y, [3])]`` returns `ZZYI` - and ``[(X, 'global')]`` returns `XIII + IXII + IIXI + IIIX` - - Args: - operations: List of tuples `(operator, qubits)`. - `operator` can be a ``qutip.Quobj`` or a string key for - ``self.op_matrix``. `qubits` is the list on which operator - will be applied. The qubits can be passed as their - index or their label in the register. - - Returns: - The final operator. - """ - op_list = [self.op_matrix["I"] for j in range(self._size)] - - if not isinstance(operations, list): - operations = [operations] - - for operator, qubits in operations: - if qubits == "global": - return sum( - self.build_operator([(operator, [q_id])]) - for q_id in self._qdict - ) - else: - qubits_set = set(qubits) - if len(qubits_set) < len(qubits): - raise ValueError("Duplicate atom ids in argument list.") - if not qubits_set.issubset(self._qdict.keys()): - raise ValueError( - "Invalid qubit names: " - f"{qubits_set - self._qdict.keys()}" - ) - if isinstance(operator, str): - try: - operator = self.op_matrix[operator] - except KeyError: - raise ValueError(f"{operator} is not a valid operator") - for qubit in qubits: - k = self._qid_index[qubit] - op_list[k] = operator - return qutip.tensor(op_list) - - def _adapt_to_sampling_rate(self, full_array: np.ndarray) -> np.ndarray: - """Adapt list to correspond to sampling rate.""" - indices = np.linspace( - 0, - len(full_array) - 1, - int(self._sampling_rate * (self._tot_duration + 1)), - dtype=int, - ) - return cast(np.ndarray, full_array[indices]) - - def _update_noise(self) -> None: - """Updates noise random parameters. - - Used at the start of each run. If SPAM isn't in chosen noises, all - atoms are set to be correctly prepared. - """ - if "SPAM" in self.config.noise and self.config.eta > 0: - dist = ( - np.random.uniform(size=len(self._qid_index)) - < self.config.spam_dict["eta"] - ) - self._bad_atoms = dict(zip(self._qid_index, dist)) - if "doppler" in self.config.noise: - detune = np.random.normal( - 0, self.config.doppler_sigma, size=len(self._qid_index) - ) - self._doppler_detune = dict(zip(self._qid_index, detune)) - - def _build_basis_and_op_matrices(self) -> None: - """Determine dimension, basis and projector operators.""" - if self._interaction == "XY": - self.basis_name = "XY" - self.dim = 2 - basis = ["u", "d"] - projectors = ["uu", "du", "ud", "dd"] - else: - if "digital" not in self.samples_obj.used_bases: - self.basis_name = "ground-rydberg" - self.dim = 2 - basis = ["r", "g"] - projectors = ["gr", "rr", "gg"] - elif "ground-rydberg" not in self.samples_obj.used_bases: - self.basis_name = "digital" - self.dim = 2 - basis = ["g", "h"] - projectors = ["hg", "hh", "gg"] - else: - self.basis_name = "all" # All three states - self.dim = 3 - basis = ["r", "g", "h"] - projectors = ["gr", "hg", "rr", "gg", "hh"] - - self.basis = {b: qutip.basis(self.dim, i) for i, b in enumerate(basis)} - self.op_matrix = {"I": qutip.qeye(self.dim)} - - for proj in projectors: - self.op_matrix["sigma_" + proj] = ( - self.basis[proj[0]] * self.basis[proj[1]].dag() - ) - - def _construct_hamiltonian(self, update: bool = True) -> None: - """Constructs the hamiltonian from the Sequence. - - Also builds qutip.Qobjs related to the Sequence if not built already, - and refreshes potential noise parameters by drawing new at random. - - Args: - update: Whether to update the noise parameters. - """ - if update: - self._update_noise() - self._extract_samples() - if not hasattr(self, "basis_name"): - self._build_basis_and_op_matrices() - - def make_vdw_term(q1: QubitId, q2: QubitId) -> qutip.Qobj: - """Construct the Van der Waals interaction Term. - - For each pair of qubits, calculate the distance between them, - then assign the local operator "sigma_rr" at each pair. - The units are given so that the coefficient includes a - 1/hbar factor. - """ - dist = np.linalg.norm(self._qdict[q1] - self._qdict[q2]) - U = 0.5 * self._device.interaction_coeff / dist**6 - return U * self.build_operator([("sigma_rr", [q1, q2])]) - - def make_xy_term(q1: QubitId, q2: QubitId) -> qutip.Qobj: - """Construct the XY interaction Term. - - For each pair of qubits, calculate the distance between them, - then assign the local operator "sigma_ud * sigma_du" at each pair. - The units are given so that the coefficient - includes a 1/hbar factor. - """ - dist = np.linalg.norm(self._qdict[q1] - self._qdict[q2]) - coords_dim = len(self._qdict[q1]) - mag_field = cast(np.ndarray, self.samples_obj._magnetic_field)[ - :coords_dim - ] - mag_norm = np.linalg.norm(mag_field) - if mag_norm < 1e-8: - cosine = 0.0 - else: - cosine = np.dot( - (self._qdict[q1] - self._qdict[q2]), - mag_field, - ) / (dist * mag_norm) - U = ( - 0.5 - * cast(float, self._device.interaction_coeff_xy) - * (1 - 3 * cosine**2) - / dist**3 - ) - return U * self.build_operator( - [("sigma_ud", [q1]), ("sigma_du", [q2])] - ) - - def make_interaction_term(masked: bool = False) -> qutip.Qobj: - if masked: - # Calculate the total number of good, unmasked qubits - effective_size = self._size - sum(self._bad_atoms.values()) - for q in self.samples_obj._slm_mask.targets: - if not self._bad_atoms[q]: - effective_size -= 1 - if effective_size < 2: - return 0 * self.build_operator([("I", "global")]) - - # make interaction term - dipole_interaction = cast(qutip.Qobj, 0) - for q1, q2 in itertools.combinations(self._qdict.keys(), r=2): - if ( - self._bad_atoms[q1] - or self._bad_atoms[q2] - or ( - masked - and self._interaction == "XY" - and ( - q1 in self.samples_obj._slm_mask.targets - or q2 in self.samples_obj._slm_mask.targets - ) - ) - ): - continue - - if self._interaction == "XY": - dipole_interaction += make_xy_term(q1, q2) - else: - dipole_interaction += make_vdw_term(q1, q2) - return dipole_interaction - - def build_coeffs_ops(basis: str, addr: str) -> list[list]: - """Build coefficients and operators for the hamiltonian QobjEvo.""" - samples = self.samples[addr][basis] - operators = self.operators[addr][basis] - # Choose operator names according to addressing: - if basis == "ground-rydberg": - op_ids = ["sigma_gr", "sigma_rr"] - elif basis == "digital": - op_ids = ["sigma_hg", "sigma_gg"] - elif basis == "XY": - op_ids = ["sigma_du", "sigma_uu"] - - terms = [] - if addr == "Global": - coeffs = [ - 0.5 * samples["amp"] * np.exp(-1j * samples["phase"]), - -0.5 * samples["det"], - ] - for op_id, coeff in zip(op_ids, coeffs): - if np.any(coeff != 0): - # Build once global operators as they are needed - if op_id not in operators: - operators[op_id] = self.build_operator( - [(op_id, "global")] - ) - terms.append( - [ - operators[op_id], - self._adapt_to_sampling_rate(coeff), - ] - ) - elif addr == "Local": - for q_id, samples_q in samples.items(): - if q_id not in operators: - operators[q_id] = {} - coeffs = [ - 0.5 - * samples_q["amp"] - * np.exp(-1j * samples_q["phase"]), - -0.5 * samples_q["det"], - ] - for coeff, op_id in zip(coeffs, op_ids): - if np.any(coeff != 0): - if op_id not in operators[q_id]: - operators[q_id][op_id] = self.build_operator( - [(op_id, [q_id])] - ) - terms.append( - [ - operators[q_id][op_id], - self._adapt_to_sampling_rate(coeff), - ] - ) - self.operators[addr][basis] = operators - return terms - - qobj_list = [] - # Time independent term: - effective_size = self._size - sum(self._bad_atoms.values()) - if self.basis_name != "digital" and effective_size > 1: - # Build time-dependent or time-independent interaction term based - # on whether an SLM mask was defined or not - if ( - self.samples_obj._slm_mask.end > 0 - and self._interaction == "XY" - ): - # Build an array of binary coefficients for the interaction - # term of unmasked qubits - coeff = np.ones(self._tot_duration) - coeff[0 : self.samples_obj._slm_mask.end] = 0 - # Build the interaction term for unmasked qubits - qobj_list = [ - [ - make_interaction_term(), - self._adapt_to_sampling_rate(coeff), - ] - ] - # Build the interaction term for masked qubits - qobj_list += [ - [ - make_interaction_term(masked=True), - self._adapt_to_sampling_rate( - np.logical_not(coeff).astype(int) - ), - ] - ] - else: - qobj_list = [make_interaction_term()] - - # Time dependent terms: - for addr in self.samples: - for basis in self.samples[addr]: - if self.samples[addr][basis]: - qobj_list += cast(list, build_coeffs_ops(basis, addr)) - - if not qobj_list: # If qobj_list ends up empty - qobj_list = [0 * self.build_operator([("I", "global")])] - - ham = qutip.QobjEvo(qobj_list, tlist=self.sampling_times) - ham = ham + ham.dag() - ham.compress() - self._hamiltonian = ham - def get_hamiltonian(self, time: float) -> qutip.Qobj: r"""Get the Hamiltonian created from the sequence at a fixed time. @@ -1156,6 +598,9 @@ def from_sequence( evaluation_times, ) + def __getattr__(self, name: str) -> Any: + return getattr(self.hamiltonian, name) + class Simulation: r"""Simulation of a pulse sequence using QuTiP. From 15d8004ef9233e9cfe3222e8398554af70cf9e33 Mon Sep 17 00:00:00 2001 From: a_corni Date: Thu, 26 Oct 2023 16:12:54 +0200 Subject: [PATCH 05/12] Simplifying refactoring --- pulser-simulation/pulser_simulation/simulation.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/pulser-simulation/pulser_simulation/simulation.py b/pulser-simulation/pulser_simulation/simulation.py index 103bb02d6..a32bc6b5f 100644 --- a/pulser-simulation/pulser_simulation/simulation.py +++ b/pulser-simulation/pulser_simulation/simulation.py @@ -310,7 +310,12 @@ def get_hamiltonian(self, time: float) -> qutip.Qobj: "less than or equal to the sequence duration " f"({self._tot_duration})." ) - return self.hamiltonian(time) + if time < 0: + raise ValueError( + f"Provided time (`time` = {time}) must be " + "greater than or equal to 0." + ) + return self._hamiltonian(time / 1000) # Creates new Qutip.Qobj # Run Simulation Evolution using Qutip def run( @@ -458,7 +463,7 @@ def _run_solver() -> CoherentResults: else: reps = 1 # At each run, new random noise: new Hamiltonian - self.hamiltonian.generate_hamiltonian(update=update_ham) + self._construct_hamiltonian(update=update_ham) # Get CoherentResults instance from sequence with added noise: cleanres_noisyseq = _run_solver() # Extract statistics at eval time: From 9f15ba77d950991715ba51decd3f2d6ac29081bb Mon Sep 17 00:00:00 2001 From: a_corni Date: Thu, 26 Oct 2023 16:46:54 +0200 Subject: [PATCH 06/12] Fix test, simplify attributes in QutipEmulator --- .../pulser_simulation/hamiltonian.py | 1 - .../pulser_simulation/simulation.py | 31 ++++++------------- tests/test_simulation.py | 4 +-- 3 files changed, 12 insertions(+), 24 deletions(-) diff --git a/pulser-simulation/pulser_simulation/hamiltonian.py b/pulser-simulation/pulser_simulation/hamiltonian.py index c7efc641c..c522f9e4e 100644 --- a/pulser-simulation/pulser_simulation/hamiltonian.py +++ b/pulser-simulation/pulser_simulation/hamiltonian.py @@ -77,7 +77,6 @@ def __init__( self.op_matrix: dict[str, qutip.Qobj] self.basis: dict[str, qutip.Qobj] self.dim: int - self._eval_times_array: np.ndarray self._bad_atoms: dict[Union[str, int], bool] = {} self._doppler_detune: dict[Union[str, int], float] = {} diff --git a/pulser-simulation/pulser_simulation/simulation.py b/pulser-simulation/pulser_simulation/simulation.py index a32bc6b5f..2193347e8 100644 --- a/pulser-simulation/pulser_simulation/simulation.py +++ b/pulser-simulation/pulser_simulation/simulation.py @@ -88,13 +88,12 @@ def __init__( if sampled_seq.max_duration == 0: raise ValueError("SequenceSamples is empty.") # Check compatibility of register and device - self._device = device - self._device.validate_register(register) + device.validate_register(register) self._register = register # Check compatibility of samples and device: if ( sampled_seq._slm_mask.end > 0 - and not self._device.supports_slm_mask + and not device.supports_slm_mask ): raise ValueError( "Samples use SLM mask but device does not have one." @@ -134,14 +133,10 @@ def __init__( ) ) _sampled_seq = replace(sampled_seq, samples_list=samples_list) - self._interaction = "XY" if _sampled_seq._in_xy else "ising" self._tot_duration = _sampled_seq.max_duration self.samples_obj = _sampled_seq.extend_duration(self._tot_duration + 1) - # Initializing qubit infos - self._qdict = self._register.qubits - - # Initializing sampling and evaluation times + # Testing sampling if not (0 < sampling_rate <= 1.0): raise ValueError( "The sampling rate (`sampling_rate` = " @@ -152,23 +147,18 @@ def __init__( raise ValueError( "`sampling_rate` is too small, less than 4 data points." ) - self._sampling_rate = sampling_rate - self.sampling_times = adapt_to_sampling_rate( - # Include extra time step for final instruction from samples: - np.arange(self._tot_duration + 1, dtype=np.double) / 1000, - self._sampling_rate, - self._tot_duration + 1, - ) - self.set_evaluation_times(evaluation_times) - # Sets the config as well as builds the hamiltonian self.hamiltonian = Hamiltonian( self.samples_obj, - self._qdict, - self._device, - self._sampling_rate, + self._register.qubits, + device, + sampling_rate, config if config else SimConfig(), ) + # Initializing evaluation times + self._eval_times_array: np.ndarray + self.set_evaluation_times(evaluation_times) + if self.samples_obj._measurement: self._meas_basis = self.samples_obj._measurement else: @@ -286,7 +276,6 @@ def set_evaluation_times( self._eval_times_array = np.union1d( eval_times, [0.0, self._tot_duration / 1000] ) - self._eval_times_instruction = value def get_hamiltonian(self, time: float) -> qutip.Qobj: r"""Get the Hamiltonian created from the sequence at a fixed time. diff --git a/tests/test_simulation.py b/tests/test_simulation.py index 27c7c9e99..f947ee43a 100644 --- a/tests/test_simulation.py +++ b/tests/test_simulation.py @@ -663,12 +663,12 @@ def test_config(): def test_noise(seq, matrices): - np.random.seed(3) + np.random.seed(2) sim2 = QutipEmulator.from_sequence( seq, sampling_rate=0.01, config=SimConfig(noise=("SPAM"), eta=0.9) ) assert sim2.run().sample_final_state() == Counter( - {"000": 857, "110": 73, "100": 70} + {"000": 955, "001": 6, "100": 39} ) with pytest.raises(NotImplementedError, match="Cannot include"): sim2.set_config(SimConfig(noise="dephasing")) From 42b1172f45acab0e45ae2b00275b4bc2c4910338 Mon Sep 17 00:00:00 2001 From: a_corni Date: Thu, 26 Oct 2023 16:59:53 +0200 Subject: [PATCH 07/12] Restore _eval_times_instruction, def device in ham --- .../pulser_simulation/hamiltonian.py | 48 +++++++------------ .../pulser_simulation/simulation.py | 8 ++-- 2 files changed, 20 insertions(+), 36 deletions(-) diff --git a/pulser-simulation/pulser_simulation/hamiltonian.py b/pulser-simulation/pulser_simulation/hamiltonian.py index c522f9e4e..bc568b0f4 100644 --- a/pulser-simulation/pulser_simulation/hamiltonian.py +++ b/pulser-simulation/pulser_simulation/hamiltonian.py @@ -31,19 +31,6 @@ from pulser_simulation.simconfig import SimConfig -def adapt_to_sampling_rate( - full_array: np.ndarray, sampling_rate: float, duration: int -) -> np.ndarray: - """Adapt list to correspond to sampling rate.""" - indices = np.linspace( - 0, - len(full_array) - 1, - int(sampling_rate * duration), - dtype=int, - ) - return cast(np.ndarray, full_array[indices]) - - class Hamiltonian: r"""Generates Hamiltonian from a sampled sequence and noise. @@ -89,11 +76,10 @@ def __init__( # Compute sampling times self.duration = self.samples_obj.max_duration - self.sampling_times = adapt_to_sampling_rate( + self.sampling_times = self._adapt_to_sampling_rate( # Include extra time step for final instruction from samples: - np.arange(self.duration, dtype=np.double) / 1000, - self._sampling_rate, - self.duration, + np.arange(self.duration, dtype=np.double) + / 1000 ) # Stores the qutip operators used in building the Hamiltonian @@ -104,6 +90,16 @@ def __init__( self.set_config(config) + def _adapt_to_sampling_rate(self, full_array: np.ndarray) -> np.ndarray: + """Adapt list to correspond to sampling rate.""" + indices = np.linspace( + 0, + len(full_array) - 1, + int(self._sampling_rate * self.duration), + dtype=int, + ) + return cast(np.ndarray, full_array[indices]) + @property def config(self) -> SimConfig: """The current configuration, as a SimConfig instance.""" @@ -568,9 +564,7 @@ def build_coeffs_ops(basis: str, addr: str) -> list[list]: terms.append( [ operators[op_id], - adapt_to_sampling_rate( - coeff, self._sampling_rate, self.duration - ), + self._adapt_to_sampling_rate(coeff), ] ) elif addr == "Local": @@ -592,11 +586,7 @@ def build_coeffs_ops(basis: str, addr: str) -> list[list]: terms.append( [ operators[q_id][op_id], - adapt_to_sampling_rate( - coeff, - self._sampling_rate, - self.duration, - ), + self._adapt_to_sampling_rate(coeff), ] ) self.operators[addr][basis] = operators @@ -620,19 +610,15 @@ def build_coeffs_ops(basis: str, addr: str) -> list[list]: qobj_list = [ [ make_interaction_term(), - adapt_to_sampling_rate( - coeff, self._sampling_rate, self.duration - ), + self._adapt_to_sampling_rate(coeff), ] ] # Build the interaction term for masked qubits qobj_list += [ [ make_interaction_term(masked=True), - adapt_to_sampling_rate( + self._adapt_to_sampling_rate( np.logical_not(coeff).astype(int), - self._sampling_rate, - self.duration, ), ] ] diff --git a/pulser-simulation/pulser_simulation/simulation.py b/pulser-simulation/pulser_simulation/simulation.py index 2193347e8..e05dbc15b 100644 --- a/pulser-simulation/pulser_simulation/simulation.py +++ b/pulser-simulation/pulser_simulation/simulation.py @@ -33,7 +33,7 @@ from pulser.result import SampledResult from pulser.sampler.samples import SequenceSamples from pulser.sequence._seq_drawer import draw_samples, draw_sequence -from pulser_simulation.hamiltonian import Hamiltonian, adapt_to_sampling_rate +from pulser_simulation.hamiltonian import Hamiltonian from pulser_simulation.qutip_result import QutipResult from pulser_simulation.simconfig import SimConfig from pulser_simulation.simresults import ( @@ -91,10 +91,7 @@ def __init__( device.validate_register(register) self._register = register # Check compatibility of samples and device: - if ( - sampled_seq._slm_mask.end > 0 - and not device.supports_slm_mask - ): + if sampled_seq._slm_mask.end > 0 and not device.supports_slm_mask: raise ValueError( "Samples use SLM mask but device does not have one." ) @@ -276,6 +273,7 @@ def set_evaluation_times( self._eval_times_array = np.union1d( eval_times, [0.0, self._tot_duration / 1000] ) + self._eval_times_instruction = value def get_hamiltonian(self, time: float) -> qutip.Qobj: r"""Get the Hamiltonian created from the sequence at a fixed time. From 0ffdb0fcaf82638f458f41004a5410007c8cc7bb Mon Sep 17 00:00:00 2001 From: a_corni Date: Tue, 7 Nov 2023 15:12:05 +0100 Subject: [PATCH 08/12] delete get_attr --- .../pulser_simulation/hamiltonian.py | 40 ++-- .../pulser_simulation/simulation.py | 158 ++++++++++++--- tests/test_simulation.py | 185 +++++++++++------- 3 files changed, 257 insertions(+), 126 deletions(-) diff --git a/pulser-simulation/pulser_simulation/hamiltonian.py b/pulser-simulation/pulser_simulation/hamiltonian.py index bc568b0f4..5abf4b84f 100644 --- a/pulser-simulation/pulser_simulation/hamiltonian.py +++ b/pulser-simulation/pulser_simulation/hamiltonian.py @@ -75,10 +75,10 @@ def __init__( self._qid_index = {qid: i for i, qid in enumerate(self._qdict)} # Compute sampling times - self.duration = self.samples_obj.max_duration + self._duration = self.samples_obj.max_duration self.sampling_times = self._adapt_to_sampling_rate( # Include extra time step for final instruction from samples: - np.arange(self.duration, dtype=np.double) + np.arange(self._duration, dtype=np.double) / 1000 ) @@ -95,7 +95,7 @@ def _adapt_to_sampling_rate(self, full_array: np.ndarray) -> np.ndarray: indices = np.linspace( 0, len(full_array) - 1, - int(self._sampling_rate * self.duration), + int(self._sampling_rate * self._duration), dtype=int, ) return cast(np.ndarray, full_array[indices]) @@ -105,19 +105,18 @@ def config(self) -> SimConfig: """The current configuration, as a SimConfig instance.""" return self._config - def _build_collapse_operators(self, prev_config: SimConfig) -> None: + def _build_collapse_operators(self, config: SimConfig) -> None: kraus_ops = [] self._collapse_ops = [] - if "dephasing" in self.config.noise: + if "dephasing" in config.noise: if self.basis_name == "digital" or self.basis_name == "all": # Go back to previous config - self.set_config(prev_config) raise NotImplementedError( "Cannot include dephasing noise in digital- or all-basis." ) # Probability of phase (Z) flip: # First order in prob - prob = self.config.dephasing_prob / 2 + prob = config.dephasing_prob / 2 n = self._size if prob > 0.1 and n > 1: warnings.warn( @@ -134,17 +133,16 @@ def _build_collapse_operators(self, prev_config: SimConfig) -> None: ] kraus_ops.append(k * qutip.sigmaz()) - if "depolarizing" in self.config.noise: + if "depolarizing" in config.noise: if self.basis_name == "digital" or self.basis_name == "all": # Go back to previous config - self.set_config(prev_config) raise NotImplementedError( "Cannot include depolarizing " + "noise in digital- or all-basis." ) # Probability of error occurrence - prob = self.config.depolarizing_prob / 4 + prob = config.depolarizing_prob / 4 n = self._size if prob > 0.1 and n > 1: warnings.warn( @@ -163,20 +161,19 @@ def _build_collapse_operators(self, prev_config: SimConfig) -> None: kraus_ops.append(k * qutip.sigmay()) kraus_ops.append(k * qutip.sigmaz()) - if "eff_noise" in self.config.noise: + if "eff_noise" in config.noise: if self.basis_name == "digital" or self.basis_name == "all": # Go back to previous config - self.set_config(prev_config) raise NotImplementedError( "Cannot include general " + "noise in digital- or all-basis." ) # Probability distribution of error occurences n = self._size - m = len(self.config.eff_noise_opers) + m = len(config.eff_noise_opers) if n > 1: for i in range(1, m): - prob_i = self.config.eff_noise_probs[i] + prob_i = config.eff_noise_probs[i] if prob_i > 0.1: warnings.warn( "The effective noise model is a first-order" @@ -186,16 +183,14 @@ def _build_collapse_operators(self, prev_config: SimConfig) -> None: ) break # Deriving Kraus operators - prob_id = self.config.eff_noise_probs[0] + prob_id = config.eff_noise_probs[0] self._collapse_ops += [ np.sqrt(prob_id**n) * qutip.tensor([self.op_matrix["I"] for _ in range(n)]) ] for i in range(1, m): - k = np.sqrt( - self.config.eff_noise_probs[i] * prob_id ** (n - 1) - ) - k_op = k * self.config.eff_noise_opers[i] + k = np.sqrt(config.eff_noise_probs[i] * prob_id ** (n - 1)) + k_op = k * config.eff_noise_opers[i] kraus_ops.append(k_op) # Building collapse operators @@ -221,7 +216,9 @@ def set_config(self, cfg: SimConfig) -> None: f"Interaction mode '{self._interaction}' does not support " f"simulation of noise types: {', '.join(not_supported)}." ) - prev_config = self.config if hasattr(self, "_config") else SimConfig() + if not hasattr(self, "basis_name"): + self._build_basis_and_op_matrices() + self._build_collapse_operators(cfg) self._config = cfg if not ("SPAM" in self.config.noise and self.config.eta > 0): self._bad_atoms = {qid: False for qid in self._qid_index} @@ -229,7 +226,6 @@ def set_config(self, cfg: SimConfig) -> None: self._doppler_detune = {qid: 0.0 for qid in self._qid_index} # Noise, samples and Hamiltonian update routine self._construct_hamiltonian() - self._build_collapse_operators(prev_config) def add_config(self, config: SimConfig) -> None: """Updates the current configuration with parameters of another one. @@ -604,7 +600,7 @@ def build_coeffs_ops(basis: str, addr: str) -> list[list]: ): # Build an array of binary coefficients for the interaction # term of unmasked qubits - coeff = np.ones(self.duration - 1) + coeff = np.ones(self._duration - 1) coeff[0 : self.samples_obj._slm_mask.end] = 0 # Build the interaction term for unmasked qubits qobj_list = [ diff --git a/pulser-simulation/pulser_simulation/simulation.py b/pulser-simulation/pulser_simulation/simulation.py index e05dbc15b..5af4c7a70 100644 --- a/pulser-simulation/pulser_simulation/simulation.py +++ b/pulser-simulation/pulser_simulation/simulation.py @@ -145,7 +145,7 @@ def __init__( "`sampling_rate` is too small, less than 4 data points." ) # Sets the config as well as builds the hamiltonian - self.hamiltonian = Hamiltonian( + self._hamiltonian = Hamiltonian( self.samples_obj, self._register.qubits, device, @@ -159,12 +159,97 @@ def __init__( if self.samples_obj._measurement: self._meas_basis = self.samples_obj._measurement else: - if self.basis_name in {"digital", "all"}: + if self._hamiltonian.basis_name in {"digital", "all"}: self._meas_basis = "digital" else: - self._meas_basis = self.basis_name + self._meas_basis = self._hamiltonian.basis_name self.set_initial_state("all-ground") + @property + def sampling_times(self) -> np.ndarray: + """The times at which hamiltonian is sampled.""" + return self._hamiltonian.sampling_times + + @property + def _sampling_rate(self) -> np.ndarray: + """The sampling rate.""" + return self._hamiltonian._sampling_rate + + @property + def dim(self) -> int: + """The dimension of the basis.""" + return self._hamiltonian.dim + + @property + def basis_name(self) -> str: + """The name of the basis.""" + return self._hamiltonian.basis_name + + @property + def basis(self) -> str: + """The basis in which result is expressed.""" + return self._hamiltonian.basis + + @property + def config(self) -> SimConfig: + """The current configuration, as a SimConfig instance.""" + return self._hamiltonian._config + + def set_config(self, cfg: SimConfig) -> None: + """Sets current config to cfg and updates simulation parameters. + + Args: + cfg: New configuration. + """ + if not isinstance(cfg, SimConfig): + raise ValueError(f"Object {cfg} is not a valid `SimConfig`.") + not_supported = ( + set(cfg.noise) + - cfg.supported_noises[self._hamiltonian._interaction] + ) + if not_supported: + raise NotImplementedError( + f"Interaction mode '{self._hamiltonian._interaction}' does not" + " support simulation of noise types:" + f"{', '.join(not_supported)}." + ) + self._hamiltonian.set_config(cfg) + + def add_config(self, config: SimConfig) -> None: + """Updates the current configuration with parameters of another one. + + Mostly useful when dealing with multiple noise types in different + configurations and wanting to merge these configurations together. + Adds simulation parameters to noises that weren't available in the + former SimConfig. Noises specified in both SimConfigs will keep + former noise parameters. + + Args: + config: SimConfig to retrieve parameters from. + """ + if not isinstance(config, SimConfig): + raise ValueError(f"Object {config} is not a valid `SimConfig`") + + not_supported = ( + set(config.noise) + - config.supported_noises[self._hamiltonian._interaction] + ) + if not_supported: + raise NotImplementedError( + f"Interaction mode '{self._hamiltonian._interaction}' does not" + " support simulation of noise types: " + f"{', '.join(not_supported)}." + ) + self._hamiltonian.add_config(config) + + def show_config(self, solver_options: bool = False) -> None: + """Shows current configuration.""" + print(self.config.__str__(solver_options)) + + def reset_config(self) -> None: + """Resets configuration to default.""" + self.set_config(SimConfig()) + @property def initial_state(self) -> qutip.Qobj: """The initial state of the simulation. @@ -187,15 +272,20 @@ def set_initial_state( if isinstance(state, str) and state == "all-ground": self._initial_state = qutip.tensor( [ - self.basis["u" if self._interaction == "XY" else "g"] - for _ in range(self._size) + self._hamiltonian.basis[ + "u" if self._hamiltonian._interaction == "XY" else "g" + ] + for _ in range(self._hamiltonian._size) ] ) else: state = cast(Union[np.ndarray, qutip.Qobj], state) shape = state.shape[0] - legal_shape = self.dim**self._size - legal_dims = [[self.dim] * self._size, [1] * self._size] + legal_shape = self._hamiltonian.dim**self._hamiltonian._size + legal_dims = [ + [self._hamiltonian.dim] * self._hamiltonian._size, + [1] * self._hamiltonian._size, + ] if shape != legal_shape: raise ValueError( "Incompatible shape of initial state." @@ -229,7 +319,7 @@ def set_evaluation_times( """Sets times at which the results of this simulation are returned.""" if isinstance(value, str): if value == "Full": - eval_times = np.copy(self.sampling_times) + eval_times = np.copy(self._hamiltonian.sampling_times) elif value == "Minimal": eval_times = np.array([]) else: @@ -245,12 +335,12 @@ def set_evaluation_times( ) indices = np.linspace( 0, - len(self.sampling_times) - 1, - int(value * len(self.sampling_times)), + len(self._hamiltonian.sampling_times) - 1, + int(value * len(self._hamiltonian.sampling_times)), dtype=int, ) # Note: if `value` is very small `eval_times` is an empty list: - eval_times = self.sampling_times[indices] + eval_times = self._hamiltonian.sampling_times[indices] elif isinstance(value, (list, tuple, np.ndarray)): if np.max(value, initial=0) > self._tot_duration / 1000: raise ValueError( @@ -302,7 +392,9 @@ def get_hamiltonian(self, time: float) -> qutip.Qobj: f"Provided time (`time` = {time}) must be " "greater than or equal to 0." ) - return self._hamiltonian(time / 1000) # Creates new Qutip.Qobj + return self._hamiltonian._hamiltonian( + time / 1000 + ) # Creates new Qutip.Qobj # Run Simulation Evolution using Qutip def run( @@ -349,7 +441,10 @@ def run( for k in ("epsilon", "epsilon_prime") } if self.config.eta > 0 and self.initial_state != qutip.tensor( - [self.basis["g"] for _ in range(self._size)] + [ + self._hamiltonian.basis["g"] + for _ in range(self._hamiltonian._size) + ] ): raise NotImplementedError( "Can't combine state preparation errors with an initial " @@ -373,16 +468,16 @@ def _run_solver() -> CoherentResults: or "eff_noise" in self.config.noise ): result = qutip.mesolve( - self._hamiltonian, + self._hamiltonian._hamiltonian, self.initial_state, self._eval_times_array, - self._collapse_ops, + self._hamiltonian._collapse_ops, progress_bar=p_bar, options=solv_ops, ) else: result = qutip.sesolve( - self._hamiltonian, + self._hamiltonian._hamiltonian, self.initial_state, self._eval_times_array, progress_bar=p_bar, @@ -390,17 +485,17 @@ def _run_solver() -> CoherentResults: ) results = [ QutipResult( - tuple(self._qdict), + tuple(self._hamiltonian._qdict), self._meas_basis, state, - self._meas_basis == self.basis_name, + self._meas_basis == self._hamiltonian.basis_name, ) for state in result.states ] return CoherentResults( results, - self._size, - self.basis_name, + self._hamiltonian._size, + self._hamiltonian.basis_name, self._eval_times_array, self._meas_basis, meas_errors, @@ -419,7 +514,9 @@ def _run_solver() -> CoherentResults: initial_configs = Counter( "".join( ( - np.random.uniform(size=len(self._qid_index)) + np.random.uniform( + size=len(self._hamiltonian._qid_index) + ) < self.config.eta ) .astype(int) @@ -441,16 +538,16 @@ def _run_solver() -> CoherentResults: if not update_ham: initial_state, reps = initial_configs[i] # We load the initial state manually - self._bad_atoms = dict( + self._hamiltonian._bad_atoms = dict( zip( - self._qid_index, + self._hamiltonian._qid_index, np.array(list(initial_state)).astype(bool), ) ) else: reps = 1 # At each run, new random noise: new Hamiltonian - self._construct_hamiltonian(update=update_ham) + self._hamiltonian._construct_hamiltonian(update=update_ham) # Get CoherentResults instance from sequence with added noise: cleanres_noisyseq = _run_solver() # Extract statistics at eval time: @@ -464,13 +561,17 @@ def _run_solver() -> CoherentResults: ) n_measures = self.config.runs * self.config.samples_per_run results = [ - SampledResult(tuple(self._qdict), self._meas_basis, total_count[t]) + SampledResult( + tuple(self._hamiltonian._qdict), + self._meas_basis, + total_count[t], + ) for t in time_indices ] return NoisyResults( results, - self._size, - self.basis_name, + self._hamiltonian._size, + self._hamiltonian.basis_name, self._eval_times_array, n_measures, ) @@ -585,9 +686,6 @@ def from_sequence( evaluation_times, ) - def __getattr__(self, name: str) -> Any: - return getattr(self.hamiltonian, name) - class Simulation: r"""Simulation of a pulse sequence using QuTiP. diff --git a/tests/test_simulation.py b/tests/test_simulation.py index f947ee43a..c41b19a1b 100644 --- a/tests/test_simulation.py +++ b/tests/test_simulation.py @@ -147,10 +147,14 @@ def test_initialization_and_construction_of_hamiltonian(seq, mod_device): for ch in sampled_seq.channels ] ) - assert sim._qdict == seq.qubit_info - assert sim._size == len(seq.qubit_info) + assert sim._hamiltonian._qdict == seq.qubit_info + assert sim._hamiltonian._size == len(seq.qubit_info) assert sim._tot_duration == 9000 # seq has 9 pulses of 1µs - assert sim._qid_index == {"control1": 0, "target": 1, "control2": 2} + assert sim._hamiltonian._qid_index == { + "control1": 0, + "target": 1, + "control2": 2, + } with pytest.raises(ValueError, match="too small, less than"): QutipEmulator.from_sequence(seq, sampling_rate=0.0001) @@ -164,12 +168,16 @@ def test_initialization_and_construction_of_hamiltonian(seq, mod_device): sim._sampling_rate * sim._tot_duration ) - assert isinstance(sim._hamiltonian, qutip.QobjEvo) + assert isinstance(sim._hamiltonian._hamiltonian, qutip.QobjEvo) # Checks adapt() method: - assert bool(set(sim._hamiltonian.tlist).intersection(sim.sampling_times)) - for qobjevo in sim._hamiltonian.ops: + assert bool( + set(sim._hamiltonian._hamiltonian.tlist).intersection( + sim.sampling_times + ) + ) + for qobjevo in sim._hamiltonian._hamiltonian.ops: for sh in qobjevo.qobj.shape: - assert sh == sim.dim**sim._size + assert sh == sim.dim**sim._hamiltonian._size assert not seq.is_parametrized() with pytest.warns(UserWarning, match="returns a copy of itself"): @@ -197,7 +205,7 @@ def test_extraction_of_sequences(seq): if addr == "Global": for slot in seq._schedule[channel]: if isinstance(slot.type, Pulse): - samples = sim.samples[addr][basis] + samples = sim._hamiltonian.samples[addr][basis] assert ( samples["amp"][slot.ti : slot.tf] == slot.type.amplitude.samples @@ -214,7 +222,7 @@ def test_extraction_of_sequences(seq): for slot in seq._schedule[channel]: if isinstance(slot.type, Pulse): for qubit in slot.targets: # TO DO: multiaddressing?? - samples = sim.samples[addr][basis][qubit] + samples = sim._hamiltonian.samples[addr][basis][qubit] assert ( samples["amp"][slot.ti : slot.tf] == slot.type.amplitude.samples @@ -240,29 +248,29 @@ def test_building_basis_and_projection_operators(seq, reg): "h": qutip.basis(3, 2), } assert ( - sim.op_matrix["sigma_rr"] + sim._hamiltonian.op_matrix["sigma_rr"] == qutip.basis(3, 0) * qutip.basis(3, 0).dag() ) assert ( - sim.op_matrix["sigma_gr"] + sim._hamiltonian.op_matrix["sigma_gr"] == qutip.basis(3, 1) * qutip.basis(3, 0).dag() ) assert ( - sim.op_matrix["sigma_hg"] + sim._hamiltonian.op_matrix["sigma_hg"] == qutip.basis(3, 2) * qutip.basis(3, 1).dag() ) # Check local operator building method: with pytest.raises(ValueError, match="Duplicate atom"): - sim.build_operator([("sigma_gg", ["target", "target"])]) + sim._hamiltonian.build_operator([("sigma_gg", ["target", "target"])]) with pytest.raises(ValueError, match="not a valid operator"): - sim.build_operator([("wrong", ["target"])]) + sim._hamiltonian.build_operator([("wrong", ["target"])]) with pytest.raises(ValueError, match="Invalid qubit names: {'wrong'}"): - sim.build_operator([("sigma_gg", ["wrong"])]) + sim._hamiltonian.build_operator([("sigma_gg", ["wrong"])]) # Check building operator with one operator - op_standard = sim.build_operator([("sigma_gg", ["target"])]) - op_one = sim.build_operator(("sigma_gg", ["target"])) + op_standard = sim._hamiltonian.build_operator([("sigma_gg", ["target"])]) + op_one = sim._hamiltonian.build_operator(("sigma_gg", ["target"])) assert np.linalg.norm(op_standard - op_one) < 1e-10 # Global ground-rydberg @@ -275,11 +283,11 @@ def test_building_basis_and_projection_operators(seq, reg): assert sim2.dim == 2 assert sim2.basis == {"r": qutip.basis(2, 0), "g": qutip.basis(2, 1)} assert ( - sim2.op_matrix["sigma_rr"] + sim2._hamiltonian.op_matrix["sigma_rr"] == qutip.basis(2, 0) * qutip.basis(2, 0).dag() ) assert ( - sim2.op_matrix["sigma_gr"] + sim2._hamiltonian.op_matrix["sigma_gr"] == qutip.basis(2, 1) * qutip.basis(2, 0).dag() ) @@ -292,11 +300,11 @@ def test_building_basis_and_projection_operators(seq, reg): assert sim2b.dim == 2 assert sim2b.basis == {"g": qutip.basis(2, 0), "h": qutip.basis(2, 1)} assert ( - sim2b.op_matrix["sigma_gg"] + sim2b._hamiltonian.op_matrix["sigma_gg"] == qutip.basis(2, 0) * qutip.basis(2, 0).dag() ) assert ( - sim2b.op_matrix["sigma_hg"] + sim2b._hamiltonian.op_matrix["sigma_hg"] == qutip.basis(2, 1) * qutip.basis(2, 0).dag() ) @@ -309,11 +317,11 @@ def test_building_basis_and_projection_operators(seq, reg): assert sim2c.dim == 2 assert sim2c.basis == {"r": qutip.basis(2, 0), "g": qutip.basis(2, 1)} assert ( - sim2c.op_matrix["sigma_rr"] + sim2c._hamiltonian.op_matrix["sigma_rr"] == qutip.basis(2, 0) * qutip.basis(2, 0).dag() ) assert ( - sim2c.op_matrix["sigma_gr"] + sim2c._hamiltonian.op_matrix["sigma_gr"] == qutip.basis(2, 1) * qutip.basis(2, 0).dag() ) @@ -332,15 +340,15 @@ def test_building_basis_and_projection_operators(seq, reg): assert sim2.dim == 2 assert sim2.basis == {"u": qutip.basis(2, 0), "d": qutip.basis(2, 1)} assert ( - sim2.op_matrix["sigma_uu"] + sim2._hamiltonian.op_matrix["sigma_uu"] == qutip.basis(2, 0) * qutip.basis(2, 0).dag() ) assert ( - sim2.op_matrix["sigma_du"] + sim2._hamiltonian.op_matrix["sigma_du"] == qutip.basis(2, 1) * qutip.basis(2, 0).dag() ) assert ( - sim2.op_matrix["sigma_ud"] + sim2._hamiltonian.op_matrix["sigma_ud"] == qutip.basis(2, 0) * qutip.basis(2, 1).dag() ) @@ -363,10 +371,12 @@ def test_empty_sequences(reg): seq.delay(100, "test") emu = QutipEmulator.from_sequence(seq, config=SimConfig(noise="SPAM")) - assert not emu.samples["Global"] - for basis in emu.samples["Local"]: - for q in emu.samples["Local"][basis]: - for qty_values in emu.samples["Local"][basis][q].values(): + assert not emu._hamiltonian.samples["Global"] + for basis in emu._hamiltonian.samples["Local"]: + for q in emu._hamiltonian.samples["Local"][basis]: + for qty_values in emu._hamiltonian.samples["Local"][basis][ + q + ].values(): np.testing.assert_equal(qty_values, 0) @@ -431,10 +441,10 @@ def test_single_atom_simulation(): ) one_sim = QutipEmulator.from_sequence(one_seq) one_res = one_sim.run() - assert one_res._size == one_sim._size + assert one_res._size == one_sim._hamiltonian._size one_sim = QutipEmulator.from_sequence(one_seq, evaluation_times="Minimal") one_resb = one_sim.run() - assert one_resb._size == one_sim._size + assert one_resb._size == one_sim._hamiltonian._size def test_add_max_step_and_delays(): @@ -463,11 +473,15 @@ def test_run(seq, patch_plt_show): with patch("matplotlib.pyplot.savefig"): sim.draw(draw_phase_area=True, fig_name="my_fig.pdf") bad_initial = np.array([1.0]) - good_initial_array = np.r_[1, np.zeros(sim.dim**sim._size - 1)] + good_initial_array = np.r_[ + 1, np.zeros(sim.dim**sim._hamiltonian._size - 1) + ] good_initial_qobj = qutip.tensor( - [qutip.basis(sim.dim, 0) for _ in range(sim._size)] + [qutip.basis(sim.dim, 0) for _ in range(sim._hamiltonian._size)] + ) + good_initial_qobj_no_dims = qutip.basis( + sim.dim**sim._hamiltonian._size, 2 ) - good_initial_qobj_no_dims = qutip.basis(sim.dim**sim._size, 2) with pytest.raises( ValueError, match="Incompatible shape of initial state" @@ -663,12 +677,12 @@ def test_config(): def test_noise(seq, matrices): - np.random.seed(2) + np.random.seed(3) sim2 = QutipEmulator.from_sequence( seq, sampling_rate=0.01, config=SimConfig(noise=("SPAM"), eta=0.9) ) assert sim2.run().sample_final_state() == Counter( - {"000": 955, "001": 6, "100": 39} + {"000": 857, "110": 73, "100": 70} ) with pytest.raises(NotImplementedError, match="Cannot include"): sim2.set_config(SimConfig(noise="dephasing")) @@ -687,14 +701,16 @@ def test_noise(seq, matrices): "epsilon": 0.01, "epsilon_prime": 0.05, } - assert sim2.samples["Global"] == {} - assert any(sim2._bad_atoms.values()) + assert sim2._hamiltonian.samples["Global"] == {} + assert any(sim2._hamiltonian._bad_atoms.values()) for basis in ("ground-rydberg", "digital"): - for t in sim2._bad_atoms: - if not sim2._bad_atoms[t]: + for t in sim2._hamiltonian._bad_atoms: + if not sim2._hamiltonian._bad_atoms[t]: continue for qty in ("amp", "det", "phase"): - assert np.all(sim2.samples["Local"][basis][t][qty] == 0.0) + assert np.all( + sim2._hamiltonian.samples["Local"][basis][t][qty] == 0.0 + ) def test_noise_with_zero_epsilons(seq, matrices): @@ -729,7 +745,7 @@ def test_dephasing(): seq, sampling_rate=0.01, config=SimConfig(noise="dephasing") ) assert sim.run().sample_final_state() == Counter({"0": 595, "1": 405}) - assert len(sim._collapse_ops) != 0 + assert len(sim._hamiltonian._collapse_ops) != 0 with pytest.warns(UserWarning, match="first-order"): reg = Register.from_coordinates([(0, 0), (0, 10)], prefix="q") seq2 = Sequence(reg, Chadoq2) @@ -756,7 +772,7 @@ def test_depolarizing(): assert sim.run().sample_final_state() == Counter({"0": 587, "1": 413}) 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 + assert len(sim._hamiltonian._collapse_ops) != 0 with pytest.warns(UserWarning, match="first-order"): reg = Register.from_coordinates([(0, 0), (0, 10)], prefix="q") seq2 = Sequence(reg, Chadoq2) @@ -790,10 +806,10 @@ def test_eff_noise(matrices): seq, sampling_rate=0.01, config=SimConfig(noise="dephasing") ) assert ( - sim._collapse_ops == sim_dph._collapse_ops + sim._hamiltonian._collapse_ops == sim_dph._hamiltonian._collapse_ops and sim.run().states[-1] == sim_dph.run().states[-1] ) - assert len(sim._collapse_ops) != 0 + assert len(sim._hamiltonian._collapse_ops) != 0 with pytest.warns(UserWarning, match="first-order"): reg = Register.from_coordinates([(0, 0), (0, 10)], prefix="q") seq2 = Sequence(reg, Chadoq2) @@ -933,9 +949,11 @@ def test_run_xy(): sim = QutipEmulator.from_sequence(simple_seq, sampling_rate=0.01) - good_initial_array = np.r_[1, np.zeros(sim.dim**sim._size - 1)] + good_initial_array = np.r_[ + 1, np.zeros(sim.dim**sim._hamiltonian._size - 1) + ] good_initial_qobj = qutip.tensor( - [qutip.basis(sim.dim, 0) for _ in range(sim._size)] + [qutip.basis(sim.dim, 0) for _ in range(sim._hamiltonian._size)] ) sim.set_initial_state(good_initial_array) assert sim.initial_state == good_initial_qobj @@ -967,7 +985,7 @@ def test_noisy_xy(): sim.set_config(SimConfig(("SPAM", "doppler"))) sim.set_config(SimConfig("SPAM", eta=0.4)) - assert sim._bad_atoms == { + assert sim._hamiltonian._bad_atoms == { "atom0": True, "atom1": False, "atom2": True, @@ -1108,44 +1126,51 @@ def test_mask_local_channel(): assert seq_._slm_mask_targets == {"q0", "q3"} sim = QutipEmulator.from_sequence(seq_) assert np.array_equal( - sim.samples["Global"]["ground-rydberg"]["amp"], + sim._hamiltonian.samples["Global"]["ground-rydberg"]["amp"], np.concatenate((pulse.amplitude.samples, [0])), ) assert np.array_equal( - sim.samples["Global"]["ground-rydberg"]["det"], + sim._hamiltonian.samples["Global"]["ground-rydberg"]["det"], np.concatenate((pulse.detuning.samples, [0])), ) - assert np.all(sim.samples["Global"]["ground-rydberg"]["phase"] == 0.0) + assert np.all( + sim._hamiltonian.samples["Global"]["ground-rydberg"]["phase"] == 0.0 + ) qubits = ["q0", "q1", "q2", "q3"] masked_qubits = ["q0", "q3"] for q in qubits: if q in masked_qubits: assert np.array_equal( - sim.samples["Local"]["ground-rydberg"][q]["det"], + sim._hamiltonian.samples["Local"]["ground-rydberg"][q]["det"], np.concatenate( (-10 / len(masked_qubits) * pulse.amplitude.samples, [0]) ), ) else: assert np.all( - sim.samples["Local"]["ground-rydberg"][q]["det"] == 0.0 + sim._hamiltonian.samples["Local"]["ground-rydberg"][q]["det"] + == 0.0 ) - assert np.all(sim.samples["Local"]["ground-rydberg"][q]["amp"] == 0.0) assert np.all( - sim.samples["Local"]["ground-rydberg"][q]["phase"] == 0.0 + sim._hamiltonian.samples["Local"]["ground-rydberg"][q]["amp"] + == 0.0 + ) + assert np.all( + sim._hamiltonian.samples["Local"]["ground-rydberg"][q]["phase"] + == 0.0 ) assert np.array_equal( - sim.samples["Local"]["digital"]["q0"]["amp"], + sim._hamiltonian.samples["Local"]["digital"]["q0"]["amp"], np.concatenate((pulse2.amplitude.samples, [0])), ) assert np.array_equal( - sim.samples["Local"]["digital"]["q0"]["det"], + sim._hamiltonian.samples["Local"]["digital"]["q0"]["det"], np.concatenate((pulse2.detuning.samples, [0])), ) assert np.all( np.isclose( - sim.samples["Local"]["digital"]["q0"]["phase"], + sim._hamiltonian.samples["Local"]["digital"]["q0"]["phase"], np.concatenate((np.pi * np.ones(1000), [0])), ) ) @@ -1163,13 +1188,13 @@ def test_effective_size_intersection(): sim = QutipEmulator.from_sequence(seq, sampling_rate=0.01) sim.set_config(SimConfig("SPAM", eta=0.4)) - assert sim._bad_atoms == { + assert sim._hamiltonian._bad_atoms == { "atom0": True, "atom1": False, "atom2": True, "atom3": False, } - assert sim.get_hamiltonian(0) != 0 * sim.build_operator( + assert sim.get_hamiltonian(0) != 0 * sim._hamiltonian.build_operator( [("I", "global")] ) @@ -1194,14 +1219,16 @@ def test_effective_size_disjoint(channel_type): assert seq._slm_mask_time == [0, 1500] sim = QutipEmulator.from_sequence(seq, sampling_rate=0.01) sim.set_config(SimConfig("SPAM", eta=0.4)) - assert sim._bad_atoms == { + assert sim._hamiltonian._bad_atoms == { "atom0": True, "atom1": False, "atom2": True, "atom3": False, } if channel_type == "mw_global": - assert sim.get_hamiltonian(0) == 0.5 * amp * sim.build_operator( + assert sim.get_hamiltonian( + 0 + ) == 0.5 * amp * sim._hamiltonian.build_operator( [(qutip.sigmax(), ["atom3"])] ) else: @@ -1209,23 +1236,29 @@ def test_effective_size_disjoint(channel_type): "ground-rydberg" if channel_type == "rydberg_global" else "digital" ) assert np.array_equal( - sim.samples["Local"][basis]["atom1"]["amp"], + sim._hamiltonian.samples["Local"][basis]["atom1"]["amp"], np.concatenate((rise.amplitude.samples, [0])), ) assert np.array_equal( - sim.samples["Local"][basis]["atom3"]["amp"], + sim._hamiltonian.samples["Local"][basis]["atom3"]["amp"], np.concatenate((rise.amplitude.samples, [0])), ) # SLM assert np.all( - sim.samples["Local"]["ground-rydberg"]["atom1"]["det"] + sim._hamiltonian.samples["Local"]["ground-rydberg"]["atom1"]["det"] == -10 * np.concatenate((rise.amplitude.samples, [0])) ) if channel_type == "raman_global": - assert np.all(sim.samples["Local"][basis]["atom1"]["det"] == 0.0) - assert np.all(sim.samples["Local"][basis]["atom3"]["det"] == 0.0) + assert np.all( + sim._hamiltonian.samples["Local"][basis]["atom1"]["det"] == 0.0 + ) + assert np.all( + sim._hamiltonian.samples["Local"][basis]["atom3"]["det"] == 0.0 + ) for q in ["atom1", "atom3"]: - assert np.all(sim.samples["Local"][basis][q]["phase"] == 0.0) + assert np.all( + sim._hamiltonian.samples["Local"][basis][q]["phase"] == 0.0 + ) def test_simulation_with_modulation(mod_device, reg, patch_plt_show): @@ -1259,8 +1292,10 @@ def test_simulation_with_modulation(mod_device, reg, patch_plt_show): seq, with_modulation=True, config=sim_config ) - assert sim.samples["Global"] == {} # All samples stored in local - raman_samples = sim.samples["Local"]["digital"] + assert ( + sim._hamiltonian.samples["Global"] == {} + ) # All samples stored in local + raman_samples = sim._hamiltonian.samples["Local"]["digital"] # Local pulses for qid, time_slice in [ ("target", slice(0, mod_dt)), @@ -1272,7 +1307,8 @@ def test_simulation_with_modulation(mod_device, reg, patch_plt_show): atol=1e-2, ) np.testing.assert_equal( - raman_samples[qid]["det"][time_slice], sim._doppler_detune[qid] + raman_samples[qid]["det"][time_slice], + sim._hamiltonian._doppler_detune[qid], ) np.testing.assert_allclose( raman_samples[qid]["phase"][time_slice], pulse1.phase @@ -1285,7 +1321,7 @@ def pos_factor(qid): # Global pulse time_slice = slice(2 * mod_dt, 3 * mod_dt) - rydberg_samples = sim.samples["Local"]["ground-rydberg"] + rydberg_samples = sim._hamiltonian.samples["Local"]["ground-rydberg"] noise_amp_base = rydberg_samples["target"]["amp"][time_slice] / ( pulse1_mod_samples * pos_factor("target") ) @@ -1295,7 +1331,8 @@ def pos_factor(qid): pulse1_mod_samples * noise_amp_base * pos_factor(qid), ) np.testing.assert_equal( - rydberg_samples[qid]["det"][time_slice], sim._doppler_detune[qid] + rydberg_samples[qid]["det"][time_slice], + sim._hamiltonian._doppler_detune[qid], ) np.testing.assert_allclose( rydberg_samples[qid]["phase"][time_slice], pulse1.phase From 5466073eb26ad00f32086a24fc5f94251c94cbb2 Mon Sep 17 00:00:00 2001 From: a_corni Date: Tue, 7 Nov 2023 18:36:27 +0100 Subject: [PATCH 09/12] using NoiseModel for config in Hamiltonian --- .../pulser_simulation/hamiltonian.py | 98 ++++++++++--------- .../pulser_simulation/simconfig.py | 35 ++++--- .../pulser_simulation/simulation.py | 44 +++++---- tests/test_sequence_sampler.py | 4 +- tests/test_simulation.py | 38 ++++--- 5 files changed, 125 insertions(+), 94 deletions(-) diff --git a/pulser-simulation/pulser_simulation/hamiltonian.py b/pulser-simulation/pulser_simulation/hamiltonian.py index 5abf4b84f..81e6a3e8c 100644 --- a/pulser-simulation/pulser_simulation/hamiltonian.py +++ b/pulser-simulation/pulser_simulation/hamiltonian.py @@ -25,10 +25,11 @@ import numpy as np import qutip +from pulser.backend.noise_model import NoiseModel from pulser.devices._device_datacls import BaseDevice from pulser.register.base_register import QubitId from pulser.sampler.samples import SequenceSamples, _PulseTargetSlot -from pulser_simulation.simconfig import SimConfig +from pulser_simulation.simconfig import SUPPORTED_NOISES, doppler_sigma class Hamiltonian: @@ -50,7 +51,7 @@ def __init__( qdict: dict[QubitId, np.ndarray], device: BaseDevice, sampling_rate: float, - config: SimConfig, + config: NoiseModel, ) -> None: """Instantiates a Hamiltonian object.""" self.samples_obj = samples_obj @@ -60,7 +61,7 @@ def __init__( # Type hints for attributes defined outside of __init__ self.basis_name: str - self._config: SimConfig + self._config: NoiseModel self.op_matrix: dict[str, qutip.Qobj] self.basis: dict[str, qutip.Qobj] self.dim: int @@ -101,14 +102,14 @@ def _adapt_to_sampling_rate(self, full_array: np.ndarray) -> np.ndarray: return cast(np.ndarray, full_array[indices]) @property - def config(self) -> SimConfig: - """The current configuration, as a SimConfig instance.""" + def config(self) -> NoiseModel: + """The current configuration, as a NoiseModel instance.""" return self._config - def _build_collapse_operators(self, config: SimConfig) -> None: + def _build_collapse_operators(self, config: NoiseModel) -> None: kraus_ops = [] self._collapse_ops = [] - if "dephasing" in config.noise: + if "dephasing" in config.noise_types: if self.basis_name == "digital" or self.basis_name == "all": # Go back to previous config raise NotImplementedError( @@ -133,7 +134,7 @@ def _build_collapse_operators(self, config: SimConfig) -> None: ] kraus_ops.append(k * qutip.sigmaz()) - if "depolarizing" in config.noise: + if "depolarizing" in config.noise_types: if self.basis_name == "digital" or self.basis_name == "all": # Go back to previous config raise NotImplementedError( @@ -161,7 +162,7 @@ def _build_collapse_operators(self, config: SimConfig) -> None: kraus_ops.append(k * qutip.sigmay()) kraus_ops.append(k * qutip.sigmaz()) - if "eff_noise" in config.noise: + if "eff_noise" in config.noise_types: if self.basis_name == "digital" or self.basis_name == "all": # Go back to previous config raise NotImplementedError( @@ -200,16 +201,16 @@ def _build_collapse_operators(self, config: SimConfig) -> None: for qid in self._qid_index ] - def set_config(self, cfg: SimConfig) -> None: + def set_config(self, cfg: NoiseModel) -> None: """Sets current config to cfg and updates simulation parameters. Args: cfg: New configuration. """ - if not isinstance(cfg, SimConfig): - raise ValueError(f"Object {cfg} is not a valid `SimConfig`.") + if not isinstance(cfg, NoiseModel): + raise ValueError(f"Object {cfg} is not a valid `NoiseModel`.") not_supported = ( - set(cfg.noise) - cfg.supported_noises[self._interaction] + set(cfg.noise_types) - SUPPORTED_NOISES[self._interaction] ) if not_supported: raise NotImplementedError( @@ -220,30 +221,33 @@ def set_config(self, cfg: SimConfig) -> None: self._build_basis_and_op_matrices() self._build_collapse_operators(cfg) self._config = cfg - if not ("SPAM" in self.config.noise and self.config.eta > 0): + if not ( + "SPAM" in self.config.noise_types + and self.config.state_prep_error > 0 + ): self._bad_atoms = {qid: False for qid in self._qid_index} - if "doppler" not in self.config.noise: + if "doppler" not in self.config.noise_types: self._doppler_detune = {qid: 0.0 for qid in self._qid_index} # Noise, samples and Hamiltonian update routine self._construct_hamiltonian() - def add_config(self, config: SimConfig) -> None: + def add_config(self, config: NoiseModel) -> None: """Updates the current configuration with parameters of another one. Mostly useful when dealing with multiple noise types in different configurations and wanting to merge these configurations together. Adds simulation parameters to noises that weren't available in the - former SimConfig. Noises specified in both SimConfigs will keep + former NoiseModel. Noises specified in both NoiseModels will keep former noise parameters. Args: - config: SimConfig to retrieve parameters from. + config: NoiseModel to retrieve parameters from. """ - if not isinstance(config, SimConfig): - raise ValueError(f"Object {config} is not a valid `SimConfig`") + if not isinstance(config, NoiseModel): + raise ValueError(f"Object {config} is not a valid `NoiseModel`") not_supported = ( - set(config.noise) - config.supported_noises[self._interaction] + set(config.noise_types) - SUPPORTED_NOISES[self._interaction] ) if not_supported: raise NotImplementedError( @@ -251,17 +255,18 @@ def add_config(self, config: SimConfig) -> None: f"simulation of noise types: {', '.join(not_supported)}." ) - old_noise_set = set(self.config.noise) - new_noise_set = old_noise_set.union(config.noise) + old_noise_set = set(self.config.noise_types) + new_noise_set = old_noise_set.union(config.noise_types) diff_noise_set = new_noise_set - old_noise_set + print(diff_noise_set) # Create temporary param_dict to add noise parameters: param_dict: dict[str, Any] = asdict(self._config) # Begin populating with added noise parameters: - param_dict["noise"] = tuple(new_noise_set) + param_dict["noise_types"] = tuple(new_noise_set) if "SPAM" in diff_noise_set: - param_dict["eta"] = config.eta - param_dict["epsilon"] = config.epsilon - param_dict["epsilon_prime"] = config.epsilon_prime + param_dict["state_prep_error"] = config.state_prep_error + param_dict["p_false_pos"] = config.p_false_pos + param_dict["p_false_neg"] = config.p_false_neg if "doppler" in diff_noise_set: param_dict["temperature"] = config.temperature if "amplitude" in diff_noise_set: @@ -273,29 +278,23 @@ def add_config(self, config: SimConfig) -> None: if "eff_noise" in diff_noise_set: param_dict["eff_noise_opers"] = config.eff_noise_opers param_dict["eff_noise_probs"] = config.eff_noise_probs - param_dict["temperature"] *= 1.0e6 # update runs: param_dict["runs"] = config.runs param_dict["samples_per_run"] = config.samples_per_run # set config with the new parameters: - self.set_config(SimConfig(**param_dict)) - - def show_config(self, solver_options: bool = False) -> None: - """Shows current configuration.""" - print(self._config.__str__(solver_options)) - - def reset_config(self) -> None: - """Resets configuration to default.""" - self.set_config(SimConfig()) + self.set_config(NoiseModel(**param_dict)) def _extract_samples(self) -> None: """Populates samples dictionary with every pulse in the sequence.""" local_noises = True - if set(self.config.noise).issubset( + if set(self.config.noise_types).issubset( {"dephasing", "SPAM", "depolarizing", "eff_noise"} ): - local_noises = "SPAM" in self.config.noise and self.config.eta > 0 + local_noises = ( + "SPAM" in self.config.noise_types + and self.config.state_prep_error > 0 + ) samples = self.samples_obj.to_nested_dict(all_local=local_noises) def add_noise( @@ -312,12 +311,12 @@ def add_noise( 0, np.random.normal(1.0, self.config.amp_sigma) ) for qid in slot.targets: - if "doppler" in self.config.noise: + if "doppler" in self.config.noise_types: noise_det = self._doppler_detune[qid] samples_dict[qid]["det"][slot.ti : slot.tf] += noise_det # Gaussian beam loss in amplitude for global pulses only # Noise is drawn at random for each pulse - if "amplitude" in self.config.noise and is_global_pulse: + if "amplitude" in self.config.noise_types and is_global_pulse: position = self._qdict[qid] r = np.linalg.norm(position) w0 = self.config.laser_waist @@ -389,7 +388,7 @@ def build_operator(self, operations: Union[list, tuple]) -> qutip.Qobj: for qubit in qubits: k = self._qid_index[qubit] op_list[k] = operator - return qutip.tensor(op_list) + return qutip.tensor(list(map(qutip.Qobj, op_list))) def _update_noise(self) -> None: """Updates noise random parameters. @@ -397,15 +396,20 @@ def _update_noise(self) -> None: Used at the start of each run. If SPAM isn't in chosen noises, all atoms are set to be correctly prepared. """ - if "SPAM" in self.config.noise and self.config.eta > 0: + if ( + "SPAM" in self.config.noise_types + and self.config.state_prep_error > 0 + ): dist = ( np.random.uniform(size=len(self._qid_index)) - < self.config.spam_dict["eta"] + < self.config.state_prep_error ) self._bad_atoms = dict(zip(self._qid_index, dist)) - if "doppler" in self.config.noise: + if "doppler" in self.config.noise_types: detune = np.random.normal( - 0, self.config.doppler_sigma, size=len(self._qid_index) + 0, + doppler_sigma(self.config.temperature * 1e-6), + size=len(self._qid_index), ) self._doppler_detune = dict(zip(self._qid_index, detune)) @@ -453,8 +457,6 @@ def _construct_hamiltonian(self, update: bool = True) -> None: if update: self._update_noise() self._extract_samples() - if not hasattr(self, "basis_name"): - self._build_basis_and_op_matrices() def make_vdw_term(q1: QubitId, q2: QubitId) -> qutip.Qobj: """Construct the Van der Waals interaction Term. diff --git a/pulser-simulation/pulser_simulation/simconfig.py b/pulser-simulation/pulser_simulation/simconfig.py index 54f469aa0..b186332a7 100644 --- a/pulser-simulation/pulser_simulation/simconfig.py +++ b/pulser-simulation/pulser_simulation/simconfig.py @@ -32,6 +32,27 @@ T = TypeVar("T", bound="SimConfig") +SUPPORTED_NOISES: dict = { + "ising": { + "dephasing", + "doppler", + "amplitude", + "SPAM", + "depolarizing", + "eff_noise", + }, + "XY": {"SPAM"}, +} + + +def doppler_sigma(temperature: float) -> float: + """Standard deviation for Doppler shifting due to thermal motion. + + Arg: + temperature: The temperature in K. + """ + return KEFF * sqrt(KB * temperature / MASS) + @dataclass(frozen=True) class SimConfig: @@ -160,7 +181,7 @@ def spam_dict(self) -> dict[str, float]: @property def doppler_sigma(self) -> float: """Standard deviation for Doppler shifting due to thermal motion.""" - return KEFF * sqrt(KB * self.temperature / MASS) + return doppler_sigma(self.temperature) def __str__(self, solver_options: bool = False) -> str: lines = [ @@ -220,14 +241,4 @@ def _check_eff_noise_opers_type(self) -> None: @property def supported_noises(self) -> dict: """Return the noises implemented on pulser.""" - return { - "ising": { - "dephasing", - "doppler", - "amplitude", - "SPAM", - "depolarizing", - "eff_noise", - }, - "XY": {"SPAM"}, - } + return SUPPORTED_NOISES diff --git a/pulser-simulation/pulser_simulation/simulation.py b/pulser-simulation/pulser_simulation/simulation.py index 5af4c7a70..ce8ca2a6a 100644 --- a/pulser-simulation/pulser_simulation/simulation.py +++ b/pulser-simulation/pulser_simulation/simulation.py @@ -28,6 +28,7 @@ import pulser.sampler as sampler from pulser import Sequence +from pulser.backend.noise_model import NoiseModel from pulser.devices._device_datacls import BaseDevice from pulser.register.base_register import BaseRegister from pulser.result import SampledResult @@ -150,7 +151,7 @@ def __init__( self._register.qubits, device, sampling_rate, - config if config else SimConfig(), + config.to_noise_model() if config else NoiseModel(), ) # Initializing evaluation times self._eval_times_array: np.ndarray @@ -171,7 +172,7 @@ def sampling_times(self) -> np.ndarray: return self._hamiltonian.sampling_times @property - def _sampling_rate(self) -> np.ndarray: + def _sampling_rate(self) -> float: """The sampling rate.""" return self._hamiltonian._sampling_rate @@ -186,14 +187,14 @@ def basis_name(self) -> str: return self._hamiltonian.basis_name @property - def basis(self) -> str: + def basis(self) -> dict[str, Any]: """The basis in which result is expressed.""" return self._hamiltonian.basis @property def config(self) -> SimConfig: """The current configuration, as a SimConfig instance.""" - return self._hamiltonian._config + return SimConfig.from_noise_model(self._hamiltonian.config) def set_config(self, cfg: SimConfig) -> None: """Sets current config to cfg and updates simulation parameters. @@ -213,7 +214,7 @@ def set_config(self, cfg: SimConfig) -> None: " support simulation of noise types:" f"{', '.join(not_supported)}." ) - self._hamiltonian.set_config(cfg) + self._hamiltonian.set_config(cfg.to_noise_model()) def add_config(self, config: SimConfig) -> None: """Updates the current configuration with parameters of another one. @@ -240,7 +241,8 @@ def add_config(self, config: SimConfig) -> None: " support simulation of noise types: " f"{', '.join(not_supported)}." ) - self._hamiltonian.add_config(config) + print("transformed", config.to_noise_model()) + self._hamiltonian.add_config(config.to_noise_model()) def show_config(self, solver_options: bool = False) -> None: """Shows current configuration.""" @@ -248,11 +250,17 @@ def show_config(self, solver_options: bool = False) -> None: def reset_config(self) -> None: """Resets configuration to default.""" - self.set_config(SimConfig()) + self._hamiltonian.set_config(NoiseModel()) @property def initial_state(self) -> qutip.Qobj: - """The initial state of the simulation. + """The initial state of the simulation.""" + return self._initial_state + + def set_initial_state( + self, state: Union[str, np.ndarray, qutip.Qobj] + ) -> None: + """Sets the initial state of the simulation. Args: state: The initial state. @@ -262,12 +270,6 @@ def initial_state(self) -> qutip.Qobj: - An ArrayLike with a shape compatible with the system - A Qobj object """ - return self._initial_state - - def set_initial_state( - self, state: Union[str, np.ndarray, qutip.Qobj] - ) -> None: - """Sets the initial state of the simulation.""" self._initial_state: qutip.Qobj if isinstance(state, str) and state == "all-ground": self._initial_state = qutip.tensor( @@ -295,7 +297,13 @@ def set_initial_state( @property def evaluation_times(self) -> np.ndarray: - """The times at which the results of this simulation are returned. + """The times at which the results of this simulation are returned.""" + return np.array(self._eval_times_array) + + def set_evaluation_times( + self, value: Union[str, ArrayLike, float] + ) -> None: + """Sets times at which the results of this simulation are returned. Args: value: Choose between: @@ -311,12 +319,6 @@ def evaluation_times(self) -> np.ndarray: - A float to act as a sampling rate for the resulting state. """ - return np.array(self._eval_times_array) - - def set_evaluation_times( - self, value: Union[str, ArrayLike, float] - ) -> None: - """Sets times at which the results of this simulation are returned.""" if isinstance(value, str): if value == "Full": eval_times = np.copy(self._hamiltonian.sampling_times) diff --git a/tests/test_sequence_sampler.py b/tests/test_sequence_sampler.py index 2d8c4afb3..0c8376784 100644 --- a/tests/test_sequence_sampler.py +++ b/tests/test_sequence_sampler.py @@ -37,7 +37,9 @@ def assert_same_samples_as_sim(seq: pulser.Sequence) -> None: """Check against the legacy sample extraction in the simulation module.""" got = sample(seq).to_nested_dict() - want = pulser_simulation.QutipEmulator.from_sequence(seq).samples.copy() + want = pulser_simulation.QutipEmulator.from_sequence( + seq + )._hamiltonian.samples.copy() def truncate_samples(samples_dict): for key, value in samples_dict.items(): diff --git a/tests/test_simulation.py b/tests/test_simulation.py index c41b19a1b..d8dda18c8 100644 --- a/tests/test_simulation.py +++ b/tests/test_simulation.py @@ -838,18 +838,20 @@ def test_add_config(matrices): ) with pytest.raises(ValueError, match="is not a valid"): sim.add_config("bad_cfg") - sim.add_config( - SimConfig( - noise=( - "SPAM", - "doppler", - "eff_noise", - ), - eff_noise_opers=[matrices["I"], matrices["X"]], - eff_noise_probs=[0.4, 0.6], - temperature=20000, - ) + with pytest.raises(ValueError, match="is not a valid"): + sim._hamiltonian.add_config("bad_cfg") + config = SimConfig( + noise=( + "SPAM", + "doppler", + "eff_noise", + ), + eff_noise_opers=[matrices["I"], matrices["X"]], + eff_noise_probs=[0.4, 0.6], + temperature=20000, ) + print("init", config) + sim.add_config(config) assert ( "doppler" in sim.config.noise and "SPAM" in sim.config.noise @@ -983,7 +985,14 @@ def test_noisy_xy(): NotImplementedError, match="mode 'XY' does not support simulation of" ): sim.set_config(SimConfig(("SPAM", "doppler"))) - + with pytest.raises(NotImplementedError, match="is not a valid"): + sim._hamiltonian.set_config(SimConfig(("SPAM", "doppler"))) + with pytest.raises( + NotImplementedError, match="mode 'XY' does not support simulation of" + ): + sim._hamiltonian.set_config( + SimConfig(("SPAM", "doppler")).to_noise_model() + ) sim.set_config(SimConfig("SPAM", eta=0.4)) assert sim._hamiltonian._bad_atoms == { "atom0": True, @@ -996,6 +1005,11 @@ def test_noisy_xy(): ): sim.add_config(SimConfig("amplitude")) + with pytest.raises( + NotImplementedError, match="simulation of noise types: amplitude" + ): + sim._hamiltonian.add_config(SimConfig("amplitude").to_noise_model()) + def test_mask_nopulses(): """Check interaction between SLM mask and a simulation with no pulses.""" From b589913e1a4e16eacaed6812c6b59b0f8dbacd32 Mon Sep 17 00:00:00 2001 From: a_corni Date: Wed, 8 Nov 2023 11:11:19 +0100 Subject: [PATCH 10/12] Fixing doc, tests --- pulser-simulation/pulser_simulation/hamiltonian.py | 2 +- pulser-simulation/pulser_simulation/simulation.py | 2 +- tests/test_simconfig.py | 3 ++- tests/test_simulation.py | 2 +- 4 files changed, 5 insertions(+), 4 deletions(-) diff --git a/pulser-simulation/pulser_simulation/hamiltonian.py b/pulser-simulation/pulser_simulation/hamiltonian.py index 81e6a3e8c..c3c58fbb5 100644 --- a/pulser-simulation/pulser_simulation/hamiltonian.py +++ b/pulser-simulation/pulser_simulation/hamiltonian.py @@ -408,7 +408,7 @@ def _update_noise(self) -> None: if "doppler" in self.config.noise_types: detune = np.random.normal( 0, - doppler_sigma(self.config.temperature * 1e-6), + doppler_sigma(self.config.temperature / 1e6), size=len(self._qid_index), ) self._doppler_detune = dict(zip(self._qid_index, detune)) diff --git a/pulser-simulation/pulser_simulation/simulation.py b/pulser-simulation/pulser_simulation/simulation.py index ce8ca2a6a..0a6d2fd98 100644 --- a/pulser-simulation/pulser_simulation/simulation.py +++ b/pulser-simulation/pulser_simulation/simulation.py @@ -79,7 +79,7 @@ def __init__( config: Optional[SimConfig] = None, evaluation_times: Union[float, str, ArrayLike] = "Full", ) -> None: - """Instantiates a Simulation object.""" + """Instantiates a QutipEmulator object.""" # Initializing the samples obj if not isinstance(sampled_seq, SequenceSamples): raise TypeError( diff --git a/tests/test_simconfig.py b/tests/test_simconfig.py index 6662ceabe..af21b92b1 100644 --- a/tests/test_simconfig.py +++ b/tests/test_simconfig.py @@ -16,7 +16,7 @@ from qutip import Qobj, qeye, sigmax, sigmaz from pulser.backend.noise_model import NoiseModel -from pulser_simulation import SimConfig +from pulser_simulation.simconfig import SimConfig, doppler_sigma @pytest.fixture @@ -57,6 +57,7 @@ def test_init(): eff_noise_probs=[0.3, 0.7], ) str_config = config.__str__(True) + assert config.doppler_sigma == doppler_sigma(50 * 1e-6) assert ( "Effective noise distribution" in str_config and "Effective noise operators" in str_config diff --git a/tests/test_simulation.py b/tests/test_simulation.py index d8dda18c8..31bb5ac30 100644 --- a/tests/test_simulation.py +++ b/tests/test_simulation.py @@ -985,7 +985,7 @@ def test_noisy_xy(): NotImplementedError, match="mode 'XY' does not support simulation of" ): sim.set_config(SimConfig(("SPAM", "doppler"))) - with pytest.raises(NotImplementedError, match="is not a valid"): + with pytest.raises(ValueError, match="is not a valid"): sim._hamiltonian.set_config(SimConfig(("SPAM", "doppler"))) with pytest.raises( NotImplementedError, match="mode 'XY' does not support simulation of" From f883996a637e86a2c64789c869ab0c4034311d1d Mon Sep 17 00:00:00 2001 From: a_corni Date: Mon, 27 Nov 2023 12:07:01 +0100 Subject: [PATCH 11/12] Move build_operator, add_config to simulation --- .../pulser_simulation/hamiltonian.py | 57 +---------------- .../pulser_simulation/simulation.py | 62 +++++++++++++++++-- tests/test_simulation.py | 22 +++---- 3 files changed, 65 insertions(+), 76 deletions(-) diff --git a/pulser-simulation/pulser_simulation/hamiltonian.py b/pulser-simulation/pulser_simulation/hamiltonian.py index c3c58fbb5..d69d2fa04 100644 --- a/pulser-simulation/pulser_simulation/hamiltonian.py +++ b/pulser-simulation/pulser_simulation/hamiltonian.py @@ -19,8 +19,7 @@ import warnings from collections import defaultdict from collections.abc import Mapping -from dataclasses import asdict -from typing import Any, Union, cast +from typing import Union, cast import numpy as np import qutip @@ -231,60 +230,6 @@ def set_config(self, cfg: NoiseModel) -> None: # Noise, samples and Hamiltonian update routine self._construct_hamiltonian() - def add_config(self, config: NoiseModel) -> None: - """Updates the current configuration with parameters of another one. - - Mostly useful when dealing with multiple noise types in different - configurations and wanting to merge these configurations together. - Adds simulation parameters to noises that weren't available in the - former NoiseModel. Noises specified in both NoiseModels will keep - former noise parameters. - - Args: - config: NoiseModel to retrieve parameters from. - """ - if not isinstance(config, NoiseModel): - raise ValueError(f"Object {config} is not a valid `NoiseModel`") - - not_supported = ( - set(config.noise_types) - SUPPORTED_NOISES[self._interaction] - ) - if not_supported: - raise NotImplementedError( - f"Interaction mode '{self._interaction}' does not support " - f"simulation of noise types: {', '.join(not_supported)}." - ) - - old_noise_set = set(self.config.noise_types) - new_noise_set = old_noise_set.union(config.noise_types) - diff_noise_set = new_noise_set - old_noise_set - print(diff_noise_set) - # Create temporary param_dict to add noise parameters: - param_dict: dict[str, Any] = asdict(self._config) - # Begin populating with added noise parameters: - param_dict["noise_types"] = tuple(new_noise_set) - if "SPAM" in diff_noise_set: - param_dict["state_prep_error"] = config.state_prep_error - param_dict["p_false_pos"] = config.p_false_pos - param_dict["p_false_neg"] = config.p_false_neg - if "doppler" in diff_noise_set: - param_dict["temperature"] = config.temperature - if "amplitude" in diff_noise_set: - param_dict["laser_waist"] = config.laser_waist - if "dephasing" in diff_noise_set: - param_dict["dephasing_prob"] = config.dephasing_prob - if "depolarizing" in diff_noise_set: - param_dict["depolarizing_prob"] = config.depolarizing_prob - if "eff_noise" in diff_noise_set: - param_dict["eff_noise_opers"] = config.eff_noise_opers - param_dict["eff_noise_probs"] = config.eff_noise_probs - # update runs: - param_dict["runs"] = config.runs - param_dict["samples_per_run"] = config.samples_per_run - - # set config with the new parameters: - self.set_config(NoiseModel(**param_dict)) - def _extract_samples(self) -> None: """Populates samples dictionary with every pulse in the sequence.""" local_noises = True diff --git a/pulser-simulation/pulser_simulation/simulation.py b/pulser-simulation/pulser_simulation/simulation.py index 0a6d2fd98..a1a9e988b 100644 --- a/pulser-simulation/pulser_simulation/simulation.py +++ b/pulser-simulation/pulser_simulation/simulation.py @@ -18,7 +18,7 @@ import warnings from collections import Counter from collections.abc import Mapping -from dataclasses import replace +from dataclasses import asdict, replace from typing import Any, Optional, Union, cast import matplotlib.pyplot as plt @@ -151,7 +151,9 @@ def __init__( self._register.qubits, device, sampling_rate, - config.to_noise_model() if config else NoiseModel(), + config.to_noise_model() + if config + else SimConfig().to_noise_model(), ) # Initializing evaluation times self._eval_times_array: np.ndarray @@ -241,8 +243,34 @@ def add_config(self, config: SimConfig) -> None: " support simulation of noise types: " f"{', '.join(not_supported)}." ) - print("transformed", config.to_noise_model()) - self._hamiltonian.add_config(config.to_noise_model()) + noise_model = config.to_noise_model() + old_noise_set = set(self._hamiltonian.config.noise_types) + new_noise_set = old_noise_set.union(noise_model.noise_types) + diff_noise_set = new_noise_set - old_noise_set + # Create temporary param_dict to add noise parameters: + param_dict: dict[str, Any] = asdict(self._hamiltonian.config) + # Begin populating with added noise parameters: + param_dict["noise_types"] = tuple(new_noise_set) + if "SPAM" in diff_noise_set: + param_dict["state_prep_error"] = noise_model.state_prep_error + param_dict["p_false_pos"] = noise_model.p_false_pos + param_dict["p_false_neg"] = noise_model.p_false_neg + if "doppler" in diff_noise_set: + param_dict["temperature"] = noise_model.temperature + if "amplitude" in diff_noise_set: + param_dict["laser_waist"] = noise_model.laser_waist + if "dephasing" in diff_noise_set: + param_dict["dephasing_prob"] = noise_model.dephasing_prob + if "depolarizing" in diff_noise_set: + param_dict["depolarizing_prob"] = noise_model.depolarizing_prob + if "eff_noise" in diff_noise_set: + param_dict["eff_noise_opers"] = noise_model.eff_noise_opers + param_dict["eff_noise_probs"] = noise_model.eff_noise_probs + # update runs: + param_dict["runs"] = noise_model.runs + param_dict["samples_per_run"] = noise_model.samples_per_run + # set config with the new parameters: + self._hamiltonian.set_config(NoiseModel(**param_dict)) def show_config(self, solver_options: bool = False) -> None: """Shows current configuration.""" @@ -250,7 +278,7 @@ def show_config(self, solver_options: bool = False) -> None: def reset_config(self) -> None: """Resets configuration to default.""" - self._hamiltonian.set_config(NoiseModel()) + self._hamiltonian.set_config(SimConfig().to_noise_model()) @property def initial_state(self) -> qutip.Qobj: @@ -367,6 +395,30 @@ def set_evaluation_times( ) self._eval_times_instruction = value + def build_operator(self, operations: Union[list, tuple]) -> qutip.Qobj: + """Creates an operator with non-trivial actions on some qubits. + + Takes as argument a list of tuples ``[(operator_1, qubits_1), + (operator_2, qubits_2)...]``. Returns the operator given by the tensor + product of {``operator_i`` applied on ``qubits_i``} and Id on the rest. + ``(operator, 'global')`` returns the sum for all ``j`` of operator + applied at ``qubit_j`` and identity elsewhere. + + Example for 4 qubits: ``[(Z, [1, 2]), (Y, [3])]`` returns `ZZYI` + and ``[(X, 'global')]`` returns `XIII + IXII + IIXI + IIIX` + + Args: + operations: List of tuples `(operator, qubits)`. + `operator` can be a ``qutip.Quobj`` or a string key for + ``self.op_matrix``. `qubits` is the list on which operator + will be applied. The qubits can be passed as their + index or their label in the register. + + Returns: + The final operator. + """ + return self._hamiltonian.build_operator(operations) + def get_hamiltonian(self, time: float) -> qutip.Qobj: r"""Get the Hamiltonian created from the sequence at a fixed time. diff --git a/tests/test_simulation.py b/tests/test_simulation.py index 31bb5ac30..614cd988b 100644 --- a/tests/test_simulation.py +++ b/tests/test_simulation.py @@ -262,15 +262,15 @@ def test_building_basis_and_projection_operators(seq, reg): # Check local operator building method: with pytest.raises(ValueError, match="Duplicate atom"): - sim._hamiltonian.build_operator([("sigma_gg", ["target", "target"])]) + sim.build_operator([("sigma_gg", ["target", "target"])]) with pytest.raises(ValueError, match="not a valid operator"): - sim._hamiltonian.build_operator([("wrong", ["target"])]) + sim.build_operator([("wrong", ["target"])]) with pytest.raises(ValueError, match="Invalid qubit names: {'wrong'}"): - sim._hamiltonian.build_operator([("sigma_gg", ["wrong"])]) + sim.build_operator([("sigma_gg", ["wrong"])]) # Check building operator with one operator - op_standard = sim._hamiltonian.build_operator([("sigma_gg", ["target"])]) - op_one = sim._hamiltonian.build_operator(("sigma_gg", ["target"])) + op_standard = sim.build_operator([("sigma_gg", ["target"])]) + op_one = sim.build_operator(("sigma_gg", ["target"])) assert np.linalg.norm(op_standard - op_one) < 1e-10 # Global ground-rydberg @@ -838,8 +838,6 @@ def test_add_config(matrices): ) with pytest.raises(ValueError, match="is not a valid"): sim.add_config("bad_cfg") - with pytest.raises(ValueError, match="is not a valid"): - sim._hamiltonian.add_config("bad_cfg") config = SimConfig( noise=( "SPAM", @@ -850,7 +848,6 @@ def test_add_config(matrices): eff_noise_probs=[0.4, 0.6], temperature=20000, ) - print("init", config) sim.add_config(config) assert ( "doppler" in sim.config.noise @@ -1005,11 +1002,6 @@ def test_noisy_xy(): ): sim.add_config(SimConfig("amplitude")) - with pytest.raises( - NotImplementedError, match="simulation of noise types: amplitude" - ): - sim._hamiltonian.add_config(SimConfig("amplitude").to_noise_model()) - def test_mask_nopulses(): """Check interaction between SLM mask and a simulation with no pulses.""" @@ -1208,7 +1200,7 @@ def test_effective_size_intersection(): "atom2": True, "atom3": False, } - assert sim.get_hamiltonian(0) != 0 * sim._hamiltonian.build_operator( + assert sim.get_hamiltonian(0) != 0 * sim.build_operator( [("I", "global")] ) @@ -1242,7 +1234,7 @@ def test_effective_size_disjoint(channel_type): if channel_type == "mw_global": assert sim.get_hamiltonian( 0 - ) == 0.5 * amp * sim._hamiltonian.build_operator( + ) == 0.5 * amp * sim.build_operator( [(qutip.sigmax(), ["atom3"])] ) else: From 5a489528267ec04053aecb57f390a8ff389d37af Mon Sep 17 00:00:00 2001 From: a_corni Date: Mon, 4 Dec 2023 16:47:18 +0100 Subject: [PATCH 12/12] Fix typing --- tests/test_simulation.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/tests/test_simulation.py b/tests/test_simulation.py index 318968281..3b5dd0097 100644 --- a/tests/test_simulation.py +++ b/tests/test_simulation.py @@ -1231,9 +1231,7 @@ def test_effective_size_disjoint(channel_type): "atom3": False, } if channel_type == "mw_global": - assert sim.get_hamiltonian( - 0 - ) == 0.5 * amp * sim.build_operator( + assert sim.get_hamiltonian(0) == 0.5 * amp * sim.build_operator( [(qutip.sigmax(), ["atom3"])] ) else: