refactoring emulator results.py file

doc/source/main-documentation/emulator.rst

@@ -0,0 +1,25 @@
+.. admonition:: Work in progress
+
+This documentation is currently in draft form and may be incomplete.
+
+Emulator
+========
+
+Qibolab provides an internal simulation instrument that enables the emulation of a variety of quantum chip configurations, including systems with single or multiple qubits, fixed- or tunable-frequency architectures, and optional couplers. This emulator allows users to execute virtual Qibolab experiments in a manner fully consistent with their execution on real quantum processing units (QPUs).
+
+For a more detailed description of how Qibolab handles different QPUs (here referred to as :class:`.Platforms`), the reader is referred to :ref:`Platform guidelines <main_doc_platform>`.
+
+Usage
+-----
+
+The emulator leverages a third-party numerical engine that solves the Master Equation for a system governed by the sum of a time-independent Hamiltonian and a time-dependent (pulse) Hamiltonian. The solver produces the system density matrix at each time step. From this solution, the emulator extracts only those density matrices corresponding to acquisition pulses defined within the pulse sequence.
+
+Since the simulator operates at the level of density matrices, it does not reproduce in-phase and quadrature (I-Q) measurement signals as in real experimental setups. Instead, it computes the measurement probabilities :math:`p_m = \bra{m} \rho \ket{m}` for each computational basis state :math:`\ket{m}`. Consequently, while the emulator supports all experiment types, in signal-based experiments (i.e., when :paramref:`AcquisitionType.INTEGRATION` is selected), the signal magnitude corresponds directly to these probabilities, whereas the signal phase carries no physical meaning.
+
+The emulator supports both :paramref:`AveragingMode.SINGLESHOT` and :paramref:`AveragingMode.CYCLIC`. In the former case, the simulator returns discrete measurement outcomes corresponding to a finite number of shots, whereas in the latter it returns expectation values derived from the density matrix, typically the probability of the :math:`\ket{1}` state. Although SINGLESHOT mode includes statistical sampling noise, it is often computationally advantageous to bypass sampling and directly return the diagonal elements of the density matrix. In CYCLIC mode, Gaussian noise is additionally introduced, with a fixed standard deviation.
+
+To accurately resolve qubit dynamics and capture all contributions from the time-dependent Hamiltonian (including control pulses), a Nyquist frequency is defined. By default, this is set to :math:`f_N = 20 \text{GHz}`, which allows accurate resolution of oscillations up to approximately :math:`10–15 \text{GHz}`. The choice of Nyquist frequency is critical, as it determines the temporal resolution of the simulation and informs adaptive tuning of the ODE solver parameters. Further details on this tuning procedure are provided in the implementation of the numerical engines.
+
+At present, state collapse is not implemented in the Qibolab emulator. As a result, this tool is not suitable for simulating mid-circuit measurements, and should be restricted to circuits in which all measurements occur simultaneously at the end of the computation. In physical systems, mid-circuit measurement induces wavefunction collapse; if the measured qubit is entangled with others, this process introduces correlations that condition the state of the remaining system. The current emulator does not account for such measurement-induced correlations, leading to intrinsically inaccurate results in these scenarios.
+
+An additional limitation arises from the handling of measurement ordering. In certain experimental protocols, the temporal order of measurements may vary across parameter sweeps, while the :paramref:`PulseSequence` object remains fixed and does not reflect such reordering. This discrepancy may introduce inconsistencies during the execution of :func:`qibolab._core.instruments.emulator.results.results`, which assumes alignment between the time-ordering structure and the acquisition pulses defined in the pulse sequence. If this alignment is violated, acquisition events may be incorrectly matched to simulation time steps, ultimately resulting in invalid outputs.

src/qibolab/_core/instruments/emulator/emulator.py

@@ -11,7 +11,7 @@
 
 from qibolab._core.components import Config
 from qibolab._core.components.configs import AcquisitionConfig
-from qibolab._core.execution_parameters import ExecutionParameters
+from qibolab._core.execution_parameters import AveragingMode, ExecutionParameters
 from qibolab._core.identifier import Result
 from qibolab._core.instruments.abstract import Controller
 from qibolab._core.pulses import (
@@ -65,21 +65,44 @@ def play(
         options: ExecutionParameters,
         sweepers: list[ParallelSweepers],
     ) -> dict[int, Result]:
+
+        if (
+            options.averaging_mode is AveragingMode.SINGLESHOT
+            and options.nshots is None
+        ):
+            raise ValueError("nshots must be specified for SINGLESHOT mode")
+
         # convert align to delays
         sequences_ = (seq.align_to_delays() for seq in sequences)
-        # just merge the results of multiple executions in a single dictionary
-        return reduce(
-            or_,
-            (
-                results(
-                    # states in computational basis
-                    self._sweep(sequence, configs, sweepers),
-                    sequence,
-                    cast(HamiltonianConfig, configs["hamiltonian"]),
-                    options,
-                )
-                for sequence in sequences_
-            ),
+
+        results_to_process = (
+            self._play_sequence(configs, sequence, options, sweepers)
+            for sequence in sequences_
+        )
+
+        return reduce(or_, results_to_process)
+
+    def _play_sequence(
+        self,
+        configs: dict[str, Config],
+        sequence: PulseSequence,
+        options: ExecutionParameters,
+        sweepers: list[ParallelSweepers],
+    ):
+        """
+        Generate results from an emulated quantum sequence execution.
+        Executes a sweep of the quantum sequence and processes the results
+        into a structured results object containing quantum states and measurement data.
+        """
+
+        sweep_results = self._sweep(sequence, configs, sweepers)
+        hamiltonian = cast(HamiltonianConfig, configs["hamiltonian"])
+        return results(
+            # states in computational basis
+            states=sweep_results,
+            sequence=sequence,
+            hamiltonian=hamiltonian,
+            options=options,
         )
 
     def _sweep(
@@ -100,7 +123,7 @@ def _sweep(
         updates = defaultdict(dict) | ({} if updates is None else updates)
 
         if len(sweepers) == 0:
-            return self._play_sequence(sequence, configs, updates)
+            return self._evolve(sequence, configs, updates)
 
         parsweep = sweepers[0]
         # collect slices of results, corresponding to the current iteration
@@ -122,13 +145,14 @@ def _sweep(
         # stack all slices in a single array, along the current outermost dimension
         return np.stack(results)
 
-    def _play_sequence(
+    def _evolve(
         self, sequence: PulseSequence, configs: dict[str, Config], updates: dict
     ) -> NDArray:
-        """Play single sequence on emulator.
+        """Evolve a pulse sequence on the quantum emulator.
 
-        The array returned by this function has a single dimension, over
-        the various measurements included in the sequence.
+        This method updates the sequence parameters, generates the time grid, constructs
+        the time-dependent Hamiltonian, evolves the initial state with optional collapse
+        operators, and returns the resulting measurement data.
         """
         sequence_ = update_sequence(sequence, updates)
         tlist_ = tlist(sequence_, self.sampling_rate, per_sample=2)

src/qibolab/_core/instruments/emulator/results.py

@@ -1,3 +1,16 @@
+"""
+In this module we often recall, for the sake of a better interpretation of the whole pipeline, the transformed dimensions of the simulation's results tensor.
+Generally this tensor will have the shape (*S, M *H_dim) (or its permutations),
+where:
+- *S is the number of iteration for each sweep in the experiment
+- M is the number of measurements applied in the pulse sequence
+- *H_dim is the complete system dimension
+
+In the results processing also, when simulating SINGLESHOTS, we'll add two dimensions:
+- Nshots, which is simply the number of shots we average on
+- M_unique, which is the number of unique measurement times
+"""
+
 from collections.abc import Iterable
 
 import numpy as np
@@ -8,7 +21,7 @@
     AveragingMode,
     ExecutionParameters,
 )
-from ...identifier import ChannelId, QubitId, Result
+from ...identifier import ChannelId, Result
 from ...pulses import Acquisition, Align, PulseId, Readout
 from ...sequence import PulseSequence
 from .engine import Operator
@@ -51,19 +64,31 @@ def shots(probabilities: NDArray, nshots: int) -> NDArray:
     return np.moveaxis(shots, -1, 0)
 
 
-def calculate_probabilities_from_density_matrix(
-    states: NDArray, subsystems: Iterable[QubitId], nsubsystems: int, d: int
-) -> NDArray:
-    """Compute probabilities from density matrix."""
-    states_ = np.reshape(states, states.shape[:-2] + 2 * nsubsystems * (d,))
-    marginal = np.einsum(
-        states_,
+def _extract_probabilities(states: NDArray) -> NDArray:
+    """
+    Calculate probabilities from a density matrix using diagonal elements.
+
+    This function extracts the diagonal elements of each density matrix,
+    which represent the probabilities of measurement outcomes.
+
+    Probabilities are normalized for fluctuations, taking the absolute value of diagonal elements.
+
+    Examples
+    --------
+    >>> dm = np.array([[0.9, 0.0], [0.0, 0.1]])
+    >>> probs = _extract_probabilities(dm)
+    >>> probs
+    array([0.9, 0.1])
+    """
+
+    diag = np.einsum(
+        states,
         # TODO: the `np.array()` wrapping call is only needed because of NumPy's type
         # annotation - in practice, it also works without
-        np.array([...] + list(range(nsubsystems)) * 2),
-        np.array([...] + list(subsystems)),
+        np.array([...] + [0, 0]),
+        np.array([...] + [0]),
     )
-    return np.abs(marginal).reshape((*states.shape[:-2], -1))
+    return np.clip(diag.real, 0, 1)
 
 
 def acquisitions(sequence: PulseSequence) -> dict[PulseId, float]:
@@ -92,17 +117,118 @@ def index(ch: ChannelId, hconfig: HamiltonianConfig) -> int:
 def select_acquisitions(
     states: list[Operator], acquisitions: Iterable[float], times: NDArray
 ) -> NDArray:
-    """Select density matrices from states.
+    """
+    Select and organize quantum state acquisitions based on acquisition times.
 
-    First, retrieve acquisitions, and locate them in the tlist, to
-    isolate the expectations related to measurements.
+    This function filters quantum states corresponding to specified acquisition times
+    and maps them to unique acquisition values. It uses binary search to find the
+    nearest state index for each acquisition time.
 
-    The return type should be rank-3 array, where the last two are the density
-    matrices dimensions, while the first one should correspond to the acquisitions.
+    It returns a NumPy array containing the full density matrices
+    of the selected quantum states, indexed by the original acquisition order.
     """
-    acq = np.array(list(acquisitions))
+    acq, index_pos = np.unique(list(acquisitions), return_inverse=True)
     samples = np.minimum(np.searchsorted(times, acq), times.size - 1)
-    return np.stack([states[n].full() for n in samples])
+    return np.stack([states[n].full() for n in samples])[index_pos]
+
+
+def _cyclic_results(
+    state_probs: NDArray,
+    sequence: PulseSequence,
+    hamiltonian: HamiltonianConfig,
+    options: ExecutionParameters,
+) -> dict[int, Result]:
+    """Process measurement results from a cyclic quantum simulation, where the output for each measurement is
+    excited state population.
+    Computes readout results by projecting quantum state probabilities onto
+    measurement subspaces and applying configured post-processing.
+    """
+
+    # Through the entire function state_probs has dimensions:
+    # (*S, M *H_dim)
+    states_computational_idx = np.stack(
+        np.unravel_index(np.arange(state_probs.shape[-1]), hamiltonian.dims)
+    )
+
+    acq_id = acquisitions(sequence).keys()
+    # from every acquisition pulse id we get the corresponding channel, and from the channel we get the
+    # corresponding qubit index, which is then used to correctly permute the rows of states_computational_idx.
+    qubit_indices = [
+        index(sequence.pulse_channels(ro_id)[0], hamiltonian) for ro_id in acq_id
+    ]
+    permuted_states_computational_idx = states_computational_idx[qubit_indices]
+
+    # applying a mask to select for each measurement the states that are outside the computational subspace, which are classified as 1
+    mask = permuted_states_computational_idx >= 1
+
+    # res is a (M, *S, ...) array
+    res = np.moveaxis(np.sum(np.where(mask, state_probs, 0), axis=-1), -1, 0)
+
+    if options.acquisition_type is AcquisitionType.INTEGRATION:
+        res = np.random.normal(res, scale=0.001)
+        zeros = np.zeros(res.shape) if np.ndim(res) != 0 else 0.0
+        res = np.stack((res, zeros), axis=-1)
+
+    return dict(zip(acq_id, res))
+
+
+def _singleshot_results(
+    state_probs: NDArray,
+    sequence: PulseSequence,
+    hamiltonian: HamiltonianConfig,
+    options: ExecutionParameters,
+) -> dict[int, Result]:
+    """Extract measurement results from simulated quantum state probabilities.
+    Performs single-shot measurement extraction from state probabilities by sampling
+    according to the specified number of shots and mapping readout operations to their
+    corresponding measurement results.
+    """
+
+    # select only unique times of measurements
+    _, direct_map, inverse_map = np.unique(
+        list(acquisitions(sequence).values()),
+        return_index=True,
+        return_inverse=True,
+    )
+
+    # here we move the -2 index of the probability array, hence it now becomes:
+    # (M, *S, *H_dim)
+    state_probs = np.moveaxis(state_probs, -2, 0)
+
+    # apply the direct mapping found in np.unique to the probability vector
+    # in order to sample at unique times and hence mantain correlation
+    # since we use the same shots for synchronous measurements.
+    unique_state_probs = state_probs[direct_map]
+
+    # shots function returns a vector of shape:
+    # (Nshots, M_unique, *S, *H_dim)
+    sampled = shots(unique_state_probs, options.nshots)
+
+    # move measurements dimension to the front, getting ready for extraction
+    # the shape now is: (M_unique, Nshots, *S, *H_dim)
+    sampled = np.moveaxis(sampled, 1, 0)
+
+    acq_id = acquisitions(sequence).keys()
+    # from every acquisition pulse id we get the corresponding channel, and from the channel we get the
+    # corresponding qubit index, which is then used to correctly permute the rows of states_computational_idx.
+    qubit_indices = [
+        index(sequence.pulse_channels(ro_id)[0], hamiltonian) for ro_id in acq_id
+    ]
+
+    # we use inverse_map to expand back the sampled results
+    # res is a (M, M, Nshots, *S, ...) array
+    res = np.stack(np.unravel_index(sampled[inverse_map], hamiltonian.dims))[
+        qubit_indices
+    ]
+    # using np.einsum, so res is a (M, Nshots, *S, ...) array
+    res = np.einsum(res, np.array([0, 0] + [...]), np.array([0] + [...]))
+    res = np.clip(res, 0, 1)
+
+    if options.acquisition_type is AcquisitionType.INTEGRATION:
+        zeros = np.zeros(res.shape) if np.ndim(res) != 0 else 0.0
+        res = np.stack((res, zeros), axis=-1)
+
+    return dict(zip(acq_id, res))
 
 
 def results(
@@ -117,34 +243,20 @@ def results(
     result for the execution of this single sequence, thus suitable
     to be returned as is.
     """
-    probabilities = calculate_probabilities_from_density_matrix(
-        states,
-        range(hamiltonian.nqubits),
-        hamiltonian.nqubits,
-        hamiltonian.transmon_levels,
-    )
-    assert options.nshots is not None
-    sampled = shots(np.moveaxis(probabilities, -2, 0), options.nshots)
-    # move measurements dimension to the front, getting ready for extraction
-    measurements = np.moveaxis(sampled, 1, 0)
-
-    results = {}
-    # introduce cached measurements to avoid losing correlations
-    cache_measurements = {}
-    for (ro_id, sample), meas in zip(acquisitions(sequence).items(), measurements):
-        i = index(sequence.pulse_channels(ro_id)[0], hamiltonian)
-        cache_measurements.setdefault(sample, meas)
-        res = np.stack(np.unravel_index(cache_measurements[sample], hamiltonian.dims))[
-            i
-        ]
 
-        if options.acquisition_type is AcquisitionType.INTEGRATION:
-            res = np.stack((res, np.zeros_like(res)), axis=-1)
-            res = np.random.normal(res, scale=0.001)
+    # probability dimensions are:
+    # (*S, M, *H_dim)
+    probabilities = _extract_probabilities(states)
 
-        if options.averaging_mode == AveragingMode.CYCLIC:
-            res = np.mean(res, axis=0)
-
-        results[ro_id] = res
+    results = (
+        _singleshot_results
+        if options.averaging_mode is AveragingMode.SINGLESHOT
+        else _cyclic_results
+    )
 
-    return results
+    return results(
+        state_probs=probabilities,
+        sequence=sequence,
+        hamiltonian=hamiltonian,
+        options=options,
+    )