From 670b22daa909fe14ae39835340c69da8ec9cf2e7 Mon Sep 17 00:00:00 2001 From: Dennis Wayo <117969019+DennisWayo@users.noreply.github.com> Date: Tue, 23 Jun 2026 15:01:16 +0500 Subject: [PATCH 1/3] Bypass cyclic result state indexing --- .../_core/instruments/emulator/results.py | 33 ++++---- tests/instruments/emulator/test_results.py | 75 +++++++++++++++++++ 2 files changed, 95 insertions(+), 13 deletions(-) diff --git a/src/qibolab/_core/instruments/emulator/results.py b/src/qibolab/_core/instruments/emulator/results.py index 4ff7ce07c6..6c2f1cfa19 100644 --- a/src/qibolab/_core/instruments/emulator/results.py +++ b/src/qibolab/_core/instruments/emulator/results.py @@ -112,6 +112,16 @@ def index(ch: ChannelId, hconfig: HamiltonianConfig) -> int: return hconfig.hilbert_space_index(target) +def _marginalize_probability( + probabilities: NDArray, dims: list[int], index: int +) -> NDArray: + """Marginalize full-system probabilities to a single Hilbert-space component.""" + leading = probabilities.shape[:-1] + probabilities = probabilities.reshape(*leading, *dims) + axes = tuple(len(leading) + i for i in range(len(dims)) if i != index) + return probabilities.sum(axis=axes) + + def select_acquisitions( states: list[Operator], acquisitions: Iterable[float], times: NDArray ) -> NDArray: @@ -142,25 +152,22 @@ def _cyclic_results( 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() + acq_id = list(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. + # corresponding Hilbert-space index to marginalize. 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) + res = np.stack( + [ + _marginalize_probability( + state_probs[..., measurement, :], hamiltonian.dims, qubit_index + )[..., 1:].sum(axis=-1) + for measurement, qubit_index in enumerate(qubit_indices) + ] + ) if options.acquisition_type is AcquisitionType.INTEGRATION: res = np.random.normal(res, scale=0.001) diff --git a/tests/instruments/emulator/test_results.py b/tests/instruments/emulator/test_results.py index 8e564f2ec0..df8b498ad6 100644 --- a/tests/instruments/emulator/test_results.py +++ b/tests/instruments/emulator/test_results.py @@ -1,5 +1,16 @@ import numpy as np +from qibolab._core.execution_parameters import ( + AcquisitionType, + AveragingMode, + ExecutionParameters, +) +from qibolab._core.instruments.emulator.results import ( + _cyclic_results, + _marginalize_probability, +) +from qibolab._core.pulses import Acquisition + def random_states(space: tuple[int, ...], sweeps: tuple[int, ...] = (), nacq: int = 1): dimension = np.prod(space, dtype=int) @@ -8,3 +19,67 @@ 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_marginalize_probability_preserves_sweep_axes(): + probabilities = np.arange(1, 25).reshape(2, 12) + + marginalized = _marginalize_probability(probabilities, [2, 3, 2], 1) + expected = probabilities.reshape(2, 2, 3, 2).sum(axis=(1, 3)) + + np.testing.assert_allclose(marginalized, expected) + + +def test_cyclic_integration_results_marginalize_probabilities(monkeypatch): + monkeypatch.setattr(np.random, "normal", lambda loc, scale: loc) + acq0 = Acquisition(duration=1) + acq1 = Acquisition(duration=1) + sequence = _Sequence( + { + "0/acquisition": [acq0], + "1/acquisition": [acq1], + } + ) + probabilities = np.array( + [ + [[0.10, 0.20, 0.30], [0.05, 0.15, 0.20]], + [[0.10, 0.20, 0.30], [0.05, 0.15, 0.20]], + ] + ) + results = _cyclic_results( + state_probs=probabilities.reshape(2, -1), + sequence=sequence, + hamiltonian=_HamiltonianConfig(dims=[2, 3]), + options=ExecutionParameters( + acquisition_type=AcquisitionType.INTEGRATION, + averaging_mode=AveragingMode.CYCLIC, + ), + ) + + np.testing.assert_allclose(results[acq0.id], [0.40, 0.0]) + np.testing.assert_allclose(results[acq1.id], [0.85, 0.0]) + + +class _Sequence: + def __init__(self, channels): + self._channels = channels + self.channels = list(channels) + self._pulse_channels = { + event.id: channel + for channel, events in channels.items() + for event in events + } + + def channel(self, channel): + return self._channels[channel] + + def pulse_channels(self, pulse_id): + return [self._pulse_channels[pulse_id]] + + +class _HamiltonianConfig: + def __init__(self, dims): + self.dims = dims + + def hilbert_space_index(self, target): + return target From 4d87a67f27fffc0c2f5ab4b65349d7a16c578c1c Mon Sep 17 00:00:00 2001 From: Dennis Wayo <117969019+DennisWayo@users.noreply.github.com> Date: Wed, 24 Jun 2026 13:14:48 +0500 Subject: [PATCH 2/3] Refine cyclic probability marginalization --- .../_core/instruments/emulator/results.py | 28 +++++++++---------- tests/instruments/emulator/test_results.py | 13 ++++++--- 2 files changed, 23 insertions(+), 18 deletions(-) diff --git a/src/qibolab/_core/instruments/emulator/results.py b/src/qibolab/_core/instruments/emulator/results.py index 6c2f1cfa19..ff71e52ebb 100644 --- a/src/qibolab/_core/instruments/emulator/results.py +++ b/src/qibolab/_core/instruments/emulator/results.py @@ -113,13 +113,20 @@ def index(ch: ChannelId, hconfig: HamiltonianConfig) -> int: def _marginalize_probability( - probabilities: NDArray, dims: list[int], index: int + probabilities: NDArray, dims: list[int], indices: Iterable[int] ) -> NDArray: - """Marginalize full-system probabilities to a single Hilbert-space component.""" - leading = probabilities.shape[:-1] - probabilities = probabilities.reshape(*leading, *dims) - axes = tuple(len(leading) + i for i in range(len(dims)) if i != index) - return probabilities.sum(axis=axes) + """Marginalize full-system probabilities to measured Hilbert-space components.""" + res = [] + for measurement, index in enumerate(indices): + distribution = probabilities[..., measurement, :] + leading = distribution.shape[:-1] + distribution = distribution.reshape(*leading, *dims) + axes = tuple(len(leading) + i for i in range(len(dims)) if i != index) + marginalized = distribution.sum(axis=axes) + # States outside the ground state are classified as 1. + res.append(marginalized[..., 1:].sum(axis=-1)) + + return np.stack(res) def select_acquisitions( @@ -160,14 +167,7 @@ def _cyclic_results( ] # res is a (M, *S, ...) array - res = np.stack( - [ - _marginalize_probability( - state_probs[..., measurement, :], hamiltonian.dims, qubit_index - )[..., 1:].sum(axis=-1) - for measurement, qubit_index in enumerate(qubit_indices) - ] - ) + res = _marginalize_probability(state_probs, hamiltonian.dims, qubit_indices) if options.acquisition_type is AcquisitionType.INTEGRATION: res = np.random.normal(res, scale=0.001) diff --git a/tests/instruments/emulator/test_results.py b/tests/instruments/emulator/test_results.py index df8b498ad6..ef95801747 100644 --- a/tests/instruments/emulator/test_results.py +++ b/tests/instruments/emulator/test_results.py @@ -22,10 +22,15 @@ def random_states(space: tuple[int, ...], sweeps: tuple[int, ...] = (), nacq: in def test_marginalize_probability_preserves_sweep_axes(): - probabilities = np.arange(1, 25).reshape(2, 12) - - marginalized = _marginalize_probability(probabilities, [2, 3, 2], 1) - expected = probabilities.reshape(2, 2, 3, 2).sum(axis=(1, 3)) + probabilities = np.arange(1, 49).reshape(2, 2, 12) + + marginalized = _marginalize_probability(probabilities, [2, 3, 2], [1, 2]) + expected = np.stack( + ( + probabilities[:, 0].reshape(2, 2, 3, 2).sum(axis=(1, 3))[..., 1:].sum(-1), + probabilities[:, 1].reshape(2, 2, 3, 2).sum(axis=(1, 2))[..., 1:].sum(-1), + ) + ) np.testing.assert_allclose(marginalized, expected) From 932ae744cd55ca855a13cd2a2003fe1e0dd84805 Mon Sep 17 00:00:00 2001 From: Dennis Wayo <117969019+DennisWayo@users.noreply.github.com> Date: Wed, 24 Jun 2026 15:16:33 +0500 Subject: [PATCH 3/3] Clarify marginalization stacking --- src/qibolab/_core/instruments/emulator/results.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/qibolab/_core/instruments/emulator/results.py b/src/qibolab/_core/instruments/emulator/results.py index ff71e52ebb..fb1f38b8f5 100644 --- a/src/qibolab/_core/instruments/emulator/results.py +++ b/src/qibolab/_core/instruments/emulator/results.py @@ -123,7 +123,8 @@ def _marginalize_probability( distribution = distribution.reshape(*leading, *dims) axes = tuple(len(leading) + i for i in range(len(dims)) if i != index) marginalized = distribution.sum(axis=axes) - # States outside the ground state are classified as 1. + # States outside the ground state are classified as 1. This sum is done + # before stacking because measured subsystems may have different dimensions. res.append(marginalized[..., 1:].sum(axis=-1)) return np.stack(res)