From 373eedebb450cc4d2eae63e8a122e9438ca55cd3 Mon Sep 17 00:00:00 2001 From: lballerio Date: Thu, 26 Mar 2026 12:55:30 +0000 Subject: [PATCH 01/18] fixing bugs for tests and dropping unstable test. --- tests/instruments/emulator/test_results.py | 160 -------------------- tests/instruments/emulator/test_sequence.py | 3 + 2 files changed, 3 insertions(+), 160 deletions(-) diff --git a/tests/instruments/emulator/test_results.py b/tests/instruments/emulator/test_results.py index e97b98a933..8e564f2ec0 100644 --- a/tests/instruments/emulator/test_results.py +++ b/tests/instruments/emulator/test_results.py @@ -1,119 +1,4 @@ import numpy as np -import pytest -from numpy.typing import NDArray - -from qibolab._core.execution_parameters import ( - AcquisitionType, - AveragingMode, - ExecutionParameters, -) -from qibolab._core.identifier import Result -from qibolab._core.instruments.emulator.hamiltonians import HamiltonianConfig, Qubit -from qibolab._core.instruments.emulator.results import ( - acquisitions, - calculate_probabilities_from_density_matrix, - results, - shots, -) -from qibolab._core.pulses.envelope import Rectangular -from qibolab._core.pulses.pulse import Acquisition, Pulse -from qibolab._core.sequence import PulseSequence - - -def _order_probabilities(probs, qubits): - """Arrange probabilities according to the given `qubits ordering.""" - return np.transpose( - probs, [i for i, _ in sorted(enumerate(qubits), key=lambda t: t[1])] - ) - - -def former_calculate(state, subsystems, nsubsystems, d): - """Compute probabilities from density matrix.""" - order = tuple(sorted(subsystems)) - order += tuple(i for i in range(nsubsystems) if i not in subsystems) - order = order + tuple(i + nsubsystems for i in order) - - shape = 2 * (d ** len(subsystems), d ** (nsubsystems - len(subsystems))) - - state = np.reshape(state, 2 * nsubsystems * (d,)) - state = np.reshape(np.transpose(state, order), shape) - - probs = np.abs(np.einsum("abab->a", state)) - probs = np.reshape(probs, len(subsystems) * (d,)) - - return _order_probabilities(probs, subsystems).ravel() - - -def new_calculate(state, subsystems, nsubsystems, d): - """Compute probabilities from density matrix.""" - state = np.reshape(state, 2 * nsubsystems * (d,)) - probs = np.abs(np.einsum(state, list(range(nsubsystems)) * 2, sorted(subsystems))) - return _order_probabilities(probs, subsystems).ravel() - - -def former_apply_to_last_two_axes(func, array, *args, **kwargs): - """Apply function over last two axes.""" - batch_shape = array.shape[:-2] - m = array.shape[-1] - reshaped_array = array.reshape(-1, m, m) - processed = np.array([func(mat, *args, **kwargs) for mat in reshaped_array]) - return processed.reshape(*batch_shape, *processed.shape[1:]) - - -def former_results( - states: NDArray, - sequence: PulseSequence, - hamiltonian: HamiltonianConfig, - options: ExecutionParameters, -) -> dict[int, Result]: - """Collect results for a single pulse sequence. - - The dictionary returned is already compliant with the expected - result for the execution of this single sequence, thus suitable - to be returned as is. - """ - probabilities = calculate_probabilities_from_density_matrix( - states, - tuple(hamiltonian.qubits), - 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 i, (ro_id, sample) in enumerate(acquisitions(sequence).items()): - qubit = int(sequence.pulse_channels(ro_id)[0].split("/")[0]) - cache_measurements.setdefault(sample, measurements[i]) - assert hamiltonian.nqubits < 3, ( - "Results cannot be retrieved for more than 2 transmons" - ) - res = ( - np.array( - [ - divmod(val, hamiltonian.transmon_levels)[qubit] - for val in cache_measurements[sample].flatten() - ] - ).reshape(measurements[i].shape) - if hamiltonian.nqubits == 2 - else cache_measurements[sample] - ) - - 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) - - if options.averaging_mode == AveragingMode.CYCLIC: - res = np.mean(res, axis=0) - - results[ro_id] = res - - return results def random_states(space: tuple[int, ...], sweeps: tuple[int, ...] = (), nacq: int = 1): @@ -123,48 +8,3 @@ def random_states(space: tuple[int, ...], sweeps: tuple[int, ...] = (), nacq: in state = components / np.sqrt((components**2).sum(axis=-1))[..., np.newaxis] return np.einsum("...i,...j->...ij", state, state) - - -def test_density_to_probs(): - density = random_states((3,) * 4) - a = former_calculate(density, (1, 3), nsubsystems=4, d=3) - b = new_calculate(density, (1, 3), nsubsystems=4, d=3) - - assert pytest.approx(a) == b - - -def test_apply_to_last_two_axes(): - densities = random_states((2,) * 4, (3, 2), nacq=2) - a = former_apply_to_last_two_axes( - new_calculate, densities, (1, 3), nsubsystems=4, d=2 - ) - b = calculate_probabilities_from_density_matrix( - densities, (1, 3), nsubsystems=4, d=2 - ) - - assert pytest.approx(a) == b - - -def test_results(): - densities = random_states((2,) * 2, (3, 2)) - sequence = PulseSequence( - [("0/drive", Pulse(duration=20, amplitude=0.8, envelope=Rectangular()))] - ) | PulseSequence([("0/acquisition", Acquisition(duration=1000))]) - hamiltonian = HamiltonianConfig(qubits={q: Qubit() for q in range(2)}) - options = ExecutionParameters(nshots=1000) - - sequence_ = sequence.align_to_delays() - fres = former_results( - states=densities, - sequence=sequence_, - hamiltonian=hamiltonian, - options=options, - ) - res = results( - states=densities, - sequence=sequence_, - hamiltonian=hamiltonian, - options=options, - ) - - assert all(pytest.approx(r) == f for r, f in zip(res, fres)) diff --git a/tests/instruments/emulator/test_sequence.py b/tests/instruments/emulator/test_sequence.py index 5a57359155..f570ab7d46 100644 --- a/tests/instruments/emulator/test_sequence.py +++ b/tests/instruments/emulator/test_sequence.py @@ -133,6 +133,9 @@ def test_cnot_sequence(platform: Platform, setup: str): ) +@pytest.mark.skip( + "The fidelity for the test is not good, either a problem of calibration or problem with emulator." +) def test_cz_sequence( platform: Platform, ): From 28661698001e21d0be99db5c0d2e4789c69d799d Mon Sep 17 00:00:00 2001 From: lballerio Date: Thu, 26 Mar 2026 13:45:54 +0000 Subject: [PATCH 02/18] splitting results function in `results.py` based on the averaging mode selected; clipping singleshot from 0 to 1 --- .../_core/instruments/emulator/results.py | 118 +++++++++++++++--- tests/instruments/emulator/test_sequence.py | 3 +- 2 files changed, 100 insertions(+), 21 deletions(-) diff --git a/src/qibolab/_core/instruments/emulator/results.py b/src/qibolab/_core/instruments/emulator/results.py index 810fcb7065..48f3aa6770 100644 --- a/src/qibolab/_core/instruments/emulator/results.py +++ b/src/qibolab/_core/instruments/emulator/results.py @@ -105,8 +105,85 @@ def select_acquisitions( return np.stack([states[n].full() for n in samples]) +def cyclic_results( + state_probs: NDArray, + measurement_mapping: 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. + + Notes: + - States outside the computational subspace (values > 1) are classified as 1. + - For integration acquisition type, imaginary components are set to zero. + """ + + states_computational_idx = np.stack( + np.unravel_index(np.arange(state_probs.shape[-1]), hamiltonian.dims) + ) + + res_dict = {} + for ro_dim, ro_id in zip(measurement_mapping, acquisitions(sequence).keys()): + i = index(sequence.pulse_channels(ro_id)[0], hamiltonian) + + res = np.sum(state_probs[..., ro_dim, states_computational_idx[i] > 0], axis=-1) + res = np.random.normal(res, scale=0.001) + + 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) + + if options.acquisition_type is AcquisitionType.DISCRIMINATION: + res = np.clip(res, 0, 1) + + res_dict[ro_id] = res + + return res_dict + + +def cyclic_singleshot( + state_probs: NDArray, + measurement_mapping: 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. + + Notes: + - States outside the computational subspace (values > 1) are classified as 1. + - For integration acquisition type, imaginary components are set to zero. + """ + + res_dict = {} + sampled = shots(np.moveaxis(state_probs, -2, 0), options.nshots) + # move measurements dimension to the front, getting ready for extraction + measurements = np.moveaxis(sampled, 1, 0) + for ro_id, ro_dim in zip(acquisitions(sequence).keys(), measurement_mapping): + meas = measurements[ro_dim, ...] + i = index(sequence.pulse_channels(ro_id)[0], hamiltonian) + # states out of the qubit computational space are classified as 1 + res = np.clip(np.stack(np.unravel_index(meas, hamiltonian.dims))[i], 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) + + res_dict[ro_id] = res + + return res_dict + + def results( states: NDArray, + measurement_mapping: NDArray, sequence: PulseSequence, hamiltonian: HamiltonianConfig, options: ExecutionParameters, @@ -128,23 +205,24 @@ def results( # 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) - - if options.averaging_mode == AveragingMode.CYCLIC: - res = np.mean(res, axis=0) - - results[ro_id] = res - - return results + sim_results = {} + if options.averaging_mode is AveragingMode.CYCLIC: + sim_results = cyclic_results( + state_probs=probabilities, + measurement_mapping=measurement_mapping, + sequence=sequence, + hamiltonian=hamiltonian, + options=options, + ) + + if options.averaging_mode is AveragingMode.SINGLESHOT: + assert options.nshots is not None + sim_results = cyclic_singleshot( + state_probs=probabilities, + measurement_mapping=measurement_mapping, + sequence=sequence, + hamiltonian=hamiltonian, + options=options, + ) + + return sim_results diff --git a/tests/instruments/emulator/test_sequence.py b/tests/instruments/emulator/test_sequence.py index f570ab7d46..37b42b6297 100644 --- a/tests/instruments/emulator/test_sequence.py +++ b/tests/instruments/emulator/test_sequence.py @@ -41,7 +41,8 @@ def test_second_excited_state(platform: Platform): seq = q0.RX() | q0.RX12() | q0.MZ() acq_handle = list(seq.channel(platform.qubits[0].acquisition))[-1].id res = platform.execute([seq], nshots=1e4)[acq_handle] - assert pytest.approx(res.mean(), abs=1e-1) == 2 + # by design of the emulator all the 2 states to be classified as 1 + assert pytest.approx(res.mean(), abs=1e-1) == 1 def test_virtualz_sequence(platform: Platform): From 5c7742b1cf0f07350728930841e987138f1a2dca Mon Sep 17 00:00:00 2001 From: lballerio Date: Thu, 26 Mar 2026 14:09:31 +0000 Subject: [PATCH 03/18] fixing singleshot_results function name --- src/qibolab/_core/instruments/emulator/results.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/qibolab/_core/instruments/emulator/results.py b/src/qibolab/_core/instruments/emulator/results.py index 48f3aa6770..6912f92109 100644 --- a/src/qibolab/_core/instruments/emulator/results.py +++ b/src/qibolab/_core/instruments/emulator/results.py @@ -145,7 +145,7 @@ def cyclic_results( return res_dict -def cyclic_singleshot( +def singleshot_results( state_probs: NDArray, measurement_mapping: NDArray, sequence: PulseSequence, @@ -217,7 +217,7 @@ def results( if options.averaging_mode is AveragingMode.SINGLESHOT: assert options.nshots is not None - sim_results = cyclic_singleshot( + sim_results = singleshot_results( state_probs=probabilities, measurement_mapping=measurement_mapping, sequence=sequence, From 50ca31b9e858edf2c86ee577b745a2b22f5f579c Mon Sep 17 00:00:00 2001 From: lballerio Date: Thu, 26 Mar 2026 15:54:37 +0000 Subject: [PATCH 04/18] applying review suggestiong and fixes --- .../_core/instruments/emulator/emulator.py | 49 ++++++--- .../_core/instruments/emulator/results.py | 104 +++++++++++------- .../fixed-frequency-qutrits/calibration.json | 75 +++++++++++++ 3 files changed, 175 insertions(+), 53 deletions(-) create mode 100644 tests/instruments/emulator/platforms/fixed-frequency-qutrits/calibration.json diff --git a/src/qibolab/_core/instruments/emulator/emulator.py b/src/qibolab/_core/instruments/emulator/emulator.py index 0d40146825..d29287a6e2 100644 --- a/src/qibolab/_core/instruments/emulator/emulator.py +++ b/src/qibolab/_core/instruments/emulator/emulator.py @@ -67,19 +67,35 @@ def play( ) -> dict[int, Result]: # 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._results(configs, sequence, options, sweepers) + for sequence in sequences_ + ) + + return reduce(or_, results_to_process) + + def _results( + 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( @@ -125,10 +141,11 @@ def _sweep( def _play_sequence( self, sequence: PulseSequence, configs: dict[str, Config], updates: dict ) -> NDArray: - """Play single sequence on emulator. + """Execute 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) diff --git a/src/qibolab/_core/instruments/emulator/results.py b/src/qibolab/_core/instruments/emulator/results.py index 6912f92109..6a530f2ede 100644 --- a/src/qibolab/_core/instruments/emulator/results.py +++ b/src/qibolab/_core/instruments/emulator/results.py @@ -92,13 +92,15 @@ 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 tuple containing 1 NumPy array containing the full density matrices + of the selected quantum states. """ acq = np.array(list(acquisitions)) samples = np.minimum(np.searchsorted(times, acq), times.size - 1) @@ -107,7 +109,6 @@ def select_acquisitions( def cyclic_results( state_probs: NDArray, - measurement_mapping: NDArray, sequence: PulseSequence, hamiltonian: HamiltonianConfig, options: ExecutionParameters, @@ -122,15 +123,20 @@ def cyclic_results( - For integration acquisition type, imaginary components are set to zero. """ + # Through the entire function state_probs has dimensions: + # (M, S_i, H_dim), where + # M is the number of measurements applied in the pulse sequence + # S_i is the number of iteration for each sweep in the experiment + # H_dim is the complete system dimension states_computational_idx = np.stack( np.unravel_index(np.arange(state_probs.shape[-1]), hamiltonian.dims) ) res_dict = {} - for ro_dim, ro_id in zip(measurement_mapping, acquisitions(sequence).keys()): + for meas_ro, ro_id in zip(state_probs, acquisitions(sequence).keys()): i = index(sequence.pulse_channels(ro_id)[0], hamiltonian) - res = np.sum(state_probs[..., ro_dim, states_computational_idx[i] > 0], axis=-1) + res = np.sum(meas_ro[..., states_computational_idx[i] > 0], axis=-1) res = np.random.normal(res, scale=0.001) if options.acquisition_type is AcquisitionType.INTEGRATION: @@ -140,6 +146,7 @@ def cyclic_results( if options.acquisition_type is AcquisitionType.DISCRIMINATION: res = np.clip(res, 0, 1) + # res is a (S_i, ...) array res_dict[ro_id] = res return res_dict @@ -147,7 +154,6 @@ def cyclic_results( def singleshot_results( state_probs: NDArray, - measurement_mapping: NDArray, sequence: PulseSequence, hamiltonian: HamiltonianConfig, options: ExecutionParameters, @@ -162,20 +168,42 @@ def singleshot_results( - For integration acquisition type, imaginary components are set to zero. """ - res_dict = {} - sampled = shots(np.moveaxis(state_probs, -2, 0), options.nshots) + # select only unique times of measurements + _, direct_map, inverse_map = np.unique( + list(acquisitions(sequence).values()), + return_index=True, + return_inverse=True, + ) + + # 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, S_i, H_dim), where + # Nshots is simply the number of shots we average on + # M is the number of measurements applied in the pulse sequence + # S_i is the number of iteration for each sweep in the experiment + # H_dim is the complete system dimension + sampled = shots(unique_state_probs, options.nshots) + # move measurements dimension to the front, getting ready for extraction - measurements = np.moveaxis(sampled, 1, 0) - for ro_id, ro_dim in zip(acquisitions(sequence).keys(), measurement_mapping): - meas = measurements[ro_dim, ...] + # the shape now is: (M, Nshots, S_i, H_dim) + sampled = np.moveaxis(sampled, 1, 0) + + res_dict = {} + for ro_id, inv_idx in zip(acquisitions(sequence).keys(), inverse_map): i = index(sequence.pulse_channels(ro_id)[0], hamiltonian) # states out of the qubit computational space are classified as 1 - res = np.clip(np.stack(np.unravel_index(meas, hamiltonian.dims))[i], 0, 1) + # here sampled has dimensions (M, Nshots, S_i, H_dim) + res = np.clip(np.unravel_index(sampled[inv_idx], hamiltonian.dims)[i], 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) + # res is a (Nshots, S_i, ...) array res_dict[ro_id] = res return res_dict @@ -183,7 +211,6 @@ def singleshot_results( def results( states: NDArray, - measurement_mapping: NDArray, sequence: PulseSequence, hamiltonian: HamiltonianConfig, options: ExecutionParameters, @@ -194,6 +221,12 @@ def results( result for the execution of this single sequence, thus suitable to be returned as is. """ + + # probability dimensions are: + # (S_i, M, H_dim), where + # S_i 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 probabilities = calculate_probabilities_from_density_matrix( states, range(hamiltonian.nqubits), @@ -205,24 +238,21 @@ def results( # move measurements dimension to the front, getting ready for extraction measurements = np.moveaxis(sampled, 1, 0) - sim_results = {} - if options.averaging_mode is AveragingMode.CYCLIC: - sim_results = cyclic_results( - state_probs=probabilities, - measurement_mapping=measurement_mapping, - sequence=sequence, - hamiltonian=hamiltonian, - options=options, - ) - - if options.averaging_mode is AveragingMode.SINGLESHOT: - assert options.nshots is not None - sim_results = singleshot_results( - state_probs=probabilities, - measurement_mapping=measurement_mapping, - sequence=sequence, - hamiltonian=hamiltonian, - options=options, - ) - - return sim_results + # here we move the -2 index of the probability array, hence it now becomes: + # (M, S_i, H_dim) + probabilities = np.moveaxis(probabilities, -2, 0) + + results = ( + singleshot_results + if options.averaging_mode is AveragingMode.SINGLESHOT + else cyclic_results + ) + assert (options.averaging_mode is not AveragingMode.SINGLESHOT) or ( + options.nshots is not None + ), "nshots must be specified for SINGLESHOT mode" + return results( + state_probs=probabilities, + sequence=sequence, + hamiltonian=hamiltonian, + options=options, + ) diff --git a/tests/instruments/emulator/platforms/fixed-frequency-qutrits/calibration.json b/tests/instruments/emulator/platforms/fixed-frequency-qutrits/calibration.json new file mode 100644 index 0000000000..96223d17f5 --- /dev/null +++ b/tests/instruments/emulator/platforms/fixed-frequency-qutrits/calibration.json @@ -0,0 +1,75 @@ +{ + "single_qubits": { + "0": { + "resonator": { + "bare_frequency": 0.0, + "dressed_frequency": 5500000000.0, + "depletion_time": 0, + "bare_frequency_amplitude": null + }, + "qubit": { + "frequency_01": 5000000000.0, + "frequency_12": 4700000000.0, + "maximum_frequency": 5000000000.0, + "asymmetry": 0.0, + "sweetspot": 0.0, + "flux_coefficients": null + }, + "readout": { + "fidelity": 0.0, + "coupling": null, + "effective_temperature": null, + "ground_state": [ + 0.0, + 1.0 + ], + "excited_state": [ + 1.0, + 0.0 + ], + "qudits_frequency": {} + }, + "t1": null, + "t2": null, + "t2_spin_echo": null, + "rb_fidelity": null + }, + "1": { + "resonator": { + "bare_frequency": 0.0, + "dressed_frequency": 5500000000.0, + "depletion_time": 0, + "bare_frequency_amplitude": null + }, + "qubit": { + "frequency_01": 5000000000.0, + "frequency_12": 4700000000.0, + "maximum_frequency": 5000000000.0, + "asymmetry": 0.0, + "sweetspot": 0.0, + "flux_coefficients": null + }, + "readout": { + "fidelity": 0.0, + "coupling": null, + "effective_temperature": null, + "ground_state": [ + 0.0, + 1.0 + ], + "excited_state": [ + 1.0, + 0.0 + ], + "qudits_frequency": {} + }, + "t1": null, + "t2": null, + "t2_spin_echo": null, + "rb_fidelity": null + } + }, + "two_qubits": {}, + "readout_mitigation_matrix": null, + "flux_crosstalk_matrix": null +} From 71fba76b5037210fb82b53bedb9641c8a90c7b67 Mon Sep 17 00:00:00 2001 From: lballerio Date: Fri, 27 Mar 2026 13:47:37 +0000 Subject: [PATCH 05/18] adding doc for qibolab emulato- focusing on emulator limitations on mid-circuits measurement simulations --- doc/source/main-documentation/emulator.rst | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) create mode 100644 doc/source/main-documentation/emulator.rst diff --git a/doc/source/main-documentation/emulator.rst b/doc/source/main-documentation/emulator.rst new file mode 100644 index 0000000000..f7a7c9a195 --- /dev/null +++ b/doc/source/main-documentation/emulator.rst @@ -0,0 +1,22 @@ +.. admonition:: Work in progress + +Emulator +========= + +Qibolab contains its own simulation instrument, using which it is possible to simulate different chips configurations (one or multiple qubits, fixed or tunable frequency, with or without couplers) and run virtual Qibolab or Qibocal experiment in the same fashion as for real QPUs. + +For a more detailed discussion on how Qibolab process different QPUs (here called :class:`.Platforms`), we recommend this page :ref:`Platform guidelines `. + + +Usage +----- + +The emulator exploits a third party engine which solves the Master Equation for the sum of a constant Hamiltonian plus the time-dependent Hamiltonian, which is the Pulse Hamiltonian. From the solver, the emulator takes the solution's density matrix of the system for each timestep and selects the only ones corresponding to the acquisition pulses in the pulse sequence. Since the initial data are density matrices, the emulator does not simulate I-Q measurement such as for a real system, but simply determines the probabilities :math:`p_m=|m>`. + +That's why, even though the emulator can simulate all kind of experiment, when simulate signal experiment (i.e. with :paramref:`AcquisitionType.INTEGRATION` selected) the signal magnitude is simply the computed probabilities while the signal phase is simply meaningless. + +The emulator simulates both :paramref:`AveragingModeSINGLESHOT` (i.e. the simulator returns a finite number of shots corresponding to the measured state) and :paramref:`AveragingMode.CYCLIC` (i.e. the simulator returns the probability of the :math:`|1>` state, from the density matrix), but even though `SINGLESHOT` experiments simulate finite-shots noise it is more convienient to directly return the diagonal entries of the solution density matrix without sampling. :paramref:`AveragingMode.CYCLIC` experiments moreover simulate gaussian noise with an hard-coded sigma value. + +In order to resolve every qubit oscillations and not miss any contribution (i.e. pulse) in the total Hamiltonian, we define the NYQUIST frequency to use, which by default is set to :math:`f_N = 20 GHz`, which will allow us to correctly resolve oscillation at most of :math:`10-15 GHz`. Setting the NYQUIST frequency is a crucial part of a correct integration since from that we can adaptively tune specific options of the ODE solver and correctly solve the state evolution. For better understanding of this tuning process please see engines' implementations. + +At the time being state collapse is not implemented in QiboLab Emulator, hence this version should not be used for simulating mid-circuit measurement, but only for circuits with synchronous measurement for all qubits at the end of the circuit. After a mid-circuit measurement, the measured qubit collapses to an eigenstate; if it is entangled with other qubits, this collapse induces correlations that condition the state of the remaining system. In the present implementation, the emulator does not account for these measurement-induced correlations and therefore produces intrinsically inaccurate results in such cases. Another limitation that prevents the use of this emulator for mid-circuit measurements is that, in certain experiments, the time ordering of measurements may vary across parameter sweeps, while the :paramref:`PulseSequence`` object (i.e., the pulse sequence being simulated) remains unchanged and therefore does not capture such reordering. This discrepancy can lead to inconsistencies when executing :func:`qibolab.src._core.instruments.emulator.results.results`, which iterates simultaneously over both the time-ordering array and the acquisition pulses defined in the :paramref:`PulseSequence`. Correct behavior relies on these two iterables remaining aligned; if either is reordered during the sweep, acquisition pulses may be associated with incorrect simulation timesteps, ultimately producing invalid results. From 5f17ea4489bd319c6a53715637cc6f4a24b602e5 Mon Sep 17 00:00:00 2001 From: lballerio Date: Mon, 30 Mar 2026 11:04:08 +0000 Subject: [PATCH 06/18] apply changes --- doc/source/main-documentation/emulator.rst | 4 +- .../_core/instruments/emulator/results.py | 44 ++++++++++++------- tests/instruments/emulator/test_sequence.py | 5 ++- 3 files changed, 32 insertions(+), 21 deletions(-) diff --git a/doc/source/main-documentation/emulator.rst b/doc/source/main-documentation/emulator.rst index f7a7c9a195..3c49fba5d8 100644 --- a/doc/source/main-documentation/emulator.rst +++ b/doc/source/main-documentation/emulator.rst @@ -15,8 +15,8 @@ The emulator exploits a third party engine which solves the Master Equation for That's why, even though the emulator can simulate all kind of experiment, when simulate signal experiment (i.e. with :paramref:`AcquisitionType.INTEGRATION` selected) the signal magnitude is simply the computed probabilities while the signal phase is simply meaningless. -The emulator simulates both :paramref:`AveragingModeSINGLESHOT` (i.e. the simulator returns a finite number of shots corresponding to the measured state) and :paramref:`AveragingMode.CYCLIC` (i.e. the simulator returns the probability of the :math:`|1>` state, from the density matrix), but even though `SINGLESHOT` experiments simulate finite-shots noise it is more convienient to directly return the diagonal entries of the solution density matrix without sampling. :paramref:`AveragingMode.CYCLIC` experiments moreover simulate gaussian noise with an hard-coded sigma value. +The emulator simulates both :paramref:`AveragingMode.SINGLESHOT` (i.e. the simulator returns a finite number of shots corresponding to the measured state) and :paramref:`AveragingMode.CYCLIC` (i.e. the simulator returns the probability of the :math:`|1>` state, from the density matrix), but even though `SINGLESHOT` experiments simulate finite-shots noise it is more convenient to directly return the diagonal entries of the solution density matrix without sampling. :paramref:`AveragingMode.CYCLIC` experiments moreover simulate gaussian noise with an hard-coded sigma value. In order to resolve every qubit oscillations and not miss any contribution (i.e. pulse) in the total Hamiltonian, we define the NYQUIST frequency to use, which by default is set to :math:`f_N = 20 GHz`, which will allow us to correctly resolve oscillation at most of :math:`10-15 GHz`. Setting the NYQUIST frequency is a crucial part of a correct integration since from that we can adaptively tune specific options of the ODE solver and correctly solve the state evolution. For better understanding of this tuning process please see engines' implementations. -At the time being state collapse is not implemented in QiboLab Emulator, hence this version should not be used for simulating mid-circuit measurement, but only for circuits with synchronous measurement for all qubits at the end of the circuit. After a mid-circuit measurement, the measured qubit collapses to an eigenstate; if it is entangled with other qubits, this collapse induces correlations that condition the state of the remaining system. In the present implementation, the emulator does not account for these measurement-induced correlations and therefore produces intrinsically inaccurate results in such cases. Another limitation that prevents the use of this emulator for mid-circuit measurements is that, in certain experiments, the time ordering of measurements may vary across parameter sweeps, while the :paramref:`PulseSequence`` object (i.e., the pulse sequence being simulated) remains unchanged and therefore does not capture such reordering. This discrepancy can lead to inconsistencies when executing :func:`qibolab.src._core.instruments.emulator.results.results`, which iterates simultaneously over both the time-ordering array and the acquisition pulses defined in the :paramref:`PulseSequence`. Correct behavior relies on these two iterables remaining aligned; if either is reordered during the sweep, acquisition pulses may be associated with incorrect simulation timesteps, ultimately producing invalid results. +At the time being state collapse is not implemented in QiboLab Emulator, hence this version should not be used for simulating mid-circuit measurement, but only for circuits with synchronous measurement for all qubits at the end of the circuit. After a mid-circuit measurement, the measured qubit collapses to an eigenstate; if it is entangled with other qubits, this collapse induces correlations that condition the state of the remaining system. In the present implementation, the emulator does not account for these measurement-induced correlations and therefore produces intrinsically inaccurate results in such cases. Another limitation that prevents the use of this emulator for mid-circuit measurements is that, in certain experiments, the time ordering of measurements may vary across parameter sweeps, while the :paramref:`PulseSequence`` object (i.e., the pulse sequence being simulated) remains unchanged and therefore does not capture such reordering. This discrepancy can lead to inconsistencies when executing :func:`qibolab._core.instruments.emulator.results.results`, which iterates simultaneously over both the time-ordering array and the acquisition pulses defined in the :paramref:`PulseSequence`. Correct behavior relies on these two iterables remaining aligned; if either is reordered during the sweep, acquisition pulses may be associated with incorrect simulation timesteps, ultimately producing invalid results. diff --git a/src/qibolab/_core/instruments/emulator/results.py b/src/qibolab/_core/instruments/emulator/results.py index 6a530f2ede..47caba1e96 100644 --- a/src/qibolab/_core/instruments/emulator/results.py +++ b/src/qibolab/_core/instruments/emulator/results.py @@ -51,19 +51,28 @@ 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 calculate_probabilities_from_density_matrix(states: NDArray) -> NDArray: + """ + Calculate probabilities from a density matrix using diagonal elements. + This function extracts the diagonal elements of a density matrix and returns + their absolute values, which represent the probabilities of measurement outcomes. + + Examples + -------- + >>> dm = np.array([[0.9, 0.0], [0.0, 0.1]]) + >>> probs = calculate_probabilities_from_density_matrix(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.abs(diag) def acquisitions(sequence: PulseSequence) -> dict[PulseId, float]: @@ -99,12 +108,12 @@ def select_acquisitions( and maps them to unique acquisition values. It uses binary search to find the nearest state index for each acquisition time. - It returns a tuple containing 1 NumPy array containing the full density matrices - of the selected quantum states. + 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( @@ -247,9 +256,10 @@ def results( if options.averaging_mode is AveragingMode.SINGLESHOT else cyclic_results ) - assert (options.averaging_mode is not AveragingMode.SINGLESHOT) or ( - options.nshots is not None - ), "nshots must be specified for SINGLESHOT mode" + + if options.averaging_mode is AveragingMode.SINGLESHOT and options.nshots is None: + raise ValueError("nshots must be specified for SINGLESHOT mode") + return results( state_probs=probabilities, sequence=sequence, diff --git a/tests/instruments/emulator/test_sequence.py b/tests/instruments/emulator/test_sequence.py index 37b42b6297..b31c3e89c0 100644 --- a/tests/instruments/emulator/test_sequence.py +++ b/tests/instruments/emulator/test_sequence.py @@ -134,8 +134,9 @@ def test_cnot_sequence(platform: Platform, setup: str): ) -@pytest.mark.skip( - "The fidelity for the test is not good, either a problem of calibration or problem with emulator." +@pytest.mark.xfail( + reason="Known CZ fidelity issue (calibration/emulator mismatch).", + strict=False, ) def test_cz_sequence( platform: Platform, From 6550ed847d22d7edb1c85da7301a4267b2327ef4 Mon Sep 17 00:00:00 2001 From: lballerio Date: Tue, 31 Mar 2026 13:55:46 +0000 Subject: [PATCH 07/18] applying requested fixes --- src/qibolab/_core/instruments/emulator/results.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/src/qibolab/_core/instruments/emulator/results.py b/src/qibolab/_core/instruments/emulator/results.py index 47caba1e96..c18287983b 100644 --- a/src/qibolab/_core/instruments/emulator/results.py +++ b/src/qibolab/_core/instruments/emulator/results.py @@ -145,16 +145,13 @@ def cyclic_results( for meas_ro, ro_id in zip(state_probs, acquisitions(sequence).keys()): i = index(sequence.pulse_channels(ro_id)[0], hamiltonian) - res = np.sum(meas_ro[..., states_computational_idx[i] > 0], axis=-1) - res = np.random.normal(res, scale=0.001) + res = np.sum(meas_ro[..., states_computational_idx[i] >= 1], axis=-1) 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) - if options.acquisition_type is AcquisitionType.DISCRIMINATION: - res = np.clip(res, 0, 1) - # res is a (S_i, ...) array res_dict[ro_id] = res From 685786aec8eeed7acb0377a8c9645e41e640509b Mon Sep 17 00:00:00 2001 From: Lorenzo Ballerio <130075247+lballerio@users.noreply.github.com> Date: Wed, 1 Apr 2026 17:23:58 +0400 Subject: [PATCH 08/18] Update tests/instruments/emulator/test_sequence.py Co-authored-by: Stefano Carrazza --- tests/instruments/emulator/test_sequence.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/instruments/emulator/test_sequence.py b/tests/instruments/emulator/test_sequence.py index b31c3e89c0..fe152f9caa 100644 --- a/tests/instruments/emulator/test_sequence.py +++ b/tests/instruments/emulator/test_sequence.py @@ -136,7 +136,7 @@ def test_cnot_sequence(platform: Platform, setup: str): @pytest.mark.xfail( reason="Known CZ fidelity issue (calibration/emulator mismatch).", - strict=False, + strict=True, ) def test_cz_sequence( platform: Platform, From 99ec3bf2c13cc7302e94c76b882ef3144e1c6dba Mon Sep 17 00:00:00 2001 From: Lorenzo Ballerio <130075247+lballerio@users.noreply.github.com> Date: Wed, 1 Apr 2026 17:25:37 +0400 Subject: [PATCH 09/18] Update doc/source/main-documentation/emulator.rst Co-authored-by: Stefano Carrazza --- doc/source/main-documentation/emulator.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/doc/source/main-documentation/emulator.rst b/doc/source/main-documentation/emulator.rst index 3c49fba5d8..0246bb6072 100644 --- a/doc/source/main-documentation/emulator.rst +++ b/doc/source/main-documentation/emulator.rst @@ -1,5 +1,7 @@ .. admonition:: Work in progress + This documentation is draft and may be incomplete. + Emulator ========= From 929f92f2f249b8172737892455db52dd711e4445 Mon Sep 17 00:00:00 2001 From: lballerio Date: Wed, 1 Apr 2026 21:09:48 +0400 Subject: [PATCH 10/18] vectorizing _cyclic_results and _singleshot_results functions --- .../_core/instruments/emulator/results.py | 115 +++++++++--------- 1 file changed, 60 insertions(+), 55 deletions(-) diff --git a/src/qibolab/_core/instruments/emulator/results.py b/src/qibolab/_core/instruments/emulator/results.py index c18287983b..be7281f3d3 100644 --- a/src/qibolab/_core/instruments/emulator/results.py +++ b/src/qibolab/_core/instruments/emulator/results.py @@ -51,16 +51,19 @@ def shots(probabilities: NDArray, nshots: int) -> NDArray: return np.moveaxis(shots, -1, 0) -def calculate_probabilities_from_density_matrix(states: NDArray) -> NDArray: +def _extract_probabilities(states: NDArray) -> NDArray: """ Calculate probabilities from a density matrix using diagonal elements. - This function extracts the diagonal elements of a density matrix and returns - their absolute values, which represent the probabilities of measurement outcomes. + + 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 = calculate_probabilities_from_density_matrix(dm) + >>> probs = _extract_probabilities(dm) >>> probs array([0.9, 0.1]) """ @@ -72,7 +75,7 @@ def calculate_probabilities_from_density_matrix(states: NDArray) -> NDArray: np.array([...] + [0, 0]), np.array([...] + [0]), ) - return np.abs(diag) + return np.clip(diag.real, 0, 1) def acquisitions(sequence: PulseSequence) -> dict[PulseId, float]: @@ -116,7 +119,7 @@ def select_acquisitions( return np.stack([states[n].full() for n in samples])[index_pos] -def cyclic_results( +def _cyclic_results( state_probs: NDArray, sequence: PulseSequence, hamiltonian: HamiltonianConfig, @@ -126,39 +129,40 @@ def cyclic_results( excited state population. Computes readout results by projecting quantum state probabilities onto measurement subspaces and applying configured post-processing. - - Notes: - - States outside the computational subspace (values > 1) are classified as 1. - - For integration acquisition type, imaginary components are set to zero. """ # Through the entire function state_probs has dimensions: - # (M, S_i, H_dim), where + # (*S, M *H_dim), where + # *S is the number of iteration for each sweep in the experiment # M is the number of measurements applied in the pulse sequence - # S_i is the number of iteration for each sweep in the experiment - # H_dim is the complete system dimension + # *H_dim is the complete system dimension states_computational_idx = np.stack( np.unravel_index(np.arange(state_probs.shape[-1]), hamiltonian.dims) ) - res_dict = {} - for meas_ro, ro_id in zip(state_probs, acquisitions(sequence).keys()): - i = index(sequence.pulse_channels(ro_id)[0], hamiltonian) + 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] - res = np.sum(meas_ro[..., states_computational_idx[i] >= 1], axis=-1) + # 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 - 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) + # res is a (M, *S, ...) array + res = np.moveaxis(np.sum(np.where(mask, state_probs, 0), axis=-1), -1, 0) - # res is a (S_i, ...) array - res_dict[ro_id] = res + 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 res_dict + return dict(zip(acq_id, res)) -def singleshot_results( +def _singleshot_results( state_probs: NDArray, sequence: PulseSequence, hamiltonian: HamiltonianConfig, @@ -168,10 +172,6 @@ def singleshot_results( 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. - - Notes: - - States outside the computational subspace (values > 1) are classified as 1. - - For integration acquisition type, imaginary components are set to zero. """ # select only unique times of measurements @@ -181,38 +181,47 @@ def singleshot_results( 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, S_i, H_dim), where + # (Nshots, M, *S, *H_dim), where # Nshots is simply the number of shots we average on # M is the number of measurements applied in the pulse sequence - # S_i is the number of iteration for each sweep in the experiment - # H_dim is the complete system dimension + # *S is the number of iteration for each sweep in the experiment + # *H_dim is the complete system dimension sampled = shots(unique_state_probs, options.nshots) # move measurements dimension to the front, getting ready for extraction - # the shape now is: (M, Nshots, S_i, H_dim) + # the shape now is: (M, Nshots, *S, *H_dim) sampled = np.moveaxis(sampled, 1, 0) - res_dict = {} - for ro_id, inv_idx in zip(acquisitions(sequence).keys(), inverse_map): - i = index(sequence.pulse_channels(ro_id)[0], hamiltonian) - # states out of the qubit computational space are classified as 1 - # here sampled has dimensions (M, Nshots, S_i, H_dim) - res = np.clip(np.unravel_index(sampled[inv_idx], hamiltonian.dims)[i], 0, 1) + 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 + ] - 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) + # res is a (M, M_unique, 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) - # res is a (Nshots, S_i, ...) array - res_dict[ro_id] = res + 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 res_dict + return dict(zip(acq_id, res)) def results( @@ -229,11 +238,11 @@ def results( """ # probability dimensions are: - # (S_i, M, H_dim), where - # S_i is the number of iteration for each sweep in the experiment + # (*S, M, *H_dim), 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 - probabilities = calculate_probabilities_from_density_matrix( + # *H_dim is the complete system dimension + probabilities = _extract_probabilities( states, range(hamiltonian.nqubits), hamiltonian.nqubits, @@ -244,14 +253,10 @@ def results( # move measurements dimension to the front, getting ready for extraction measurements = np.moveaxis(sampled, 1, 0) - # here we move the -2 index of the probability array, hence it now becomes: - # (M, S_i, H_dim) - probabilities = np.moveaxis(probabilities, -2, 0) - results = ( - singleshot_results + _singleshot_results if options.averaging_mode is AveragingMode.SINGLESHOT - else cyclic_results + else _cyclic_results ) if options.averaging_mode is AveragingMode.SINGLESHOT and options.nshots is None: From 40bfc0b6e52071d7706e1b8510be51e4fcbd0140 Mon Sep 17 00:00:00 2001 From: lballerio Date: Wed, 1 Apr 2026 21:19:54 +0400 Subject: [PATCH 11/18] restored test_cz_sequence for emulator test --- tests/instruments/emulator/test_sequence.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/tests/instruments/emulator/test_sequence.py b/tests/instruments/emulator/test_sequence.py index fe152f9caa..3b7695549c 100644 --- a/tests/instruments/emulator/test_sequence.py +++ b/tests/instruments/emulator/test_sequence.py @@ -134,10 +134,6 @@ def test_cnot_sequence(platform: Platform, setup: str): ) -@pytest.mark.xfail( - reason="Known CZ fidelity issue (calibration/emulator mismatch).", - strict=True, -) def test_cz_sequence( platform: Platform, ): From f276b72321fdfa324b7aa985cf4dd0b30e4b9694 Mon Sep 17 00:00:00 2001 From: lballerio Date: Wed, 1 Apr 2026 21:31:53 +0400 Subject: [PATCH 12/18] renaming of two functions in emulator.py --- src/qibolab/_core/instruments/emulator/emulator.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/qibolab/_core/instruments/emulator/emulator.py b/src/qibolab/_core/instruments/emulator/emulator.py index d29287a6e2..0456d009e0 100644 --- a/src/qibolab/_core/instruments/emulator/emulator.py +++ b/src/qibolab/_core/instruments/emulator/emulator.py @@ -69,13 +69,13 @@ def play( sequences_ = (seq.align_to_delays() for seq in sequences) results_to_process = ( - self._results(configs, sequence, options, sweepers) + self._play_sequence(configs, sequence, options, sweepers) for sequence in sequences_ ) return reduce(or_, results_to_process) - def _results( + def _play_sequence( self, configs: dict[str, Config], sequence: PulseSequence, @@ -116,7 +116,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 @@ -138,10 +138,10 @@ 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: - """Execute a pulse sequence on the quantum emulator. + """Evolve a pulse sequence on the quantum emulator. This method updates the sequence parameters, generates the time grid, constructs the time-dependent Hamiltonian, evolves the initial state with optional collapse From 0b2e9602ea2f8c8e172272fdb8a9ca5967d1d6c5 Mon Sep 17 00:00:00 2001 From: lballerio Date: Wed, 1 Apr 2026 21:47:09 +0400 Subject: [PATCH 13/18] rewriting emulator doc --- doc/source/main-documentation/emulator.rst | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/doc/source/main-documentation/emulator.rst b/doc/source/main-documentation/emulator.rst index 0246bb6072..650706093c 100644 --- a/doc/source/main-documentation/emulator.rst +++ b/doc/source/main-documentation/emulator.rst @@ -1,24 +1,25 @@ .. admonition:: Work in progress - This documentation is draft and may be incomplete. +This documentation is currently in draft form and may be incomplete. Emulator -========= +======== -Qibolab contains its own simulation instrument, using which it is possible to simulate different chips configurations (one or multiple qubits, fixed or tunable frequency, with or without couplers) and run virtual Qibolab or Qibocal experiment in the same fashion as for real QPUs. - -For a more detailed discussion on how Qibolab process different QPUs (here called :class:`.Platforms`), we recommend this page :ref:`Platform guidelines `. +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 or Qibocal 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 . Usage ----- -The emulator exploits a third party engine which solves the Master Equation for the sum of a constant Hamiltonian plus the time-dependent Hamiltonian, which is the Pulse Hamiltonian. From the solver, the emulator takes the solution's density matrix of the system for each timestep and selects the only ones corresponding to the acquisition pulses in the pulse sequence. Since the initial data are density matrices, the emulator does not simulate I-Q measurement such as for a real system, but simply determines the probabilities :math:`p_m=|m>`. +The emulator relies on 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 = \langle m | \rho | m \rangle` for each computational basis state :math:`|m\rangle`. 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. -That's why, even though the emulator can simulate all kind of experiment, when simulate signal experiment (i.e. with :paramref:`AcquisitionType.INTEGRATION` selected) the signal magnitude is simply the computed probabilities while the signal phase is simply meaningless. +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:`|1\rangle` 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. -The emulator simulates both :paramref:`AveragingMode.SINGLESHOT` (i.e. the simulator returns a finite number of shots corresponding to the measured state) and :paramref:`AveragingMode.CYCLIC` (i.e. the simulator returns the probability of the :math:`|1>` state, from the density matrix), but even though `SINGLESHOT` experiments simulate finite-shots noise it is more convenient to directly return the diagonal entries of the solution density matrix without sampling. :paramref:`AveragingMode.CYCLIC` experiments moreover simulate gaussian noise with an hard-coded sigma value. +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. -In order to resolve every qubit oscillations and not miss any contribution (i.e. pulse) in the total Hamiltonian, we define the NYQUIST frequency to use, which by default is set to :math:`f_N = 20 GHz`, which will allow us to correctly resolve oscillation at most of :math:`10-15 GHz`. Setting the NYQUIST frequency is a crucial part of a correct integration since from that we can adaptively tune specific options of the ODE solver and correctly solve the state evolution. For better understanding of this tuning process please see engines' implementations. +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. -At the time being state collapse is not implemented in QiboLab Emulator, hence this version should not be used for simulating mid-circuit measurement, but only for circuits with synchronous measurement for all qubits at the end of the circuit. After a mid-circuit measurement, the measured qubit collapses to an eigenstate; if it is entangled with other qubits, this collapse induces correlations that condition the state of the remaining system. In the present implementation, the emulator does not account for these measurement-induced correlations and therefore produces intrinsically inaccurate results in such cases. Another limitation that prevents the use of this emulator for mid-circuit measurements is that, in certain experiments, the time ordering of measurements may vary across parameter sweeps, while the :paramref:`PulseSequence`` object (i.e., the pulse sequence being simulated) remains unchanged and therefore does not capture such reordering. This discrepancy can lead to inconsistencies when executing :func:`qibolab._core.instruments.emulator.results.results`, which iterates simultaneously over both the time-ordering array and the acquisition pulses defined in the :paramref:`PulseSequence`. Correct behavior relies on these two iterables remaining aligned; if either is reordered during the sweep, acquisition pulses may be associated with incorrect simulation timesteps, ultimately producing invalid results. +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. From 4080923e83ee74c78b532d402eaab27e56521ec2 Mon Sep 17 00:00:00 2001 From: lballerio Date: Thu, 2 Apr 2026 09:55:17 +0400 Subject: [PATCH 14/18] changing some comments in results.py file or the emulator --- src/qibolab/_core/instruments/emulator/results.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/qibolab/_core/instruments/emulator/results.py b/src/qibolab/_core/instruments/emulator/results.py index be7281f3d3..e0bb430101 100644 --- a/src/qibolab/_core/instruments/emulator/results.py +++ b/src/qibolab/_core/instruments/emulator/results.py @@ -191,15 +191,15 @@ def _singleshot_results( unique_state_probs = state_probs[direct_map] # shots function returns a vector of shape: - # (Nshots, M, *S, *H_dim), where + # (Nshots, M_unique, *S, *H_dim), where # Nshots is simply the number of shots we average on - # M is the number of measurements applied in the pulse sequence + # M_unique is the number of unique measurement times # *S is the number of iteration for each sweep in the experiment # *H_dim is the complete system dimension sampled = shots(unique_state_probs, options.nshots) # move measurements dimension to the front, getting ready for extraction - # the shape now is: (M, Nshots, *S, *H_dim) + # the shape now is: (M_unique, Nshots, *S, *H_dim) sampled = np.moveaxis(sampled, 1, 0) acq_id = acquisitions(sequence).keys() @@ -209,7 +209,8 @@ def _singleshot_results( index(sequence.pulse_channels(ro_id)[0], hamiltonian) for ro_id in acq_id ] - # res is a (M, M_unique, Nshots, *S, ...) array + # 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 ] From 8350b6fbb80161512fc93d20f21ec6740ebbe03b Mon Sep 17 00:00:00 2001 From: Lorenzo Ballerio <130075247+lballerio@users.noreply.github.com> Date: Tue, 14 Apr 2026 11:21:10 +0400 Subject: [PATCH 15/18] Update doc/source/main-documentation/emulator.rst Co-authored-by: Alessandro Candido --- doc/source/main-documentation/emulator.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/source/main-documentation/emulator.rst b/doc/source/main-documentation/emulator.rst index 650706093c..f0213a580f 100644 --- a/doc/source/main-documentation/emulator.rst +++ b/doc/source/main-documentation/emulator.rst @@ -5,7 +5,7 @@ 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 or Qibocal experiments in a manner fully consistent with their execution on real quantum processing units (QPUs). +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 . From ec140bb81f10c773e0bf8f1183cb6472e812a486 Mon Sep 17 00:00:00 2001 From: lballerio Date: Fri, 3 Apr 2026 14:07:27 +0400 Subject: [PATCH 16/18] fixing wrong equation in emulatgor.rst documentation --- doc/source/main-documentation/emulator.rst | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/doc/source/main-documentation/emulator.rst b/doc/source/main-documentation/emulator.rst index f0213a580f..1c73ec5886 100644 --- a/doc/source/main-documentation/emulator.rst +++ b/doc/source/main-documentation/emulator.rst @@ -7,16 +7,16 @@ 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 . +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 `. Usage ----- -The emulator relies on 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. +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 = \langle m | \rho | m \rangle` for each computational basis state :math:`|m\rangle`. 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. +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:`|1\rangle` 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. +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. From 77f9ea2fa104b0198c2307d89e62b7e3b347e17b Mon Sep 17 00:00:00 2001 From: lballerio Date: Tue, 14 Apr 2026 11:35:23 +0400 Subject: [PATCH 17/18] addressing @alecandido final comments in PR --- .../_core/instruments/emulator/emulator.py | 9 +++++- .../_core/instruments/emulator/results.py | 32 +++++++++---------- 2 files changed, 24 insertions(+), 17 deletions(-) diff --git a/src/qibolab/_core/instruments/emulator/emulator.py b/src/qibolab/_core/instruments/emulator/emulator.py index 0456d009e0..cfef30d6d3 100644 --- a/src/qibolab/_core/instruments/emulator/emulator.py +++ b/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,6 +65,13 @@ 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) diff --git a/src/qibolab/_core/instruments/emulator/results.py b/src/qibolab/_core/instruments/emulator/results.py index e0bb430101..05f8ebe237 100644 --- a/src/qibolab/_core/instruments/emulator/results.py +++ b/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 @@ -132,10 +145,7 @@ def _cyclic_results( """ # Through the entire function state_probs has dimensions: - # (*S, M *H_dim), 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 + # (*S, M *H_dim) states_computational_idx = np.stack( np.unravel_index(np.arange(state_probs.shape[-1]), hamiltonian.dims) ) @@ -191,11 +201,7 @@ def _singleshot_results( unique_state_probs = state_probs[direct_map] # shots function returns a vector of shape: - # (Nshots, M_unique, *S, *H_dim), where - # Nshots is simply the number of shots we average on - # M_unique is the number of unique measurement times - # *S is the number of iteration for each sweep in the experiment - # *H_dim is the complete system dimension + # (Nshots, M_unique, *S, *H_dim) sampled = shots(unique_state_probs, options.nshots) # move measurements dimension to the front, getting ready for extraction @@ -239,10 +245,7 @@ def results( """ # probability dimensions are: - # (*S, M, *H_dim), 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 + # (*S, M, *H_dim) probabilities = _extract_probabilities( states, range(hamiltonian.nqubits), @@ -260,9 +263,6 @@ def results( else _cyclic_results ) - if options.averaging_mode is AveragingMode.SINGLESHOT and options.nshots is None: - raise ValueError("nshots must be specified for SINGLESHOT mode") - return results( state_probs=probabilities, sequence=sequence, From 086b06d5ae6b4bc06087538c6f94d4ae16ce4daf Mon Sep 17 00:00:00 2001 From: lballerio Date: Thu, 23 Apr 2026 13:54:39 +0400 Subject: [PATCH 18/18] untracked rebase modifications --- src/qibolab/_core/instruments/emulator/results.py | 13 ++----------- 1 file changed, 2 insertions(+), 11 deletions(-) diff --git a/src/qibolab/_core/instruments/emulator/results.py b/src/qibolab/_core/instruments/emulator/results.py index 05f8ebe237..56bc803d94 100644 --- a/src/qibolab/_core/instruments/emulator/results.py +++ b/src/qibolab/_core/instruments/emulator/results.py @@ -21,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 @@ -246,16 +246,7 @@ def results( # probability dimensions are: # (*S, M, *H_dim) - probabilities = _extract_probabilities( - 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) + probabilities = _extract_probabilities(states) results = ( _singleshot_results