From 6f821a586dd6ca6e0a8782fe04be76a1df36e83d Mon Sep 17 00:00:00 2001 From: Tobias Hangleiter Date: Fri, 1 Nov 2024 08:48:21 +0100 Subject: [PATCH 01/26] Add util.adot --- filter_functions/pulse_sequence.py | 1 + 1 file changed, 1 insertion(+) diff --git a/filter_functions/pulse_sequence.py b/filter_functions/pulse_sequence.py index 6f674fd..ff59e38 100644 --- a/filter_functions/pulse_sequence.py +++ b/filter_functions/pulse_sequence.py @@ -1824,6 +1824,7 @@ def concatenate( control_matrix = numeric.calculate_control_matrix_from_atomic( phases, control_matrix_atomic, propagators_liouville, show_progressbar, which='correlations' if calc_pulse_correlation_FF else 'total', + return_accumulated=calc_second_order_FF ) # Set the attribute and calculate filter function (if the pulse correlation From 57ff3c4db96f0de997d4425d668a8cb99190b835 Mon Sep 17 00:00:00 2001 From: Tobias Hangleiter Date: Fri, 1 Nov 2024 08:49:03 +0100 Subject: [PATCH 02/26] Add first draft for second order concatenation To do: - edge cases - document - test --- filter_functions/numeric.py | 61 ++++++++++++++++++++++++++++++ filter_functions/pulse_sequence.py | 30 ++++++++++++++- tests/test_sequencing.py | 22 +++++++++++ 3 files changed, 111 insertions(+), 2 deletions(-) diff --git a/filter_functions/numeric.py b/filter_functions/numeric.py index e123516..3c2df84 100644 --- a/filter_functions/numeric.py +++ b/filter_functions/numeric.py @@ -1773,6 +1773,67 @@ def calculate_pulse_correlation_filter_function(control_matrix: ndarray, return np.einsum(subscripts, control_matrix.conj(), control_matrix) +def calculate_second_order_from_atomic( + basis: Basis, + filter_function_atomic: ndarray, + control_matrix_atomic: ndarray, + control_matrix_accumulated: ndarray, + phases: ndarray, + propagators: ndarray, + propagators_liouville: ndarray, + intermediates: Sequence[Mapping[str, ndarray]] +): + G, n_nops, n_basis, n_omega = control_matrix_atomic.shape + d = propagators.shape[-1] + + tmp_1 = np.empty((n_nops, n_basis, n_omega), dtype=complex) + tmp_2 = np.empty((n_nops, n_nops, n_basis, n_basis, n_omega), dtype=complex) + basis_transformed = np.empty((n_basis, d, d), dtype=complex) + n_opers_basis = np.empty((n_nops, n_basis, d, d), dtype=complex) + + expr_1 = oe.contract_expression('ajo,jk->ako', (n_nops, n_basis, n_omega), (n_basis, n_basis)) + expr_2 = oe.contract_expression('abpqo,pk,ql->abklo', + (n_nops, n_nops, n_basis, n_basis, n_omega), + (n_basis, n_basis), (n_basis, n_basis)) + expr_3 = oe.contract_expression('oijmn,akij,blmn->abklo', (n_omega, d, d, d, d), + (n_nops, n_basis, d, d), (n_nops, n_basis, d, d)) + + def _incomplete_time_step(h, out): + _transform_by_unitary(eigvecs_propagated[h], basis, out=basis_transformed) + np.multiply(n_opers_transformed[:, h, None], basis_transformed.swapaxes(-1, -2), + out=n_opers_basis) + return expr_3(second_order_integral[h], n_opers_basis, n_opers_basis, out=out) + + result = filter_function_atomic.copy() + + for g in range(1, G): + eigvecs_propagated = _propagate_eigenvectors(propagators[g-1:g], + intermediates[g]['eigvecs_propagated']) + n_opers_transformed = intermediates[g]['n_opers_transformed'] + second_order_integral = intermediates[g]['second_order_integral'] + second_order_complete_steps = intermediates[g]['second_order_complete_steps'] + + # B'_(g-1)(ω) + tmp_1[:] = control_matrix_accumulated[g-1] + tmp_1 *= phases[g-1].conj() + tmp_1 = expr_1(tmp_1, propagators_liouville[g-1].T, out=tmp_1) + + # B^(g)*(ω) B'_(g-1)(ω) + tmp_2 = np.multiply(control_matrix_atomic[g, :, None, :, None].conj(), + tmp_1[None, :, None, :], out=tmp_2) + + # N_(g)(ω) + tmp_2 += second_order_complete_steps + + result += expr_2(tmp_2, propagators_liouville[g-1], propagators_liouville[g-1], out=tmp_2) + + for h in range(len(eigvecs_propagated)): + # J'_(g)(ω) + result += _incomplete_time_step(h, out=tmp_2) + + return result + + def diagonalize(hamiltonian: ndarray, dt: Coefficients) -> Tuple[ndarray, ndarray, ndarray]: r"""Diagonalize a Hamiltonian. diff --git a/filter_functions/pulse_sequence.py b/filter_functions/pulse_sequence.py index ff59e38..52d674c 100644 --- a/filter_functions/pulse_sequence.py +++ b/filter_functions/pulse_sequence.py @@ -1657,6 +1657,7 @@ def concatenate( pulses: Iterable[PulseSequence], calc_pulse_correlation_FF: bool = False, calc_filter_function: Optional[bool] = None, + calc_second_order_FF: Optional[bool] = None, which: str = 'fidelity', omega: Optional[Coefficients] = None, show_progressbar: bool = False @@ -1691,7 +1692,9 @@ def concatenate( carried out or not. Overrides the automatic behavior of calculating it if at least one pulse has a cached control matrix. If ``True`` and no pulse has a cached control matrix, a - list of frequencies must be supplied as *omega*. + list of frequencies must be supplied as *omega*. Overridden if + either *calc_pulse_correlation_FF* or *calc_second_order_FF* are + true. which: str, optional Which filter function to compute. Either 'fidelity' (default) or 'generalized' (see :meth:`PulseSequence.get_filter_function` and @@ -1720,7 +1723,10 @@ def concatenate( if all(pls.is_cached('total_propagator') for pls in pulses): newpulse.total_propagator = util.mdot([pls.total_propagator for pls in pulses][::-1]) - if calc_filter_function is False and not calc_pulse_correlation_FF: + if calc_pulse_correlation_FF or calc_second_order_FF is True: + calc_filter_function = True + + if calc_filter_function is False: return newpulse # If the pulses have different noise operators, we cannot reuse cached @@ -1742,6 +1748,11 @@ def concatenate( if identifier in pulse_identifier: n_opers_present[i, j] = True + if calc_second_order_FF and not n_opers_present.all(): + warn('Second order FF requested but not all pulses have the same n_opers. ' + 'Not implemented.', UserWarning) + calc_second_order_FF = False + # If at least two pulses have the same noise operators, we gain an # advantage when concatenating the filter functions over calculating them # from scratch at a later point @@ -1827,6 +1838,21 @@ def concatenate( return_accumulated=calc_second_order_FF ) + if calc_second_order_FF: + control_matrix, control_matrix_accumulated = control_matrix + filter_function = numeric.calculate_second_order_from_atomic( + basis=newpulse.basis, + filter_function_atomic=pulses[0].get_filter_function(omega, order=2), + control_matrix_atomic=control_matrix_atomic, + control_matrix_accumulated=control_matrix_accumulated, + phases=phases, + propagators=util.adot([pulse.total_propagator for pulse in pulses[:-1]]), + propagators_liouville=propagators_liouville, + intermediates=[pulse.intermediates for pulse in pulses] + ) + + newpulse.cache_filter_function(omega, filter_function=filter_function, order=2) + # Set the attribute and calculate filter function (if the pulse correlation # FF has been calculated, this is a little overhead but negligible) newpulse.cache_filter_function(omega, control_matrix, which=which) diff --git a/tests/test_sequencing.py b/tests/test_sequencing.py index 5fd38ee..b8e6ebd 100644 --- a/tests/test_sequencing.py +++ b/tests/test_sequencing.py @@ -448,6 +448,23 @@ def test_concatenate_split_cnot(self): cnot_concatenated.frequency_data['filter_function'], rtol, atol) + def test_second_order(self): + for d, n_dt in zip(rng.integers(2, 5, 10), rng.integers(1, 11, 10)): + pulse = testutil.rand_pulse_sequence(d, n_dt) + omega = util.get_sample_frequencies(pulse, 11) + pulses = list(pulse) + pulse.cache_filter_function(omega, order=1, cache_intermediates=True) + pulse.cache_filter_function(omega, order=2, cache_intermediates=True) + for pls in pulses: + pls.cache_filter_function(omega, order=1, cache_intermediates=True) + pls.cache_filter_function(omega, order=2, cache_intermediates=True) + + concat_pulse = ff.concatenate(pulses, calc_second_order_FF=True) + self.assertArrayAlmostEqual(concat_pulse.get_filter_function(omega, order=2), + pulse.get_filter_function(omega, order=2), + atol=1e-13) + + def test_different_n_opers(self): """Test behavior when concatenating with different n_opers.""" for d, n_dt in zip(rng.integers(2, 5, 20), @@ -541,6 +558,11 @@ def test_different_n_opers(self): # Test forcibly caching self.assertTrue(pulse_13_2.is_cached('filter_function')) + with self.assertWarns(UserWarning): + # Issues a warning and disables second order calculation if unequal nopers + result = ff.concatenate([pulse_1, pulse_3], calc_second_order_FF=True) + self.assertFalse(result.is_cached('filter_function_2')) + def test_concatenate_periodic(self): """Test concatenation for periodic Hamiltonians""" X, Y, Z = util.paulis[1:] From aab9dc6bed1c857dc5013c2781d173883c274fbc Mon Sep 17 00:00:00 2001 From: Tobias Hangleiter Date: Fri, 1 Nov 2024 16:18:49 +0100 Subject: [PATCH 03/26] Add progressbar --- filter_functions/numeric.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/filter_functions/numeric.py b/filter_functions/numeric.py index 3c2df84..e82114d 100644 --- a/filter_functions/numeric.py +++ b/filter_functions/numeric.py @@ -1781,7 +1781,8 @@ def calculate_second_order_from_atomic( phases: ndarray, propagators: ndarray, propagators_liouville: ndarray, - intermediates: Sequence[Mapping[str, ndarray]] + intermediates: Sequence[Mapping[str, ndarray]], + show_progressbar: bool = False, ): G, n_nops, n_basis, n_omega = control_matrix_atomic.shape d = propagators.shape[-1] @@ -1806,7 +1807,8 @@ def _incomplete_time_step(h, out): result = filter_function_atomic.copy() - for g in range(1, G): + for g in util.progressbar_range(1, G, disable=not show_progressbar, + desc='Calculating second order FF'): eigvecs_propagated = _propagate_eigenvectors(propagators[g-1:g], intermediates[g]['eigvecs_propagated']) n_opers_transformed = intermediates[g]['n_opers_transformed'] From 1c82b4149d8a0e4dbd2659bc250f4e8187396481 Mon Sep 17 00:00:00 2001 From: Tobias Hangleiter Date: Fri, 1 Nov 2024 18:20:25 +0100 Subject: [PATCH 04/26] Drop unused items at index 0 The "cumulative control matrix" for instance is zero for g=0. --- filter_functions/pulse_sequence.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/filter_functions/pulse_sequence.py b/filter_functions/pulse_sequence.py index 52d674c..e68468c 100644 --- a/filter_functions/pulse_sequence.py +++ b/filter_functions/pulse_sequence.py @@ -1844,7 +1844,7 @@ def concatenate( basis=newpulse.basis, filter_function_atomic=pulses[0].get_filter_function(omega, order=2), control_matrix_atomic=control_matrix_atomic, - control_matrix_accumulated=control_matrix_accumulated, + control_matrix_accumulated=control_matrix_accumulated[:-1], phases=phases, propagators=util.adot([pulse.total_propagator for pulse in pulses[:-1]]), propagators_liouville=propagators_liouville, From 8e43c25cf7408c47c9adc336b659e58db25c2bfe Mon Sep 17 00:00:00 2001 From: Tobias Hangleiter Date: Fri, 1 Nov 2024 18:35:45 +0100 Subject: [PATCH 05/26] Align notation --- filter_functions/numeric.py | 24 ++++++++++++------------ filter_functions/pulse_sequence.py | 6 +++--- tests/test_sequencing.py | 2 +- 3 files changed, 16 insertions(+), 16 deletions(-) diff --git a/filter_functions/numeric.py b/filter_functions/numeric.py index e82114d..8e44448 100644 --- a/filter_functions/numeric.py +++ b/filter_functions/numeric.py @@ -626,7 +626,7 @@ def calculate_control_matrix_from_atomic( propagators_liouville: ndarray, show_progressbar: bool = False, which: str = 'total', - return_accumulated: bool = False + return_cumulative: bool = False ) -> Union[ndarray, Tuple[ndarray, ndarray]]: r""" Calculate the control matrix from the control matrices of atomic @@ -649,7 +649,7 @@ def calculate_control_matrix_from_atomic( Compute the total control matrix (the sum of all time steps) or the correlation control matrix (first axis holds each pulses' contribution) - return_accumulated: bool, optional + return_cumulative: bool, optional Also return the accumulated sum, that is, an array that holds the control matrix of the sequence up to position *g* in element *g* of its first axis. If each atomic unit is a single segment, @@ -691,8 +691,8 @@ def calculate_control_matrix_from_atomic( else: # which == 'total' control_matrix = np.zeros(control_matrix_atomic.shape[1:], dtype=complex) - if return_accumulated: - control_matrix_accumulated = np.empty_like(control_matrix_atomic) + if return_cumulative: + control_matrix_atomic_cumulative = np.empty_like(control_matrix_atomic) else: # A buffer for intermediate terms in the calculation. tmp = np.empty_like(control_matrix) @@ -701,8 +701,8 @@ def calculate_control_matrix_from_atomic( desc='Calculating control matrix'): if which == 'correlations': tmp = control_matrix[g] - elif return_accumulated: - tmp = control_matrix_accumulated[g] + elif return_cumulative: + tmp = control_matrix_atomic_cumulative[g] # else: defined outside the loop if g > 0: @@ -714,11 +714,11 @@ def calculate_control_matrix_from_atomic( if which == 'total': control_matrix += tmp - if return_accumulated: - control_matrix_accumulated[g] = control_matrix + if return_cumulative: + control_matrix_atomic_cumulative[g] = control_matrix - if return_accumulated: - return control_matrix, control_matrix_accumulated + if return_cumulative: + return control_matrix, control_matrix_atomic_cumulative return control_matrix @@ -1777,7 +1777,7 @@ def calculate_second_order_from_atomic( basis: Basis, filter_function_atomic: ndarray, control_matrix_atomic: ndarray, - control_matrix_accumulated: ndarray, + control_matrix_atomic_cumulative: ndarray, phases: ndarray, propagators: ndarray, propagators_liouville: ndarray, @@ -1816,7 +1816,7 @@ def _incomplete_time_step(h, out): second_order_complete_steps = intermediates[g]['second_order_complete_steps'] # B'_(g-1)(ω) - tmp_1[:] = control_matrix_accumulated[g-1] + tmp_1[:] = control_matrix_atomic_cumulative[g - 1] tmp_1 *= phases[g-1].conj() tmp_1 = expr_1(tmp_1, propagators_liouville[g-1].T, out=tmp_1) diff --git a/filter_functions/pulse_sequence.py b/filter_functions/pulse_sequence.py index e68468c..71f4985 100644 --- a/filter_functions/pulse_sequence.py +++ b/filter_functions/pulse_sequence.py @@ -1835,16 +1835,16 @@ def concatenate( control_matrix = numeric.calculate_control_matrix_from_atomic( phases, control_matrix_atomic, propagators_liouville, show_progressbar, which='correlations' if calc_pulse_correlation_FF else 'total', - return_accumulated=calc_second_order_FF + return_cumulative=calc_second_order_FF ) if calc_second_order_FF: - control_matrix, control_matrix_accumulated = control_matrix + control_matrix, control_matrix_atomic_cumulative = control_matrix filter_function = numeric.calculate_second_order_from_atomic( basis=newpulse.basis, filter_function_atomic=pulses[0].get_filter_function(omega, order=2), control_matrix_atomic=control_matrix_atomic, - control_matrix_accumulated=control_matrix_accumulated[:-1], + control_matrix_atomic_cumulative=control_matrix_atomic_cumulative[:-1], phases=phases, propagators=util.adot([pulse.total_propagator for pulse in pulses[:-1]]), propagators_liouville=propagators_liouville, diff --git a/tests/test_sequencing.py b/tests/test_sequencing.py index b8e6ebd..0a82be7 100644 --- a/tests/test_sequencing.py +++ b/tests/test_sequencing.py @@ -131,7 +131,7 @@ def test_caching(self): np.array([pulse[:i].get_total_phases(omega) for i in range(1, len(pulse))]), np.array([p.get_control_matrix(omega) for p in pulse]), np.array([pulse[:i].total_propagator_liouville for i in range(1, len(pulse))]), - return_accumulated=True + return_cumulative=True ) self.assertArrayAlmostEqual(ctrlmat_cumulative[:-1], pulse.intermediates['control_matrix_step_cumulative'], From 041da8abc657cde0a1d6c60179c4c1d8a242128c Mon Sep 17 00:00:00 2001 From: Tobias Hangleiter Date: Sat, 2 Nov 2024 11:35:29 +0100 Subject: [PATCH 06/26] Cache cumulative second order --- filter_functions/numeric.py | 2 +- tests/test_sequencing.py | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/filter_functions/numeric.py b/filter_functions/numeric.py index 8e44448..8c0f38a 100644 --- a/filter_functions/numeric.py +++ b/filter_functions/numeric.py @@ -1816,7 +1816,7 @@ def _incomplete_time_step(h, out): second_order_complete_steps = intermediates[g]['second_order_complete_steps'] # B'_(g-1)(ω) - tmp_1[:] = control_matrix_atomic_cumulative[g - 1] + tmp_1[:] = control_matrix_atomic_cumulative[g-1] tmp_1 *= phases[g-1].conj() tmp_1 = expr_1(tmp_1, propagators_liouville[g-1].T, out=tmp_1) diff --git a/tests/test_sequencing.py b/tests/test_sequencing.py index 0a82be7..d29d93b 100644 --- a/tests/test_sequencing.py +++ b/tests/test_sequencing.py @@ -464,7 +464,6 @@ def test_second_order(self): pulse.get_filter_function(omega, order=2), atol=1e-13) - def test_different_n_opers(self): """Test behavior when concatenating with different n_opers.""" for d, n_dt in zip(rng.integers(2, 5, 20), From 022f4c2a4e71befcbc6fa7b0ac74bc53d7b5ac9f Mon Sep 17 00:00:00 2001 From: Tobias Hangleiter Date: Sat, 21 Dec 2024 09:30:50 +0100 Subject: [PATCH 07/26] Commentary and documentation --- filter_functions/numeric.py | 65 ++++++++++++++++++++++++++++-- filter_functions/pulse_sequence.py | 2 +- 2 files changed, 63 insertions(+), 4 deletions(-) diff --git a/filter_functions/numeric.py b/filter_functions/numeric.py index 8c0f38a..b7c1584 100644 --- a/filter_functions/numeric.py +++ b/filter_functions/numeric.py @@ -639,7 +639,7 @@ def calculate_control_matrix_from_atomic( :math:`g=0`, they are unity. control_matrix_atomic: array_like, shape (n_dt, n_nops, d**2, n_omega) The pulse control matrices for :math:`g\in\{1, 2, \dots, G\}`. - propagators_liouville: array_like, shape (n_dt-1, n_nops, d**2, d**2) + propagators_liouville: array_like, shape (n_dt-1, d**2, d**2) The transfer matrices of the cumulative propagators for :math:`g\in\{1, 2, \dots, G-1\}`. For :math:`g=0` it is the identity. @@ -1551,8 +1551,9 @@ def calculate_second_order_filter_function( ------- second_order_filter_function : ndarray, shape (n_nops, n_nops, d**2, d**2, n_omega) The second order filter function. - filter_function_cumulative_step : ndarray, shape (n_dt, n_nops, n_nops, d**2, d**2, n_omega) - The accumulated filter function. + intermediates : dict[str, ndarray] + Intermediate results of the calculation. Only if + cache_intermediates is True. .. _notes: @@ -1784,6 +1785,63 @@ def calculate_second_order_from_atomic( intermediates: Sequence[Mapping[str, ndarray]], show_progressbar: bool = False, ): + r"""Calculate the second-order filter function from those of atomic + segments. + + Parameters + ---------- + basis: Basis, shape (d**2, d, d) + The operator basis for the filter function. + filter_function_atomic: ndarray, shape (n_dt, n_nops, n_nops, d**2, d**2, n_omega) + The filter functions for segments :math:`g\in\{1, 2, \dots, G\}`. + control_matrix_atomic: ndarray, shape (n_dt, n_nops, d**2, n_omega) + The pulse control matrices for :math:`g\in\{1, 2, \dots, G\}`. + control_matrix_atomic_cumulative: ndarray, shape (n_dt, n_nops, d**2, n_omega) + The accumulated sum of atomic control matrices as returned by + :func:`calculate_control_matrix_from_atomic` with + ``return_cumulative=True``. + phases: ndarray, shape (n_dt-1, n_omega) + The phase factors for :math:`g\in\{1, 2, \dots, G-1\}`. For + :math:`g=0`, they are unity. + propagators: ndarray, shape (n_dt-1, d, d) + propagators_liouville: ndarray, shape (n_dt-1, d**2, d**2) + The transfer matrices of the cumulative propagators for + :math:`g\in\{1, 2, \dots, G-1\}`. For :math:`g=0` it is the + identity. + intermediates: Sequence[Dict[str, ndarray]} + Intermediate terms of the calculation of the control matrix and + second-order filter function that can be reused here. Each entry + of the sequence of length `n_pls` should be a dictionary with + the following required keys: + - 'eigvecs_propagated' + - 'n_opers_transformed' + - 'second_order_integral' + - 'second_order_complete_steps' + The first two are populated by + :func:`calculate_control_matrix_from_scratch`, the last two by + :func:`calculate_second_order_filter_function`. + show_progressbar: bool, optional + Show a progress bar for the calculation. + + Returns + ------- + second_order_FF: ndarray, shape (n_nops, n_nops, d**2, d**2, n_omega) + The second-order filter function + :math:`\mathcal{F}^{(2)}(\omega)`. + + See Also + -------- + calculate_control_matrix_from_scratch: Control matrix from scratch. + calculate_control_matrix_from_atomic: Similar function for the first order FF. + calculate_second_order_filter_function: Second-order FF from scratch. + """ + required_intermediates = {'eigvecs_propagated', 'n_opers_transformed', + 'second_order_integral', 'second_order_complete_steps'} + for required_key in required_intermediates: + if not all(required_key in intermediate for intermediate in intermediates): + raise ValueError(f"Required intermediate term {required_key} not found in all " + "intermediates.") + G, n_nops, n_basis, n_omega = control_matrix_atomic.shape d = propagators.shape[-1] @@ -1827,6 +1885,7 @@ def _incomplete_time_step(h, out): # N_(g)(ω) tmp_2 += second_order_complete_steps + # Q^(g-1) N_(g)(ω) Q^(g-1) result += expr_2(tmp_2, propagators_liouville[g-1], propagators_liouville[g-1], out=tmp_2) for h in range(len(eigvecs_propagated)): diff --git a/filter_functions/pulse_sequence.py b/filter_functions/pulse_sequence.py index 71f4985..9ed9700 100644 --- a/filter_functions/pulse_sequence.py +++ b/filter_functions/pulse_sequence.py @@ -1769,7 +1769,7 @@ def concatenate( if not equal_omega: if calc_filter_function: - raise ValueError("Calculation of filter function forced but not all pulses " + raise ValueError("Calculation of filter function forced but not all pulses " + "have the same frequencies cached and none were supplied!") if calc_pulse_correlation_FF: raise ValueError("Cannot compute the pulse correlation filter functions; do not " From 2d0c72f3504102f137142d25bd5e882479090aae Mon Sep 17 00:00:00 2001 From: Tobias Hangleiter Date: Sat, 21 Dec 2024 09:32:07 +0100 Subject: [PATCH 08/26] Rename to be consistent with control matrix functions --- filter_functions/numeric.py | 11 +++++++---- filter_functions/pulse_sequence.py | 4 ++-- tests/test_core.py | 2 +- 3 files changed, 10 insertions(+), 7 deletions(-) diff --git a/filter_functions/numeric.py b/filter_functions/numeric.py index b7c1584..246d6ac 100644 --- a/filter_functions/numeric.py +++ b/filter_functions/numeric.py @@ -84,7 +84,9 @@ 'calculate_control_matrix_periodic', 'calculate_cumulant_function', 'calculate_decay_amplitudes', 'calculate_filter_function', 'calculate_frequency_shifts', 'calculate_noise_operators_from_atomic', 'calculate_noise_operators_from_scratch', - 'calculate_pulse_correlation_filter_function', 'calculate_second_order_filter_function', + 'calculate_pulse_correlation_filter_function', + 'calculate_second_order_filter_function_from_scratch', + 'calculate_second_order_filter_function_from_atomic', 'diagonalize', 'error_transfer_matrix', 'infidelity'] @@ -1486,7 +1488,7 @@ def calculate_filter_function(control_matrix: ndarray, which: str = 'fidelity') return np.einsum(subscripts, control_matrix.conj(), control_matrix) -def calculate_second_order_filter_function( +def calculate_second_order_filter_function_from_scratch( eigvals: ndarray, eigvecs: ndarray, propagators: ndarray, @@ -1591,6 +1593,7 @@ def calculate_second_order_filter_function( infidelity: Compute the infidelity directly. pulse_sequence.concatenate: Concatenate ``PulseSequence`` objects. calculate_pulse_correlation_filter_function + calculate_second_order_filter_function_from_atomic: Second-order FF from atomic. """ G, d = eigvals.shape n_nops = len(n_coeffs) @@ -1774,7 +1777,7 @@ def calculate_pulse_correlation_filter_function(control_matrix: ndarray, return np.einsum(subscripts, control_matrix.conj(), control_matrix) -def calculate_second_order_from_atomic( +def calculate_second_order_filter_function_from_atomic( basis: Basis, filter_function_atomic: ndarray, control_matrix_atomic: ndarray, @@ -1833,7 +1836,7 @@ def calculate_second_order_from_atomic( -------- calculate_control_matrix_from_scratch: Control matrix from scratch. calculate_control_matrix_from_atomic: Similar function for the first order FF. - calculate_second_order_filter_function: Second-order FF from scratch. + calculate_second_order_filter_function_from_scratch: Second-order FF from scratch. """ required_intermediates = {'eigvecs_propagated', 'n_opers_transformed', 'second_order_integral', 'second_order_complete_steps'} diff --git a/filter_functions/pulse_sequence.py b/filter_functions/pulse_sequence.py index 9ed9700..8c5089c 100644 --- a/filter_functions/pulse_sequence.py +++ b/filter_functions/pulse_sequence.py @@ -868,7 +868,7 @@ def cache_filter_function( filter_function = numeric.calculate_filter_function(control_matrix, which) else: # order == 2 - filter_function = numeric.calculate_second_order_filter_function( + filter_function = numeric.calculate_second_order_filter_function_from_scratch( self.eigvals, self.eigvecs, self.propagators, self.omega, self.basis, self.n_opers, self.n_coeffs, self.dt, self._intermediates, show_progressbar, cache_intermediates, cache_second_order_cumulative @@ -1840,7 +1840,7 @@ def concatenate( if calc_second_order_FF: control_matrix, control_matrix_atomic_cumulative = control_matrix - filter_function = numeric.calculate_second_order_from_atomic( + filter_function = numeric.calculate_second_order_filter_function_from_atomic( basis=newpulse.basis, filter_function_atomic=pulses[0].get_filter_function(omega, order=2), control_matrix_atomic=control_matrix_atomic, diff --git a/tests/test_core.py b/tests/test_core.py index eccb0f7..31173ac 100644 --- a/tests/test_core.py +++ b/tests/test_core.py @@ -792,7 +792,7 @@ def test_second_order_filter_function(self): F_1 = pulse.get_filter_function(omega, order=2) # Test caching F_2 = pulse.get_filter_function(omega, order=2) - F_3 = numeric.calculate_second_order_filter_function( + F_3 = numeric.calculate_second_order_filter_function_from_scratch( pulse.eigvals, pulse.eigvecs, pulse.propagators, omega, pulse.basis, pulse.n_opers, pulse.n_coeffs, pulse.dt, show_progressbar=False, intermediates=None ) From c637a56c6f7bd40f8bfc43203f370901f1836c78 Mon Sep 17 00:00:00 2001 From: Tobias Hangleiter Date: Sat, 21 Dec 2024 09:34:29 +0100 Subject: [PATCH 09/26] Move function --- filter_functions/numeric.py | 130 ++++++++++++++++++------------------ 1 file changed, 65 insertions(+), 65 deletions(-) diff --git a/filter_functions/numeric.py b/filter_functions/numeric.py index 246d6ac..38d78ed 100644 --- a/filter_functions/numeric.py +++ b/filter_functions/numeric.py @@ -1712,71 +1712,6 @@ def _incomplete_time_step(g, out): return result -@util.parse_optional_parameters(which=('fidelity', 'generalized')) -def calculate_pulse_correlation_filter_function(control_matrix: ndarray, - which: str = 'fidelity') -> ndarray: - r"""Compute pulse correlation filter function from control matrix. - - Parameters - ---------- - control_matrix: array_like, shape (n_pulses, n_nops, d**2, n_omega) - The control matrix. - which : str, optional - Which filter function to return. Either 'fidelity' (default) or - 'generalized' (see :ref:`Notes `). - - Returns - ------- - filter_function_pc: ndarray, shape (n_pls, n_pls, n_nops, n_nops, [d**2, d**2,] n_omega) - The pulse correlation filter functions for each pulse and noise - operator correlations. The first two axes hold the pulse - correlations, the second two the noise correlations. - which : str, optional - Which filter function to return. Either 'fidelity' (default) or - 'generalized' (see :ref:`Notes `). - - .. _notes: - - Notes - ----- - The generalized pulse correlation filter function is given by - - .. math:: - - F_{\alpha\beta,kl}^{(gg')}(\omega) = \bigl[ - \mathcal{Q}^{(g'-1)\dagger} - \tilde{\mathcal{B}}^{(g')\dagger}(\omega) - \bigr]_{k\alpha} \bigl[ - \tilde{\mathcal{B}}^{(g)}(\omega)\mathcal{Q}^{(g-1)} - \bigr]_{\beta l} e^{i\omega(t_{g-1} - t_{g'-1})}, - - with :math:`\tilde{\mathcal{B}}^{(g)}` the control matrix of the - :math:`g`-th pulse. The fidelity pulse correlation function is - obtained by tracing out the basis indices, - - .. math:: - - F_{\alpha\beta}^{(gg')}(\omega) = - \sum_{k} F_{\alpha\beta,kk}^{(gg')}(\omega) - - See Also - -------- - calculate_control_matrix_from_scratch: Control matrix from scratch. - calculate_control_matrix_from_atomic: Control matrix from concatenation. - calculate_filter_function: Regular filter function. - """ - if control_matrix.ndim != 4: - raise ValueError('Expected control_matrix.ndim == 4.') - - if which == 'fidelity': - subscripts = 'gako,hbko->ghabo' - else: - # which == 'generalized' - subscripts = 'gako,hblo->ghabklo' - - return np.einsum(subscripts, control_matrix.conj(), control_matrix) - - def calculate_second_order_filter_function_from_atomic( basis: Basis, filter_function_atomic: ndarray, @@ -1898,6 +1833,71 @@ def _incomplete_time_step(h, out): return result +@util.parse_optional_parameters(which=('fidelity', 'generalized')) +def calculate_pulse_correlation_filter_function(control_matrix: ndarray, + which: str = 'fidelity') -> ndarray: + r"""Compute pulse correlation filter function from control matrix. + + Parameters + ---------- + control_matrix: array_like, shape (n_pulses, n_nops, d**2, n_omega) + The control matrix. + which : str, optional + Which filter function to return. Either 'fidelity' (default) or + 'generalized' (see :ref:`Notes `). + + Returns + ------- + filter_function_pc: ndarray, shape (n_pls, n_pls, n_nops, n_nops, [d**2, d**2,] n_omega) + The pulse correlation filter functions for each pulse and noise + operator correlations. The first two axes hold the pulse + correlations, the second two the noise correlations. + which : str, optional + Which filter function to return. Either 'fidelity' (default) or + 'generalized' (see :ref:`Notes `). + + .. _notes: + + Notes + ----- + The generalized pulse correlation filter function is given by + + .. math:: + + F_{\alpha\beta,kl}^{(gg')}(\omega) = \bigl[ + \mathcal{Q}^{(g'-1)\dagger} + \tilde{\mathcal{B}}^{(g')\dagger}(\omega) + \bigr]_{k\alpha} \bigl[ + \tilde{\mathcal{B}}^{(g)}(\omega)\mathcal{Q}^{(g-1)} + \bigr]_{\beta l} e^{i\omega(t_{g-1} - t_{g'-1})}, + + with :math:`\tilde{\mathcal{B}}^{(g)}` the control matrix of the + :math:`g`-th pulse. The fidelity pulse correlation function is + obtained by tracing out the basis indices, + + .. math:: + + F_{\alpha\beta}^{(gg')}(\omega) = + \sum_{k} F_{\alpha\beta,kk}^{(gg')}(\omega) + + See Also + -------- + calculate_control_matrix_from_scratch: Control matrix from scratch. + calculate_control_matrix_from_atomic: Control matrix from concatenation. + calculate_filter_function: Regular filter function. + """ + if control_matrix.ndim != 4: + raise ValueError('Expected control_matrix.ndim == 4.') + + if which == 'fidelity': + subscripts = 'gako,hbko->ghabo' + else: + # which == 'generalized' + subscripts = 'gako,hblo->ghabklo' + + return np.einsum(subscripts, control_matrix.conj(), control_matrix) + + def diagonalize(hamiltonian: ndarray, dt: Coefficients) -> Tuple[ndarray, ndarray, ndarray]: r"""Diagonalize a Hamiltonian. From 71b526c54e5616c26b9aa499b733a9b51d7034bf Mon Sep 17 00:00:00 2001 From: Tobias Hangleiter Date: Thu, 28 Aug 2025 09:19:46 +0200 Subject: [PATCH 10/26] Typo --- filter_functions/numeric.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/filter_functions/numeric.py b/filter_functions/numeric.py index 38d78ed..a473567 100644 --- a/filter_functions/numeric.py +++ b/filter_functions/numeric.py @@ -1757,7 +1757,7 @@ def calculate_second_order_filter_function_from_atomic( - 'second_order_complete_steps' The first two are populated by :func:`calculate_control_matrix_from_scratch`, the last two by - :func:`calculate_second_order_filter_function`. + :func:`calculate_second_order_filter_function_from_scratch`. show_progressbar: bool, optional Show a progress bar for the calculation. From c1a338767b663c7a740ed6ac0f4934c7f7cc4b29 Mon Sep 17 00:00:00 2001 From: Tobias Hangleiter Date: Thu, 28 Aug 2025 09:19:55 +0200 Subject: [PATCH 11/26] Increase test tolerance --- tests/test_sequencing.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/test_sequencing.py b/tests/test_sequencing.py index d29d93b..11c1c4d 100644 --- a/tests/test_sequencing.py +++ b/tests/test_sequencing.py @@ -135,7 +135,7 @@ def test_caching(self): ) self.assertArrayAlmostEqual(ctrlmat_cumulative[:-1], pulse.intermediates['control_matrix_step_cumulative'], - atol=1e-15) + atol=1e-14) def test_concatenate_without_filter_function(self): """Concatenate two Spin Echos without filter functions.""" @@ -1402,3 +1402,4 @@ def test_accuracy(self): # Test the eigenvalues and -vectors by the characteristic equation self.assertCorrectDiagonalization(remapped_pulse, atol=1e-14) + From 2fcb70b4a071dbecba68238f84fd77e6d4b6045e Mon Sep 17 00:00:00 2001 From: Tobias Hangleiter Date: Thu, 4 Sep 2025 09:59:38 +0200 Subject: [PATCH 12/26] Simplify second-order concatenation --- filter_functions/numeric.py | 75 +++++++++++++++--------------- filter_functions/pulse_sequence.py | 5 +- 2 files changed, 40 insertions(+), 40 deletions(-) diff --git a/filter_functions/numeric.py b/filter_functions/numeric.py index a473567..fe7a2cf 100644 --- a/filter_functions/numeric.py +++ b/filter_functions/numeric.py @@ -1677,15 +1677,18 @@ def _incomplete_time_step(g, out): if g > 0: ctrlmat_step = ctrlmat_step_cache[g] ctrlmat_step_cumulative = ctrlmat_step_cumulative_cache[g-1] + + # Assign references; writing to int_buf writes to int_cache[g]! if cache_intermediates: int_buf = int_cache[g] if cache_cumulative: step_buf = filter_function_step_cumulative[g] + # else: defined outside the loop # We use step_buf as a buffer for the last interval (with nested time # dependence) and afterwards the intervals up to the last (where the # time dependence separates and we can use previous result for the - # control matrix). opt_einsum seems to be faster than numpy here. + # control matrix). if g > 0: # all (complete) intervals up to last # αko,βlo->αβklo @@ -1716,8 +1719,8 @@ def calculate_second_order_filter_function_from_atomic( basis: Basis, filter_function_atomic: ndarray, control_matrix_atomic: ndarray, + control_matrix_atomic_step: ndarray, control_matrix_atomic_cumulative: ndarray, - phases: ndarray, propagators: ndarray, propagators_liouville: ndarray, intermediates: Sequence[Mapping[str, ndarray]], @@ -1730,19 +1733,18 @@ def calculate_second_order_filter_function_from_atomic( ---------- basis: Basis, shape (d**2, d, d) The operator basis for the filter function. - filter_function_atomic: ndarray, shape (n_dt, n_nops, n_nops, d**2, d**2, n_omega) - The filter functions for segments :math:`g\in\{1, 2, \dots, G\}`. control_matrix_atomic: ndarray, shape (n_dt, n_nops, d**2, n_omega) + filter_function_atomic: ndarray, shape (n_nops, n_nops, d**2, d**2, n_omega) + The filter function for the first segment, :math:`g = 1`. The pulse control matrices for :math:`g\in\{1, 2, \dots, G\}`. - control_matrix_atomic_cumulative: ndarray, shape (n_dt, n_nops, d**2, n_omega) - The accumulated sum of atomic control matrices as returned by - :func:`calculate_control_matrix_from_atomic` with - ``return_cumulative=True``. - phases: ndarray, shape (n_dt-1, n_omega) - The phase factors for :math:`g\in\{1, 2, \dots, G-1\}`. For - :math:`g=0`, they are unity. - propagators: ndarray, shape (n_dt-1, d, d) - propagators_liouville: ndarray, shape (n_dt-1, d**2, d**2) + control_matrix_atomic_step: ndarray, shape (G-1, n_nops, d**2, n_omega) + The pulse correlation control matrix, i.e., the summands of the + first-order concatenated control matrix of the pulse sequence, + :math:`\mathcal{G}^{(g)}(\omega)`. + control_matrix_atomic_cumulative: ndarray, shape (G, n_nops, d**2, n_omega) + propagators: ndarray, shape (G-1, d, d) + The cumulative propagators of the *G* pulses. + propagators_liouville: ndarray, shape (G-1, d**2, d**2) The transfer matrices of the cumulative propagators for :math:`g\in\{1, 2, \dots, G-1\}`. For :math:`g=0` it is the identity. @@ -1783,26 +1785,27 @@ def calculate_second_order_filter_function_from_atomic( G, n_nops, n_basis, n_omega = control_matrix_atomic.shape d = propagators.shape[-1] - tmp_1 = np.empty((n_nops, n_basis, n_omega), dtype=complex) - tmp_2 = np.empty((n_nops, n_nops, n_basis, n_basis, n_omega), dtype=complex) + # A temporary array to throw intermediate results into + tmp = np.empty((n_nops, n_nops, n_basis, n_basis, n_omega), dtype=complex) basis_transformed = np.empty((n_basis, d, d), dtype=complex) n_opers_basis = np.empty((n_nops, n_basis, d, d), dtype=complex) - expr_1 = oe.contract_expression('ajo,jk->ako', (n_nops, n_basis, n_omega), (n_basis, n_basis)) - expr_2 = oe.contract_expression('abpqo,pk,ql->abklo', + expr_N = oe.contract_expression('pk,abpqo,ql->abklo', + (n_basis, n_basis), (n_nops, n_nops, n_basis, n_basis, n_omega), - (n_basis, n_basis), (n_basis, n_basis)) - expr_3 = oe.contract_expression('oijmn,akij,blmn->abklo', (n_omega, d, d, d, d), - (n_nops, n_basis, d, d), (n_nops, n_basis, d, d)) + (n_basis, n_basis)) + expr_J = oe.contract_expression('akij,oijmn,blmn->abklo', + (n_nops, n_basis, d, d), + (n_omega, d, d, d, d), + (n_nops, n_basis, d, d)) def _incomplete_time_step(h, out): _transform_by_unitary(eigvecs_propagated[h], basis, out=basis_transformed) np.multiply(n_opers_transformed[:, h, None], basis_transformed.swapaxes(-1, -2), out=n_opers_basis) - return expr_3(second_order_integral[h], n_opers_basis, n_opers_basis, out=out) + return expr_J(n_opers_basis, second_order_integral[h], n_opers_basis, out=out) result = filter_function_atomic.copy() - for g in util.progressbar_range(1, G, disable=not show_progressbar, desc='Calculating second order FF'): eigvecs_propagated = _propagate_eigenvectors(propagators[g-1:g], @@ -1811,24 +1814,22 @@ def _incomplete_time_step(h, out): second_order_integral = intermediates[g]['second_order_integral'] second_order_complete_steps = intermediates[g]['second_order_complete_steps'] - # B'_(g-1)(ω) - tmp_1[:] = control_matrix_atomic_cumulative[g-1] - tmp_1 *= phases[g-1].conj() - tmp_1 = expr_1(tmp_1, propagators_liouville[g-1].T, out=tmp_1) - - # B^(g)*(ω) B'_(g-1)(ω) - tmp_2 = np.multiply(control_matrix_atomic[g, :, None, :, None].conj(), - tmp_1[None, :, None, :], out=tmp_2) - - # N_(g)(ω) - tmp_2 += second_order_complete_steps + # First, we compute N_(g|g-1..1)(ω), the term including all complete time steps: + # G*_(g)(ω) B_(g-1)(ω), just an outer product: abo,klo->abklo + result += np.multiply(control_matrix_atomic_step[g, :, None, :, None].conj(), + control_matrix_atomic_cumulative[g-1, None, :, None], + out=tmp) - # Q^(g-1) N_(g)(ω) Q^(g-1) - result += expr_2(tmp_2, propagators_liouville[g-1], propagators_liouville[g-1], out=tmp_2) + # Q^T_(g-1..1) N_(g)(ω) Q_(g-1..1) + result += expr_N(propagators_liouville[g-1], + second_order_complete_steps, + propagators_liouville[g-1], + out=tmp) + # Lastly, all incomplete timesteps for h in range(len(eigvecs_propagated)): - # J'_(g)(ω) - result += _incomplete_time_step(h, out=tmp_2) + # J_(g|g-1..1)(ω) + result += _incomplete_time_step(h, out=tmp) return result diff --git a/filter_functions/pulse_sequence.py b/filter_functions/pulse_sequence.py index 8c5089c..404d4ad 100644 --- a/filter_functions/pulse_sequence.py +++ b/filter_functions/pulse_sequence.py @@ -1844,13 +1844,12 @@ def concatenate( basis=newpulse.basis, filter_function_atomic=pulses[0].get_filter_function(omega, order=2), control_matrix_atomic=control_matrix_atomic, - control_matrix_atomic_cumulative=control_matrix_atomic_cumulative[:-1], - phases=phases, + control_matrix_atomic_step=control_matrix_atomic_step, + control_matrix_atomic_cumulative=control_matrix_atomic_cumulative, propagators=util.adot([pulse.total_propagator for pulse in pulses[:-1]]), propagators_liouville=propagators_liouville, intermediates=[pulse.intermediates for pulse in pulses] ) - newpulse.cache_filter_function(omega, filter_function=filter_function, order=2) # Set the attribute and calculate filter function (if the pulse correlation From 11f37f1b467960bfa85bf64540bd05b213a13124 Mon Sep 17 00:00:00 2001 From: Tobias Hangleiter Date: Thu, 4 Sep 2025 10:07:25 +0200 Subject: [PATCH 13/26] Improve performance, remove return_cumulative arg The latter is equivalent to setting which='correlations' and performing the accumulation by hand. --- filter_functions/numeric.py | 67 +++++++++++------------------- filter_functions/pulse_sequence.py | 16 ++++--- 2 files changed, 35 insertions(+), 48 deletions(-) diff --git a/filter_functions/numeric.py b/filter_functions/numeric.py index fe7a2cf..97db282 100644 --- a/filter_functions/numeric.py +++ b/filter_functions/numeric.py @@ -627,8 +627,7 @@ def calculate_control_matrix_from_atomic( control_matrix_atomic: ndarray, propagators_liouville: ndarray, show_progressbar: bool = False, - which: str = 'total', - return_cumulative: bool = False + which: str = 'total' ) -> Union[ndarray, Tuple[ndarray, ndarray]]: r""" Calculate the control matrix from the control matrices of atomic @@ -651,15 +650,6 @@ def calculate_control_matrix_from_atomic( Compute the total control matrix (the sum of all time steps) or the correlation control matrix (first axis holds each pulses' contribution) - return_cumulative: bool, optional - Also return the accumulated sum, that is, an array that holds - the control matrix of the sequence up to position *g* in element - *g* of its first axis. If each atomic unit is a single segment, - corresponds to the intermediate 'control_matrix_step_cumulative' - returned by :func:`calculate_control_matrix_from_scratch` if - *cache_intermediates* is True. - Only if *which* is 'total' (otherwise corresponds to - ``control_matrix.cumsum(axis=0)``). Returns ------- @@ -681,48 +671,39 @@ def calculate_control_matrix_from_atomic( liouville_representation: Liouville representation for a given basis. """ G = len(control_matrix_atomic) - # Set up a reusable contraction expression. In some cases it is faster to - # also contract the time dimension in the same expression instead of - # looping over it, but we don't distinguish here for readability. - expr = oe.contract_expression('ijo,jk->iko', - control_matrix_atomic.shape[1:], - propagators_liouville.shape[1:]) + if control_matrix_atomic.flags.c_contiguous: + def restore_memory_layout(x): return np.ascontiguousarray(x.swapaxes(-1, -2)) + elif control_matrix_atomic.flags.f_contigous: + def restore_memory_layout(x): return np.asfortranarray(x.swapaxes(-1, -2)) + else: + def restore_memory_layout(x): return x.swapaxes(-1, -2) + + # It is quite a bit faster to work with frequencies on the second-to-last axis + control_matrix_atomic = np.ascontiguousarray(control_matrix_atomic.swapaxes(-1, -2)) if which == 'correlations': control_matrix = np.empty_like(control_matrix_atomic) - else: - # which == 'total' - control_matrix = np.zeros(control_matrix_atomic.shape[1:], dtype=complex) - if return_cumulative: - control_matrix_atomic_cumulative = np.empty_like(control_matrix_atomic) - else: - # A buffer for intermediate terms in the calculation. - tmp = np.empty_like(control_matrix) - - for g in util.progressbar_range(G, show_progressbar=show_progressbar, + control_matrix[0] = control_matrix_atomic[0] + elif which == 'total': + # First time step is simply the first atomic control matrix + control_matrix = control_matrix_atomic[0].copy() + # A buffer for intermediate terms in the calculation. + step = np.empty_like(control_matrix) + + for g in util.progressbar_range(1, G, show_progressbar=show_progressbar, desc='Calculating control matrix'): if which == 'correlations': - tmp = control_matrix[g] - elif return_cumulative: - tmp = control_matrix_atomic_cumulative[g] + step = control_matrix[g] # else: defined outside the loop - if g > 0: - # For the first time step phases and propagators are 1 - tmp = np.multiply(phases[g-1], control_matrix_atomic[g], out=tmp) - tmp = expr(tmp, propagators_liouville[g-1], out=tmp) - else: - tmp[:] = control_matrix_atomic[g] + step = np.multiply(phases[g-1, :, None], control_matrix_atomic[g], out=step) + step @= propagators_liouville[g-1] if which == 'total': - control_matrix += tmp - if return_cumulative: - control_matrix_atomic_cumulative[g] = control_matrix - - if return_cumulative: - return control_matrix, control_matrix_atomic_cumulative + control_matrix += step - return control_matrix + # Make sure we give back copies just the way they were before + return restore_memory_layout(control_matrix) def calculate_control_matrix_from_scratch( diff --git a/filter_functions/pulse_sequence.py b/filter_functions/pulse_sequence.py index 404d4ad..4845740 100644 --- a/filter_functions/pulse_sequence.py +++ b/filter_functions/pulse_sequence.py @@ -1824,22 +1824,28 @@ def concatenate( pulse.dt, t=pulse.t, show_progressbar=show_progressbar, cache_intermediates=False ) - # Set the total propagator for possible future concatenations (if not done - # so above) + # Set the total propagator for possible future concatenations (if not done so above) if not newpulse.is_cached('total_propagator'): newpulse.total_propagator = util.mdot([pls.total_propagator for pls in pulses][::-1]) newpulse.cache_total_phases(omega) newpulse.total_propagator_liouville = liouville_representation(newpulse.total_propagator, newpulse.basis) + # 'correlations' corresponds to returning each summand of the sum over g, which is exactly what + # is also needed for the second-order filter function. control_matrix = numeric.calculate_control_matrix_from_atomic( phases, control_matrix_atomic, propagators_liouville, show_progressbar, - which='correlations' if calc_pulse_correlation_FF else 'total', - return_cumulative=calc_second_order_FF + which='correlations' if calc_pulse_correlation_FF or calc_second_order_FF else 'total', ) if calc_second_order_FF: - control_matrix, control_matrix_atomic_cumulative = control_matrix + # control_matrix is the step-wise one. + control_matrix_atomic_step = control_matrix + control_matrix_atomic_cumulative = control_matrix_atomic_step.cumsum(axis=0) + if which == 'total': + # restore the total control matrix for below + control_matrix = control_matrix_atomic_cumulative[-1] + filter_function = numeric.calculate_second_order_filter_function_from_atomic( basis=newpulse.basis, filter_function_atomic=pulses[0].get_filter_function(omega, order=2), From 2ea32125bd4859dd864c9c94f8f24aa13106679f Mon Sep 17 00:00:00 2001 From: Tobias Hangleiter Date: Thu, 4 Sep 2025 10:08:38 +0200 Subject: [PATCH 14/26] Simplify loop body --- filter_functions/numeric.py | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/filter_functions/numeric.py b/filter_functions/numeric.py index 97db282..09980dd 100644 --- a/filter_functions/numeric.py +++ b/filter_functions/numeric.py @@ -439,19 +439,16 @@ def calculate_noise_operators_from_atomic( calculate_control_matrix_from_atomic: Same calculation in Liouville space. """ G = len(noise_operators_atomic) - noise_operators = np.zeros(noise_operators_atomic.shape[1:], dtype=complex) + noise_operators = noise_operators_atomic[0].copy() tmp = np.empty_like(noise_operators) - for g in util.progressbar_range(G, show_progressbar=show_progressbar, + for g in util.progressbar_range(1, G, show_progressbar=show_progressbar, desc='Calculating noise operators'): - - if g > 0: - tmp = np.multiply(noise_operators_atomic[g], phases[g-1, :, None, None, None], out=tmp) - tmp = _transform_by_unitary(propagators[g-1], tmp, out=tmp) - else: - tmp = noise_operators_atomic[g] - - noise_operators += tmp + noise_operators += _transform_by_unitary( + propagators[g-1], + np.multiply(noise_operators_atomic[g], phases[g-1, :, None, None, None], out=tmp), + out=tmp + ) return noise_operators From abab99279250908f75074425c8d22d9e16430266 Mon Sep 17 00:00:00 2001 From: Tobias Hangleiter Date: Thu, 4 Sep 2025 10:08:44 +0200 Subject: [PATCH 15/26] Update tests --- tests/test_sequencing.py | 33 +++++++++++++++++++++------------ 1 file changed, 21 insertions(+), 12 deletions(-) diff --git a/tests/test_sequencing.py b/tests/test_sequencing.py index 11c1c4d..71f5f53 100644 --- a/tests/test_sequencing.py +++ b/tests/test_sequencing.py @@ -127,16 +127,6 @@ def test_caching(self): self.assertArrayEqual(ctrlmat, slc.get_control_matrix(omega)) self.assertArrayEqual(FF, slc.get_filter_function(omega, order=2)) - _, ctrlmat_cumulative = numeric.calculate_control_matrix_from_atomic( - np.array([pulse[:i].get_total_phases(omega) for i in range(1, len(pulse))]), - np.array([p.get_control_matrix(omega) for p in pulse]), - np.array([pulse[:i].total_propagator_liouville for i in range(1, len(pulse))]), - return_cumulative=True - ) - self.assertArrayAlmostEqual(ctrlmat_cumulative[:-1], - pulse.intermediates['control_matrix_step_cumulative'], - atol=1e-14) - def test_concatenate_without_filter_function(self): """Concatenate two Spin Echos without filter functions.""" tau = 10 @@ -449,9 +439,24 @@ def test_concatenate_split_cnot(self): rtol, atol) def test_second_order(self): - for d, n_dt in zip(rng.integers(2, 5, 10), rng.integers(1, 11, 10)): + for d, n_dt in zip(rng.integers(2, 5, 10), rng.integers(3, 11, 10)): pulse = testutil.rand_pulse_sequence(d, n_dt) omega = util.get_sample_frequencies(pulse, 11) + + # Split at random index + idx = rng.integers(1, n_dt-1) + pulses = pulse[:idx], pulse[idx:] + for pls in pulses: + pls.cleanup('frequency dependent') + pls.cache_filter_function(omega, order=1, cache_intermediates=True) + pls.cache_filter_function(omega, order=2, cache_intermediates=True) + + concat_pulse = ff.concatenate(pulses, calc_second_order_FF=True) + self.assertArrayAlmostEqual(concat_pulse.get_filter_function(omega, order=2), + pulse.get_filter_function(omega, order=2), + atol=1e-13) + + # Split into pulses with single time segments pulses = list(pulse) pulse.cache_filter_function(omega, order=1, cache_intermediates=True) pulse.cache_filter_function(omega, order=2, cache_intermediates=True) @@ -464,6 +469,11 @@ def test_second_order(self): pulse.get_filter_function(omega, order=2), atol=1e-13) + # Assert some assertions + pulses[rng.integers(0, len(pulses))]._intermediates.clear() + with self.assertRaises(ValueError): + ff.concatenate(pulses, calc_second_order_FF=True) + def test_different_n_opers(self): """Test behavior when concatenating with different n_opers.""" for d, n_dt in zip(rng.integers(2, 5, 20), @@ -1402,4 +1412,3 @@ def test_accuracy(self): # Test the eigenvalues and -vectors by the characteristic equation self.assertCorrectDiagonalization(remapped_pulse, atol=1e-14) - From 07e6e426927e0d3b76f4e33fb2c24aac95306dd7 Mon Sep 17 00:00:00 2001 From: Tobias Hangleiter Date: Thu, 4 Sep 2025 10:09:04 +0200 Subject: [PATCH 16/26] make doc more consistent --- filter_functions/numeric.py | 48 +++++++++++++++--------------- filter_functions/pulse_sequence.py | 13 ++++++-- 2 files changed, 35 insertions(+), 26 deletions(-) diff --git a/filter_functions/numeric.py b/filter_functions/numeric.py index 09980dd..ec702b6 100644 --- a/filter_functions/numeric.py +++ b/filter_functions/numeric.py @@ -385,13 +385,13 @@ def calculate_noise_operators_from_atomic( Parameters ---------- - phases: array_like, shape (n_dt-1, n_omega) + phases: array_like, shape (G-1, n_omega) The phase factors for :math:`g\in\{1, 2, \dots, G-1\}`. For :math:`g=0` they are unity. - noise_operators_atomic: array_like, shape (n_dt, n_nops, d, d, n_omega) + noise_operators_atomic: array_like, shape (G, n_nops, d, d, n_omega) The noise operators in the interaction picture of the g-th pulse, i.e. for :math:`g\in\{1, 2, \dots, G\}`. - propagators: array_like, shape (n_dt-1, d, d) + propagators: array_like, shape (G-1, d, d) The cumulative propagators of the pulses :math:`g\in\{1, 2, \dots, G-1\}`. For :math:`g=0` it is the identity. @@ -489,10 +489,10 @@ def calculate_noise_operators_from_scratch( n_coeffs: array_like, shape (n_nops, n_dt) The sensitivities of the system to the noise operators given by *n_opers* at the given time step. - dt: array_like, shape (n_dt) + dt: array_like, shape (n_dt,) Sequence duration, i.e. for the :math:`g`-th pulse :math:`t_g - t_{g-1}`. - t: array_like, shape (n_dt+1), optional + t: array_like, shape (n_dt+1,), optional The absolute times of the different segments. Can also be computed from *dt*. show_progressbar: bool, optional @@ -632,12 +632,12 @@ def calculate_control_matrix_from_atomic( Parameters ---------- - phases: array_like, shape (n_dt-1, n_omega) + phases: array_like, shape (G-1, n_omega) The phase factors for :math:`g\in\{1, 2, \dots, G-1\}`. For :math:`g=0`, they are unity. - control_matrix_atomic: array_like, shape (n_dt, n_nops, d**2, n_omega) + control_matrix_atomic: array_like, shape (G, n_nops, d**2, n_omega) The pulse control matrices for :math:`g\in\{1, 2, \dots, G\}`. - propagators_liouville: array_like, shape (n_dt-1, d**2, d**2) + propagators_liouville: array_like, shape (G-1, d**2, d**2) The transfer matrices of the cumulative propagators for :math:`g\in\{1, 2, \dots, G-1\}`. For :math:`g=0` it is the identity. @@ -650,7 +650,7 @@ def calculate_control_matrix_from_atomic( Returns ------- - control_matrix: ndarray, shape ([n_pls,] n_nops, d**2, n_omega) + control_matrix: ndarray, shape ([G,] n_nops, d**2, n_omega) The control matrix :math:`\tilde{\mathcal{B}}(\omega)`. Notes @@ -745,10 +745,10 @@ def calculate_control_matrix_from_scratch( n_coeffs: array_like, shape (n_nops, n_dt) The sensitivities of the system to the noise operators given by *n_opers* at the given time step. - dt: array_like, shape (n_dt) + dt: array_like, shape (n_dt,) Sequence duration, i.e. for the :math:`g`-th pulse :math:`t_g - t_{g-1}`. - t: array_like, shape (n_dt+1), optional + t: array_like, shape (n_dt+1,), optional The absolute times of the different segments. Can also be computed from *dt*. show_progressbar: bool, optional @@ -1003,10 +1003,10 @@ def calculate_cumulant_function( Also take into account the frequency shifts :math:`\Delta` that correspond to second order Magnus expansion and constitute unitary terms. Default ``False``. - decay_amplitudes, array_like, shape ([[n_pls, n_pls,] n_nops,] n_nops, d**2, d**2), optional + decay_amplitudes, array_like, shape ([[G, G,] n_nops,] n_nops, d**2, d**2), optional A precomputed cumulant function. If given, *spectrum*, *omega* are not required. - frequency_shifts, array_like, shape ([[n_pls, n_pls,] n_nops,] n_nops, d**2, d**2), optional + frequency_shifts, array_like, shape ([[G, G,] n_nops,] n_nops, d**2, d**2), optional A precomputed frequency shift. If given, *spectrum*, *omega* are not required for second order terms. show_progressbar: bool, optional @@ -1023,7 +1023,7 @@ def calculate_cumulant_function( Returns ------- - cumulant_function: ndarray, shape ([[n_pls, n_pls,] n_nops,] n_nops, d**2, d**2) + cumulant_function: ndarray, shape ([[G, G,] n_nops,] n_nops, d**2, d**2) The cumulant function. The individual noise operator contributions chosen by ``n_oper_identifiers`` are on the third to last axis / axes, depending on whether the noise is @@ -1251,7 +1251,7 @@ def calculate_decay_amplitudes( Returns ------- - decay_amplitudes: ndarray, shape ([[n_pls, n_pls,] n_nops,] n_nops, d**2, d**2) + decay_amplitudes: ndarray, shape ([[G, G,] n_nops,] n_nops, d**2, d**2) The decay amplitudes. .. _notes: @@ -1506,7 +1506,7 @@ def calculate_second_order_filter_function_from_scratch( n_coeffs: array_like, shape (n_nops, n_dt) The sensitivities of the system to the noise operators given by *n_opers* at the given time step. - dt: array_like, shape (n_dt) + dt: array_like, shape (n_dt,) Sequence duration, i.e. for the :math:`l`-th pulse :math:`t_l - t_{l-1}`. intermediates: Dict[str, ndarray], optional @@ -1711,9 +1711,9 @@ def calculate_second_order_filter_function_from_atomic( ---------- basis: Basis, shape (d**2, d, d) The operator basis for the filter function. - control_matrix_atomic: ndarray, shape (n_dt, n_nops, d**2, n_omega) filter_function_atomic: ndarray, shape (n_nops, n_nops, d**2, d**2, n_omega) The filter function for the first segment, :math:`g = 1`. + control_matrix_atomic: ndarray, shape (G, n_nops, d**2, n_omega) The pulse control matrices for :math:`g\in\{1, 2, \dots, G\}`. control_matrix_atomic_step: ndarray, shape (G-1, n_nops, d**2, n_omega) The pulse correlation control matrix, i.e., the summands of the @@ -1729,7 +1729,7 @@ def calculate_second_order_filter_function_from_atomic( intermediates: Sequence[Dict[str, ndarray]} Intermediate terms of the calculation of the control matrix and second-order filter function that can be reused here. Each entry - of the sequence of length `n_pls` should be a dictionary with + of the sequence of length `G` should be a dictionary with the following required keys: - 'eigvecs_propagated' - 'n_opers_transformed' @@ -1827,7 +1827,7 @@ def calculate_pulse_correlation_filter_function(control_matrix: ndarray, Returns ------- - filter_function_pc: ndarray, shape (n_pls, n_pls, n_nops, n_nops, [d**2, d**2,] n_omega) + filter_function_pc: ndarray, shape (G, G, n_nops, n_nops, [d**2, d**2,] n_omega) The pulse correlation filter functions for each pulse and noise operator correlations. The first two axes hold the pulse correlations, the second two the noise correlations. @@ -1895,9 +1895,9 @@ def diagonalize(hamiltonian: ndarray, dt: Coefficients) -> Tuple[ndarray, ndarra ---------- hamiltonian: array_like, shape (n_dt, d, d) Hamiltonian of shape (n_dt, d, d) with d the dimensionality of - the system - dt: array_like - The time differences + the system. + dt: array_like, shape (n_dt,) + The time differences. Returns ------- @@ -1970,7 +1970,7 @@ def error_transfer_matrix( Also take into account the frequency shifts :math:`\Delta` that correspond to second order Magnus expansion and constitute unitary terms. Default ``False``. - cumulant_function: ndarray, shape ([[n_pls, n_pls,] n_nops,] n_nops, d**2, d**2) + cumulant_function: ndarray, shape ([[G, G,] n_nops,] n_nops, d**2, d**2) A precomputed cumulant function. If given, *pulse*, *spectrum*, *omega* are not required. show_progressbar: bool, optional @@ -2137,7 +2137,7 @@ def infidelity( Returns ------- - infid: ndarray, shape ([[n_pls, n_pls,], n_nops,] n_nops) + infid: ndarray, shape ([[G, G,], n_nops,] n_nops) Array with the infidelity contributions for each spectrum *spectrum* on the last axis or axes, depending on the shape of *spectrum* and *which*. If ``which`` is ``correlations``, the diff --git a/filter_functions/pulse_sequence.py b/filter_functions/pulse_sequence.py index 4845740..1d4a408 100644 --- a/filter_functions/pulse_sequence.py +++ b/filter_functions/pulse_sequence.py @@ -647,7 +647,7 @@ def cache_control_matrix(self, omega: Coefficients, ---------- omega: array_like, shape (n_omega,) The frequencies for which to cache the filter function. - control_matrix: array_like, shape ([n_nops,] n_nops, d**2, n_omega), optional + control_matrix: array_like, shape ([G,] n_nops, d**2, n_omega), optional The control matrix for the frequencies *omega*. If ``None``, it is computed. show_progressbar: bool @@ -814,7 +814,7 @@ def cache_filter_function( ---------- omega: array_like, shape (n_omega,) The frequencies for which to cache the filter function. - control_matrix: array_like, shape ([n_nops,] n_nops, d**2, n_omega), optional + control_matrix: array_like, shape ([G,] n_nops, d**2, n_omega), optional The control matrix for the frequencies *omega*. If ``None``, it is computed and the filter function derived from it. filter_function: array_like, shape (n_nops, n_nops, [d**2, d**2,] n_omega), optional @@ -1695,6 +1695,15 @@ def concatenate( list of frequencies must be supplied as *omega*. Overridden if either *calc_pulse_correlation_FF* or *calc_second_order_FF* are true. + calc_second_order_FF : bool, optional + Compute the second-order filter function. Requires atomic pulses + to have retained `intermediates` during the calculation of the + control matrix. + + .. warning:: + This is an experimental feature and might have unexpected + bugs. + which: str, optional Which filter function to compute. Either 'fidelity' (default) or 'generalized' (see :meth:`PulseSequence.get_filter_function` and From 4b175651136d7ef8e306f3bd9751904fd4e3fcf3 Mon Sep 17 00:00:00 2001 From: Tobias Hangleiter Date: Thu, 4 Sep 2025 10:49:23 +0200 Subject: [PATCH 17/26] Fix conditionals --- filter_functions/numeric.py | 5 +++-- filter_functions/pulse_sequence.py | 2 +- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/filter_functions/numeric.py b/filter_functions/numeric.py index ec702b6..af6fa39 100644 --- a/filter_functions/numeric.py +++ b/filter_functions/numeric.py @@ -678,11 +678,12 @@ def restore_memory_layout(x): return x.swapaxes(-1, -2) # It is quite a bit faster to work with frequencies on the second-to-last axis control_matrix_atomic = np.ascontiguousarray(control_matrix_atomic.swapaxes(-1, -2)) + # First time step is simply the first atomic control matrix if which == 'correlations': control_matrix = np.empty_like(control_matrix_atomic) control_matrix[0] = control_matrix_atomic[0] - elif which == 'total': - # First time step is simply the first atomic control matrix + else: + # which == 'total' control_matrix = control_matrix_atomic[0].copy() # A buffer for intermediate terms in the calculation. step = np.empty_like(control_matrix) diff --git a/filter_functions/pulse_sequence.py b/filter_functions/pulse_sequence.py index 1d4a408..46056b8 100644 --- a/filter_functions/pulse_sequence.py +++ b/filter_functions/pulse_sequence.py @@ -1851,7 +1851,7 @@ def concatenate( # control_matrix is the step-wise one. control_matrix_atomic_step = control_matrix control_matrix_atomic_cumulative = control_matrix_atomic_step.cumsum(axis=0) - if which == 'total': + if not calc_pulse_correlation_FF: # restore the total control matrix for below control_matrix = control_matrix_atomic_cumulative[-1] From 46ec5ea65507ef756b970189184dd9dfbcd763a1 Mon Sep 17 00:00:00 2001 From: Tobias Hangleiter Date: Thu, 4 Sep 2025 11:03:15 +0200 Subject: [PATCH 18/26] Fix typo --- filter_functions/numeric.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/filter_functions/numeric.py b/filter_functions/numeric.py index af6fa39..5c31a92 100644 --- a/filter_functions/numeric.py +++ b/filter_functions/numeric.py @@ -670,7 +670,7 @@ def calculate_control_matrix_from_atomic( G = len(control_matrix_atomic) if control_matrix_atomic.flags.c_contiguous: def restore_memory_layout(x): return np.ascontiguousarray(x.swapaxes(-1, -2)) - elif control_matrix_atomic.flags.f_contigous: + elif control_matrix_atomic.flags.f_contiguous: def restore_memory_layout(x): return np.asfortranarray(x.swapaxes(-1, -2)) else: def restore_memory_layout(x): return x.swapaxes(-1, -2) From 35564e5238a6f60c16d42b98f4c7485f5bdb2260 Mon Sep 17 00:00:00 2001 From: Tobias Hangleiter Date: Thu, 4 Sep 2025 11:03:24 +0200 Subject: [PATCH 19/26] Add test for memory layout --- tests/test_sequencing.py | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/tests/test_sequencing.py b/tests/test_sequencing.py index 71f5f53..c25d4e0 100644 --- a/tests/test_sequencing.py +++ b/tests/test_sequencing.py @@ -62,6 +62,29 @@ def test_concatenate_base(self): pulse_2.omega = [3, 4] ff.concatenate([pulse_1, pulse_2], calc_filter_function=True) + # Make sure memory layout is not changed + omega = np.arange(1, 10) + phases = np.stack([pulse_1.get_total_phases(omega), + pulse_2.get_total_phases(omega)]) + control_matrix = np.stack([pulse_1.get_control_matrix(omega), + pulse_2.get_control_matrix(omega)]) + propagators = np.stack([pulse_1.total_propagator_liouville, + pulse_2.total_propagator_liouville]) + + concatenated = ff.numeric.calculate_control_matrix_from_atomic( + phases, + np.asfortranarray(control_matrix), + propagators + ) + self.assertTrue(concatenated.flags.f_contiguous) + + concatenated = ff.numeric.calculate_control_matrix_from_atomic( + phases, + np.ascontiguousarray(control_matrix.swapaxes(-1, -2)).swapaxes(-1, -2), + propagators + ) + self.assertFalse(concatenated.flags.contiguous) + def test_slicing(self): """Tests _getitem__.""" for d, n in zip(rng.integers(2, 5, 20), rng.integers(3, 51, 20)): From 993c72afa06f54be9d6776987d9aa21bb9a6fbe0 Mon Sep 17 00:00:00 2001 From: Tobias Hangleiter Date: Thu, 4 Sep 2025 14:20:56 +0200 Subject: [PATCH 20/26] Also test for c-major layout --- tests/test_sequencing.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/tests/test_sequencing.py b/tests/test_sequencing.py index c25d4e0..45ce192 100644 --- a/tests/test_sequencing.py +++ b/tests/test_sequencing.py @@ -71,6 +71,13 @@ def test_concatenate_base(self): propagators = np.stack([pulse_1.total_propagator_liouville, pulse_2.total_propagator_liouville]) + concatenated = ff.numeric.calculate_control_matrix_from_atomic( + phases, + np.ascontiguousarray(control_matrix), + propagators + ) + self.assertTrue(concatenated.flags.c_contiguous) + concatenated = ff.numeric.calculate_control_matrix_from_atomic( phases, np.asfortranarray(control_matrix), From 6a47c529fd6bf84823bfd212d36175ac0f78026d Mon Sep 17 00:00:00 2001 From: Tobias Hangleiter Date: Thu, 4 Sep 2025 15:24:49 +0200 Subject: [PATCH 21/26] Clean up second_order_from_scratch() This should hopefully be more readable and explicit --- filter_functions/numeric.py | 55 ++++++++++++++++++++----------------- 1 file changed, 30 insertions(+), 25 deletions(-) diff --git a/filter_functions/numeric.py b/filter_functions/numeric.py index 5c31a92..49d57bc 100644 --- a/filter_functions/numeric.py +++ b/filter_functions/numeric.py @@ -1588,8 +1588,8 @@ def calculate_second_order_filter_function_from_scratch( frc_bufs = (np.empty((n_omega, d, d), dtype=complex), np.empty((d, d, d, d), dtype=complex)) msk_bufs = np.empty((4, n_omega, d, d, d, d), dtype=bool) n_opers_basis_buf = np.empty((n_nops, n_basis, d, d), dtype=complex) - complete_step_buf = np.zeros((n_nops, n_nops, n_basis, n_basis, n_omega), dtype=complex) - result = np.zeros_like(complete_step_buf) + complete_steps = np.zeros((n_nops, n_nops, n_basis, n_basis, n_omega), dtype=complex) + incomplete_steps = np.zeros_like(complete_steps) # intermediate results from calculation of control matrix if intermediates is None: @@ -1600,7 +1600,6 @@ def calculate_second_order_filter_function_from_scratch( basis_transformed_cache = intermediates['basis_transformed'] ctrlmat_step_cache = intermediates['control_matrix_step'] ctrlmat_step_cumulative_cache = intermediates['control_matrix_step_cumulative'] - have_intermediates = True except KeyError: have_intermediates = False # No cache. Precompute some things and perform the costly computations @@ -1611,25 +1610,26 @@ def calculate_second_order_filter_function_from_scratch( basis_transformed = np.empty(basis.shape, dtype=complex) ctrlmat_step = np.zeros((n_nops, n_basis, n_omega), dtype=complex) ctrlmat_step_cumulative = np.zeros((n_nops, n_basis, n_omega), dtype=complex) + else: + have_intermediates = True if cache_intermediates: int_cache = np.empty((G, n_omega, d, d, d, d), dtype=complex) else: int_buf = np.empty((n_omega, d, d, d, d), dtype=complex) if cache_cumulative := cache_intermediates and cache_cumulative: - filter_function_step_cumulative = np.empty((G,) + (n_nops,)*2 + (n_basis,)*2 + (n_omega,), - dtype=complex) + cumulative_cache = np.empty((G,) + (n_nops,)*2 + (n_basis,)*2 + (n_omega,), dtype=complex) else: - step_buf = np.empty_like(complete_step_buf) + step_buf = np.empty_like(complete_steps) step_expr = oe.contract_expression('oijmn,akij,blmn->abklo', (n_omega, d, d, d, d), *[(n_nops, n_basis, d, d)]*2, optimize=[(0, 1), (0, 1)]) - def _incomplete_time_step(g, out): - _second_order_integral(omega, eigvals[g], dt[g], int_buf, frc_bufs, dE_bufs, msk_bufs) + def _incomplete_time_step(h, out): + _second_order_integral(omega, eigvals[h], dt[h], int_buf, frc_bufs, dE_bufs, msk_bufs) # αij,kji->αkij - np.multiply(n_opers_transformed[:, g, None], basis_transformed.swapaxes(-1, -2), + np.multiply(n_opers_transformed[:, h, None], basis_transformed.swapaxes(-1, -2), out=n_opers_basis_buf) return step_expr(int_buf, n_opers_basis_buf, n_opers_basis_buf, out=out) @@ -1661,34 +1661,39 @@ def _incomplete_time_step(g, out): if cache_intermediates: int_buf = int_cache[g] if cache_cumulative: - step_buf = filter_function_step_cumulative[g] + step_buf = cumulative_cache[g] # else: defined outside the loop # We use step_buf as a buffer for the last interval (with nested time # dependence) and afterwards the intervals up to the last (where the # time dependence separates and we can use previous result for the # control matrix). - if g > 0: - # all (complete) intervals up to last - # αko,βlo->αβklo - complete_step_buf += np.multiply(ctrlmat_step[:, None, :, None].conj(), - ctrlmat_step_cumulative[None, :, None], - out=step_buf) - # last (incomplete) interval - result += _incomplete_time_step(g, out=step_buf) - if cache_cumulative: - filter_function_step_cumulative[g] = result + complete_step_buf + # last (incomplete) interval + incomplete_steps += _incomplete_time_step(g, out=step_buf) + # if cache_cumulative, cumulative_cache[g] holds the gth incomplete step - if not cache_cumulative: - # if True, already done above - result += complete_step_buf + if g > 0: + # all (complete) intervals up to last. For g=0, this is zero. + # αko,βlo->αβklo + complete_steps += np.multiply(ctrlmat_step[:, None, :, None].conj(), + ctrlmat_step_cumulative[None, :, None], + out=step_buf) + if cache_cumulative: + # cumulative_cache[g] holds the gth complete step, so we add + # incomplete_steps, which holds everything else up to g. + cumulative_cache[g] = incomplete_steps + complete_steps + + if cache_cumulative: + result = cumulative_cache[-1] + else: + result = incomplete_steps + complete_steps if cache_intermediates: intermediates['second_order_integral'] = int_cache - intermediates['second_order_complete_steps'] = complete_step_buf + intermediates['second_order_complete_steps'] = complete_steps if cache_cumulative: - intermediates['filter_function_2_step_cumulative'] = filter_function_step_cumulative + intermediates['filter_function_2_step_cumulative'] = cumulative_cache return result, intermediates return result From dd5400b1545c1f1db72d8a476ce7e8bdabd164fc Mon Sep 17 00:00:00 2001 From: Tobias Hangleiter Date: Thu, 4 Sep 2025 15:24:56 +0200 Subject: [PATCH 22/26] Add doc --- filter_functions/pulse_sequence.py | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/filter_functions/pulse_sequence.py b/filter_functions/pulse_sequence.py index 46056b8..711f315 100644 --- a/filter_functions/pulse_sequence.py +++ b/filter_functions/pulse_sequence.py @@ -720,7 +720,20 @@ def get_filter_function( required by other computations. cache_second_order_cumulative: bool, optional Also cache the accumulated filter function for each time - step. Only if order is 2. + step. Only if order is 2. This can take up a lot of memory, + but is useful when interested in the time evolution. In that + case, compute the filter function for the entire pulse and + slice it afterwards:: + + pulse.cache_filter_function( + omega, order=2, cache_second_order_cumulative=True + ) + pulses_t = [] + for i in range(len(pulse)): + pulses_t.append(pulse[:i]) + + ``pulses_t`` will contain the filter functions for each time + step. Returns ------- From 7e661abbf088e56b14350bc29555d71d0ea8325d Mon Sep 17 00:00:00 2001 From: Tobias Hangleiter Date: Thu, 4 Sep 2025 16:17:57 +0200 Subject: [PATCH 23/26] Use opt_einsum, much faster --- filter_functions/superoperator.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/filter_functions/superoperator.py b/filter_functions/superoperator.py index 8dd3c47..02aabb9 100644 --- a/filter_functions/superoperator.py +++ b/filter_functions/superoperator.py @@ -40,6 +40,7 @@ from typing import Optional, Tuple, Union import numpy as np +import opt_einsum as oe from numpy import linalg as nla from numpy import ndarray from scipy import linalg as sla @@ -78,8 +79,8 @@ def liouville_representation(U: ndarray, basis: _b.Basis) -> ndarray: Hilbert space. """ U = np.asanyarray(U) - conjugated_basis = np.einsum('...ba,ibc,...cd->...iad', U.conj(), basis, U, - optimize=['einsum_path', (1, 2), (0, 1)]) + conjugated_basis = oe.contract('...ba,ibc,...cd->...iad', U.conj(), basis, U, + optimize=[(0, 1), (0, 1)]) return basis.expand(conjugated_basis, hermitian=basis.isherm) From f76a9f9ddccfdb034941efd10227d8aabf3a83ef Mon Sep 17 00:00:00 2001 From: Tobias Hangleiter Date: Thu, 4 Sep 2025 16:44:29 +0200 Subject: [PATCH 24/26] Doc --- .../examples/advanced_concatenation.ipynb | 76 ++++++++++++++++++- 1 file changed, 72 insertions(+), 4 deletions(-) diff --git a/doc/source/examples/advanced_concatenation.ipynb b/doc/source/examples/advanced_concatenation.ipynb index 76fbb25..13b5012 100644 --- a/doc/source/examples/advanced_concatenation.ipynb +++ b/doc/source/examples/advanced_concatenation.ipynb @@ -239,7 +239,9 @@ " for i in range(3):\n", " for j in range(3):\n", " t = ax[i, j].get_title()[7:]\n", - " ax[i, j].set_title('$F^{(' + pulses[i] + pulses[j] + t)" + " ax[i, j].set_title('$F^{(' + pulses[i] + pulses[j] + t)\n", + " fig.suptitle(f'{gate_type} gates')\n", + " fig.tight_layout()" ] }, { @@ -254,12 +256,78 @@ "\n", "While the imaginary part cancels out when calculating fidelities, $\\mathcal{I}\\propto\\sum_{ij} \\int\\mathrm{d}\\omega S(\\omega)F^{(ij)}(\\omega)$, the real part does not and the offdiagonals therefore lead to corrections in the total fidelity of a composite pulse, $\\mathcal{I}_\\text{tot}\\neq\\sum_g\\mathcal{I}^{(g)}$ with $\\mathcal{I}^{(g)}$ the infidelities of the individual pulses. These corrections can thus in principle also be negative, leading to improved fidelities for composite pulses." ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Second-order filter functions\n", + "We can also compute the second-order filter function as the concatenation of individual pulses. Although the calculation is not completely independent of the internal dynamics of the individual pulses as in the first-order case due to the nested time integral, a significant chunk can still be obtained from previously computed quantities. \n", + "\n", + "To do this, we need to make sure that some intermediate results are cached when computing the first- and second-order filter functions of the individual pulses, setting the `cache_intermediates` flag. We then just set the `calc_second_order_FF` argument to `True`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Clean up from before to have a clean slate for demonstration purposes\n", + "from itertools import chain\n", + "for key, pulse in chain(X2.items(), Y2.items()):\n", + " pulse.cleanup('frequency dependent')\n", + " pulse.cache_filter_function(omega[key], cache_intermediates=True,\n", + " order=1)\n", + " pulse.cache_filter_function(omega[key], cache_intermediates=True,\n", + " order=2)\n", + "\n", + "H = {key: ff.concatenate((Y2, X2, X2), calc_second_order_FF=True, which='generalized')\n", + " for (key, X2), (key, Y2) in zip(X2.items(), Y2.items())}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Visualize the FFs:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for key, pulse in H.items():\n", + " fig, axes = plt.subplots(3, 3, sharex=True, sharey=False,\n", + " figsize=(9, 6))\n", + " for i in range(1, 4):\n", + " for j in range(1, 4):\n", + " FF1 = pulse.get_filter_function(omega[key], order=1,\n", + " which='generalized')\n", + " FF2 = pulse.get_filter_function(omega[key], order=2,\n", + " which='generalized')\n", + " axes[i-1, j-1].semilogx(omega[key], FF1[0,0,i,j].real,\n", + " color='tab:blue', ls='-')\n", + " axes[i-1, j-1].semilogx(omega[key], FF1[0,0,i,j].imag,\n", + " color='tab:blue', ls='--')\n", + " axes[i-1, j-1].semilogx(omega[key], FF2[0,0,i,j].real,\n", + " color='tab:orange', ls='-')\n", + " axes[i-1, j-1].semilogx(omega[key], FF2[0,0,i,j].imag,\n", + " color='tab:orange', ls='--')\n", + " axes[i-1, j-1].margins(x=0)\n", + " fig.suptitle(f'{key} gates')\n", + " fig.supxlabel(r'$\\omega$')\n", + " fig.supylabel(r'$\\mathcal{F}(\\omega)$')\n", + " fig.tight_layout()" + ] } ], "metadata": { "kernelspec": { - "display_name": "Python 3 (Spyder)", - "language": "python3", + "display_name": "Python 3 (ipykernel)", + "language": "python", "name": "python3" }, "language_info": { @@ -272,7 +340,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.12.3" + "version": "3.13.5" }, "widgets": { "application/vnd.jupyter.widget-state+json": { From a02223b50b428f83c004b5873dd000b03c0e34a7 Mon Sep 17 00:00:00 2001 From: Tobias Hangleiter Date: Thu, 4 Sep 2025 17:00:40 +0200 Subject: [PATCH 25/26] Fix flaky test --- tests/test_sequencing.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/tests/test_sequencing.py b/tests/test_sequencing.py index 45ce192..5ff4092 100644 --- a/tests/test_sequencing.py +++ b/tests/test_sequencing.py @@ -530,7 +530,10 @@ def test_different_n_opers(self): n_coeffs[permutation], n_ids[permutation])), dt) - more_n_idx = sample(range(10), rng.integers(2, 5)) + + # draw more noise indices, but make sure they're not exactly the same + while sorted(more_n_idx := sample(range(10), rng.integers(2, 5))) == sorted(n_idx): + continue more_n_opers = opers[more_n_idx] more_n_coeffs = np.ones((more_n_opers.shape[0], n_dt)) more_n_coeffs *= np.abs(rng.standard_normal( From 2aebca1b948f23ab7874a228bebe96eb748abac6 Mon Sep 17 00:00:00 2001 From: Tobias Hangleiter Date: Thu, 4 Sep 2025 18:36:30 +0200 Subject: [PATCH 26/26] Import matplotlib.. --- doc/source/examples/advanced_concatenation.ipynb | 3 +++ 1 file changed, 3 insertions(+) diff --git a/doc/source/examples/advanced_concatenation.ipynb b/doc/source/examples/advanced_concatenation.ipynb index 13b5012..f42b3bb 100644 --- a/doc/source/examples/advanced_concatenation.ipynb +++ b/doc/source/examples/advanced_concatenation.ipynb @@ -275,6 +275,7 @@ "source": [ "# Clean up from before to have a clean slate for demonstration purposes\n", "from itertools import chain\n", + "\n", "for key, pulse in chain(X2.items(), Y2.items()):\n", " pulse.cleanup('frequency dependent')\n", " pulse.cache_filter_function(omega[key], cache_intermediates=True,\n", @@ -299,6 +300,8 @@ "metadata": {}, "outputs": [], "source": [ + "import matplotlib.pyplot as plt\n", + "\n", "for key, pulse in H.items():\n", " fig, axes = plt.subplots(3, 3, sharex=True, sharey=False,\n", " figsize=(9, 6))\n",