From 9bbac0fe4c50c1633dc8a21c2736155bb6ff626f Mon Sep 17 00:00:00 2001 From: Tobias Hangleiter Date: Tue, 2 Mar 2021 14:32:49 +0100 Subject: [PATCH 01/11] Rename and make private trafo functions more versatile --- filter_functions/numeric.py | 85 ++++++++++++++++++++++--------------- tests/test_core.py | 6 +-- 2 files changed, 54 insertions(+), 37 deletions(-) diff --git a/filter_functions/numeric.py b/filter_functions/numeric.py index 6114275..224dc3d 100644 --- a/filter_functions/numeric.py +++ b/filter_functions/numeric.py @@ -65,7 +65,7 @@ frequencies """ from collections import deque -from itertools import accumulate, repeat +from itertools import accumulate, repeat, zip_longest from typing import Any, Callable, Dict, Optional, Sequence, Tuple, Union from warnings import warn @@ -93,39 +93,49 @@ def _propagate_eigenvectors(propagators, eigvecs): return propagators.transpose(0, 2, 1).conj() @ eigvecs -def _transform_noise_operators(n_coeffs, n_opers, eigvecs): - r""" - Transform noise operators into the eigenspaces spanned by eigvecs. +def _transform_hamiltonian(eigvecs, opers, coeffs=None): + r"""Transform a Hamiltonian into the eigenspaces spanned by eigvecs. + I.e., the following transformation is performed: .. math:: - B_\alpha\rightarrow s_\alpha^{(g)}V^{(g)}B_\alpha V^{(g)\dagger} + s_\alpha^{(g)} B_\alpha\rightarrow + s_\alpha^{(g)} V^{(g)}B_\alpha V^{(g)\dagger} + + where :math:`s_\alpha^{(g)}` are coefficients of the operator + :math:`B_\alpha`. """ - assert len(n_opers) == len(n_coeffs) - n_opers_transformed = np.empty((len(n_opers), *eigvecs.shape), dtype=complex) - for j, (n_coeff, n_oper) in enumerate(zip(n_coeffs, n_opers)): - n_opers_transformed[j] = n_oper @ eigvecs - n_opers_transformed[j] = eigvecs.conj().transpose(0, 2, 1) @ n_opers_transformed[j] - n_opers_transformed[j] *= n_coeff[:, None, None] + if coeffs is None: + coeffs = [] + else: + assert len(opers) == len(coeffs) - return n_opers_transformed + opers_transformed = np.empty((len(opers), *eigvecs.shape), dtype=complex) + for j, (coeff, oper) in enumerate(zip_longest(coeffs, opers, fillvalue=None)): + opers_transformed[j] = _transform_by_unitary(eigvecs, oper, out=opers_transformed[j]) + if coeff is not None: + opers_transformed[j] *= coeff[:, None, None] + return opers_transformed -def _transform_basis(basis, eigvecs_propagated, out=None): - r""" - Transform the basis into the eigenspace spanned by V propagated by Q + +def _transform_by_unitary(unitary, oper, out=None): + r"""Transform the operators by a unitary. Uses broadcasting. I.e., the following transformation is performed: .. math:: - C_k\rightarrow Q_{g-1}V^{(g)\dagger}C_k V^{(g)}Q_{g-1}^\dagger. + C_k\rightarrow U C_k U^\dagger. """ - out = np.matmul(basis, eigvecs_propagated, out=out) - out = np.matmul(eigvecs_propagated.conj().T, out, out=out) + if out is None: + out = np.empty(oper.shape, dtype=oper.dtype) + + out = np.matmul(oper, unitary, out=out) + out = np.matmul(unitary.conj().swapaxes(-1, -2), out, out=out) return out @@ -715,7 +725,7 @@ def calculate_control_matrix_from_scratch( cache_intermediates: bool, optional Keep and return intermediate terms :math:`\mathcal{G}^{(g)}(\omega)` of the sum so that - :math:`\mathcal{R}(\omega)=\sum_g\mathcal{G}^{(g)}(\omega)`. + :math:`\mathcal{B}(\omega)=\sum_g\mathcal{G}^{(g)}(\omega)`. Otherwise the sum is performed in-place. out: ndarray, optional A location into which the result is stored. See @@ -771,28 +781,29 @@ def calculate_control_matrix_from_scratch( # Precompute noise opers transformed to eigenbasis of each pulse segment # and Q^\dagger @ V eigvecs_propagated = _propagate_eigenvectors(propagators[:-1], eigvecs) - n_opers_transformed = _transform_noise_operators(n_coeffs, n_opers, eigvecs) + n_opers_transformed = _transform_hamiltonian(eigvecs, n_opers, n_coeffs) # Allocate result and buffers for intermediate arrays - exp_buf, int_buf = np.empty((2, len(omega), d, d), dtype=complex) - + exp_buf = np.empty((len(omega), d, d), dtype=complex) if out is None: - control_matrix = np.zeros((len(n_opers), len(basis), len(omega)), dtype=complex) - else: - control_matrix = out + out = np.zeros((len(n_opers), len(basis), len(omega)), dtype=complex) if cache_intermediates: basis_transformed_cache = np.empty((len(dt), *basis.shape), dtype=complex) + phase_factors_cache = np.empty((len(dt), len(omega)), dtype=complex) + int_cache = np.empty((len(dt), len(omega), d, d), dtype=complex) sum_cache = np.empty((len(dt), len(n_opers), len(basis), len(omega)), dtype=complex) else: basis_transformed = np.empty(basis.shape, dtype=complex) + phase_factors = np.empty(len(omega), dtype=complex) + int_buf = np.empty((len(omega), d, d), dtype=complex) sum_buf = np.empty((len(n_opers), len(basis), len(omega)), dtype=complex) # Optimize the contraction path dynamically since it differs for different # values of d expr = oe.contract_expression('o,jmn,omn,knm->jko', omega.shape, n_opers_transformed[:, 0].shape, - int_buf.shape, basis.shape, optimize=True) + exp_buf.shape, basis.shape, optimize=True) for g in util.progressbar_range(len(dt), show_progressbar=show_progressbar, desc='Calculating control matrix'): @@ -800,22 +811,28 @@ def calculate_control_matrix_from_scratch( # Assign references to the locations in the cache for the quantities # that should be stored basis_transformed = basis_transformed_cache[g] + phase_factors = phase_factors_cache[g] + int_buf = int_cache[g] sum_buf = sum_cache[g] - basis_transformed = _transform_basis(basis, eigvecs_propagated[g], out=basis_transformed) + basis_transformed = _transform_by_unitary(eigvecs_propagated[g], basis, + out=basis_transformed) + phase_factors = util.cexp(omega*t[g], out=phase_factors) int_buf = _first_order_integral(omega, eigvals[g], dt[g], exp_buf, int_buf) - sum_buf = expr(util.cexp(omega*t[g]), n_opers_transformed[:, g], int_buf, + sum_buf = expr(phase_factors, n_opers_transformed[:, g], int_buf, basis_transformed, out=sum_buf) - control_matrix += sum_buf + out += sum_buf if cache_intermediates: intermediates = dict(n_opers_transformed=n_opers_transformed, basis_transformed=basis_transformed_cache, + phase_factors=phase_factors_cache, + first_order_integral=int_cache, control_matrix_step=sum_cache) - return control_matrix, intermediates + return out, intermediates - return control_matrix + return out def calculate_control_matrix_periodic(phases: ndarray, control_matrix: ndarray, @@ -1478,7 +1495,7 @@ def calculate_second_order_filter_function( t = np.concatenate(([0], np.asarray(dt).cumsum())) # Cheap to precompute as these don't use a lot of memory eigvecs_propagated = _propagate_eigenvectors(propagators[:-1], eigvecs) - n_opers_transformed = _transform_noise_operators(n_coeffs, n_opers, eigvecs) + n_opers_transformed = _transform_hamiltonian(eigvecs, n_opers, n_coeffs) # These are populated anew during every iteration, so there is no need # to keep every time step basis_transformed = np.empty(basis.shape, dtype=complex) @@ -1494,8 +1511,8 @@ def calculate_second_order_filter_function( for g in util.progressbar_range(len(dt), show_progressbar=show_progressbar, desc='Calculating second order FF'): if not intermediates: - basis_transformed = _transform_basis(basis, eigvecs_propagated[g], - out=basis_transformed) + basis_transformed = _transform_by_unitary(eigvecs_propagated[g], basis, + out=basis_transformed) # Need to compute G^(g) since no cache given. First initialize # buffer to zero. There is a probably lots of overhead computing # this individually for every time step. diff --git a/tests/test_core.py b/tests/test_core.py index 8874ddd..0e1cca4 100644 --- a/tests/test_core.py +++ b/tests/test_core.py @@ -556,9 +556,9 @@ def test_cache_intermediates(self): self.assertArrayAlmostEqual(pulse._intermediates['control_matrix_step'].sum(0), ctrlmat) self.assertArrayAlmostEqual(numeric.calculate_filter_function(ctrlmat), filtfun) self.assertArrayAlmostEqual(pulse._intermediates['n_opers_transformed'], - numeric._transform_noise_operators(pulse.n_coeffs, - pulse.n_opers, - pulse.eigvecs)) + numeric._transform_hamiltonian(pulse.eigvecs, + pulse.n_opers, + pulse.n_coeffs)) eigvecs_prop = numeric._propagate_eigenvectors(pulse.propagators[:-1], pulse.eigvecs) basis_transformed = np.einsum('gba,kbc,gcd->gkad', eigvecs_prop.conj(), pulse.basis, eigvecs_prop) From 96f9ae6bea94ac9978c34be7620552b7d2cff532 Mon Sep 17 00:00:00 2001 From: Tobias Hangleiter Date: Tue, 2 Mar 2021 14:35:13 +0100 Subject: [PATCH 02/11] get_indices_from_identifiers finds indices of subset in a larger set --- filter_functions/numeric.py | 6 +++--- filter_functions/plotting.py | 9 +++++---- filter_functions/pulse_sequence.py | 2 +- filter_functions/util.py | 29 +++++++++-------------------- tests/test_util.py | 14 +++++++------- 5 files changed, 25 insertions(+), 35 deletions(-) diff --git a/filter_functions/numeric.py b/filter_functions/numeric.py index 224dc3d..09d33b4 100644 --- a/filter_functions/numeric.py +++ b/filter_functions/numeric.py @@ -1197,7 +1197,7 @@ def calculate_decay_amplitudes( """ # TODO: Replace infidelity() by this? # Noise operator indices - idx = util.get_indices_from_identifiers(pulse, n_oper_identifiers, 'noise') + idx = util.get_indices_from_identifiers(pulse.n_oper_identifiers, n_oper_identifiers) if which == 'total': # Faster to use filter function instead of control matrix if pulse.is_cached('filter_function_gen'): @@ -1314,7 +1314,7 @@ def calculate_frequency_shifts( pulse_sequence.concatenate: Concatenate ``PulseSequence`` objects. calculate_pulse_correlation_filter_function """ - idx = util.get_indices_from_identifiers(pulse, n_oper_identifiers, 'noise') + idx = util.get_indices_from_identifiers(pulse.n_oper_identifiers, n_oper_identifiers) filter_function_2 = pulse.get_filter_function(omega, order=2, show_progressbar=show_progressbar) integrand = _get_integrand(spectrum, omega, idx, which_pulse='total', which_FF='generalized', @@ -1947,7 +1947,7 @@ def infidelity(pulse: 'PulseSequence', spectrum: Union[Coefficients, Callable], plotting.plot_infidelity_convergence: Convenience function to plot results. """ # Noise operator indices - idx = util.get_indices_from_identifiers(pulse, n_oper_identifiers, 'noise') + idx = util.get_indices_from_identifiers(pulse.n_oper_identifiers, n_oper_identifiers) if test_convergence: if not callable(spectrum): diff --git a/filter_functions/plotting.py b/filter_functions/plotting.py index 60a0b21..786bb98 100644 --- a/filter_functions/plotting.py +++ b/filter_functions/plotting.py @@ -315,7 +315,7 @@ def plot_pulse_train( ValueError If an invalid number of c_oper_labels were given """ - c_oper_inds = util.get_indices_from_identifiers(pulse, c_oper_identifiers, 'control') + c_oper_inds = util.get_indices_from_identifiers(pulse.c_oper_identifiers, c_oper_identifiers) c_oper_identifiers = pulse.c_oper_identifiers[c_oper_inds] if fig is None and axes is None: @@ -425,7 +425,7 @@ def plot_filter_function( else: omega = pulse.omega - n_oper_inds = util.get_indices_from_identifiers(pulse, n_oper_identifiers, 'noise') + n_oper_inds = util.get_indices_from_identifiers(pulse.n_oper_identifiers, n_oper_identifiers) n_oper_identifiers = pulse.n_oper_identifiers[n_oper_inds] if fig is None and axes is None: @@ -548,7 +548,7 @@ def plot_pulse_correlation_filter_function( If the pulse correlation filter function was not computed during concatenation. """ - n_oper_inds = util.get_indices_from_identifiers(pulse, n_oper_identifiers, 'noise') + n_oper_inds = util.get_indices_from_identifiers(pulse.n_oper_identifiers, n_oper_identifiers) n_oper_identifiers = pulse.n_oper_identifiers[n_oper_inds] diag_idx = np.arange(len(pulse.n_opers)) F_pc = pulse.get_pulse_correlation_filter_function() @@ -777,7 +777,8 @@ def plot_cumulant_function( raise ValueError('Require either precomputed cumulant function ' + 'or pulse, spectrum, and omega as arguments.') - n_oper_inds = util.get_indices_from_identifiers(pulse, n_oper_identifiers, 'noise') + n_oper_inds = util.get_indices_from_identifiers(pulse.n_oper_identifiers, + n_oper_identifiers) n_oper_identifiers = pulse.n_oper_identifiers[n_oper_inds] K = numeric.calculate_cumulant_function(pulse, spectrum, omega, n_oper_identifiers, 'total', second_order) diff --git a/filter_functions/pulse_sequence.py b/filter_functions/pulse_sequence.py index cb1f0dd..8c4a832 100644 --- a/filter_functions/pulse_sequence.py +++ b/filter_functions/pulse_sequence.py @@ -2439,7 +2439,7 @@ def extend( if additional_noise_Hamiltonian is not None: newpulse_n_oper_inds = util.get_indices_from_identifiers( - newpulse, n_oper_identifiers[n_ops_counter:], 'noise' + newpulse.n_oper_identifiers, n_oper_identifiers[n_ops_counter:] ) control_matrix[n_ops_counter:] = numeric.calculate_control_matrix_from_scratch( newpulse.eigvals, newpulse.eigvecs, newpulse.propagators, omega, newpulse.basis, diff --git a/filter_functions/util.py b/filter_functions/util.py index 8ec12a5..3f0bb26 100644 --- a/filter_functions/util.py +++ b/filter_functions/util.py @@ -26,8 +26,8 @@ :func:`abs2` Absolute value squared :func:`get_indices_from_identifiers` - The the indices of control or noise operators with given identifiers - as they are saved in a ``PulseSequence``. + The the indices of a subset of identifiers within a list of + identifiers. :func:`tensor` Fast, flexible tensor product of an arbitrary number of inputs using :func:`~numpy.einsum` @@ -221,32 +221,21 @@ def _parse_dims_arg(name: str, dims: Sequence[Sequence[int]], rank: int) -> None raise ValueError(f'Require all lists in {name}_dims to be of same length!') -@parse_optional_parameters({'kind': ('noise', 'control')}) -def get_indices_from_identifiers(pulse: 'PulseSequence', - identifiers: Union[None, str, Sequence[str]], - kind: str) -> Tuple[Sequence[int], - Sequence[str]]: +def get_indices_from_identifiers(all_identifiers: Sequence[str], + identifiers: Union[None, str, Sequence[str]]) -> Sequence[int]: """Get the indices of operators for given identifiers. Parameters ---------- - pulse: PulseSequence - The PulseSequence instance for which to get the indices. + all_identifiers: sequence of str + All available identifiers. identifiers: str or sequence of str The identifiers whose indices to get. - kind: str - Whether to get 'control' or 'noise' operator indices. """ - if kind == 'noise': - pulse_identifiers = pulse.n_oper_identifiers - else: - # kind == 'control' - pulse_identifiers = pulse.c_oper_identifiers - identifier_to_index_table = {identifier: index for index, identifier - in enumerate(pulse_identifiers)} + in enumerate(all_identifiers)} if identifiers is None: - inds = np.arange(len(pulse_identifiers)) + inds = np.arange(len(all_identifiers)) else: try: if isinstance(identifiers, str): @@ -256,7 +245,7 @@ def get_indices_from_identifiers(pulse: 'PulseSequence', for identifier in identifiers]) except KeyError: raise ValueError('Invalid identifiers given. All available ones ' + - f'are: {pulse_identifiers}') + f'are: {all_identifiers}') return inds diff --git a/tests/test_util.py b/tests/test_util.py index 7294754..c23cb1b 100644 --- a/tests/test_util.py +++ b/tests/test_util.py @@ -55,26 +55,26 @@ def test_get_indices_from_identifiers(self): [[util.paulis[2], [2]]], [1] ) - idx = util.get_indices_from_identifiers(pulse, ['X'], 'control') + idx = util.get_indices_from_identifiers(pulse.c_oper_identifiers, ['X']) self.assertArrayEqual(idx, [0]) - idx = util.get_indices_from_identifiers(pulse, 'X', 'control') + idx = util.get_indices_from_identifiers(pulse.c_oper_identifiers, 'X') self.assertArrayEqual(idx, [0]) - idx = util.get_indices_from_identifiers(pulse, ['Z', 'X'], 'control') + idx = util.get_indices_from_identifiers(pulse.c_oper_identifiers, ['Z', 'X']) self.assertArrayEqual(idx, [1, 0]) - idx = util.get_indices_from_identifiers(pulse, None, 'control') + idx = util.get_indices_from_identifiers(pulse.c_oper_identifiers, None) self.assertArrayEqual(idx, [0, 1]) - idx = util.get_indices_from_identifiers(pulse, ['B_0'], 'noise') + idx = util.get_indices_from_identifiers(pulse.n_oper_identifiers, ['B_0']) self.assertArrayEqual(idx, [0]) - idx = util.get_indices_from_identifiers(pulse, 'B_0', 'noise') + idx = util.get_indices_from_identifiers(pulse.n_oper_identifiers, 'B_0') self.assertArrayEqual(idx, [0]) with self.assertRaises(ValueError): - util.get_indices_from_identifiers(pulse, ['foobar'], 'noise') + util.get_indices_from_identifiers(pulse.n_oper_identifiers, ['foobar']) def test_tensor(self): shapes = [(1, 2, 3, 4, 5), (5, 4, 3, 2, 1)] From a74d950936722c1bebfdaebaa690d46efe5761b7 Mon Sep 17 00:00:00 2001 From: Tobias Hangleiter Date: Tue, 2 Mar 2021 14:39:50 +0100 Subject: [PATCH 03/11] Add cache_intermediates functionality to calc_noise_opers... --- filter_functions/numeric.py | 35 +++++++++++++++++++++++++---------- 1 file changed, 25 insertions(+), 10 deletions(-) diff --git a/filter_functions/numeric.py b/filter_functions/numeric.py index 09d33b4..29609c5 100644 --- a/filter_functions/numeric.py +++ b/filter_functions/numeric.py @@ -157,7 +157,7 @@ def _first_order_integral(E: ndarray, eigvals: ndarray, dt: float, int_buf.imag = np.add.outer(E, dE, out=int_buf.imag) # Catch zero-division - mask = (int_buf.imag != 0) + mask = (np.abs(int_buf.imag) > 1e-7) exp_buf = util.cexp(int_buf.imag*dt, out=exp_buf, where=mask) exp_buf = np.subtract(exp_buf, 1, out=exp_buf, where=mask) int_buf = np.divide(exp_buf, int_buf, out=int_buf, where=mask) @@ -473,7 +473,8 @@ def calculate_noise_operators_from_scratch( n_coeffs: Sequence[Coefficients], dt: Coefficients, t: Optional[Coefficients] = None, - show_progressbar: bool = False + show_progressbar: bool = False, + cache_intermediates: bool = False ) -> ndarray: r""" Calculate the noise operators in interaction picture from scratch. @@ -573,30 +574,44 @@ def calculate_noise_operators_from_scratch( n_coeffs = np.asarray(n_coeffs) # Precompute noise opers transformed to eigenbasis of each pulse - # segment and Q^\dagger @ V + # segment and V^\dagger @ Q eigvecs_propagated = _propagate_eigenvectors(eigvecs, propagators[:-1]) - n_opers_transformed = _transform_noise_operators(n_coeffs, n_opers, eigvecs) + n_opers_transformed = _transform_hamiltonian(eigvecs, n_opers, n_coeffs) # Allocate memory exp_buf, int_buf = np.empty((2, len(omega), d, d), dtype=complex) - intermediate = np.empty((len(omega), len(n_opers), d, d), dtype=complex) noise_operators = np.zeros((len(omega), len(n_opers), d, d), dtype=complex) + if cache_intermediates: + sum_cache = np.empty((len(dt), len(omega), len(n_opers), d, d), dtype=complex) + else: + sum_buf = np.empty((len(omega), len(n_opers), d, d), dtype=complex) + # Set up reusable expressions expr_1 = oe.contract_expression('akl,okl->oakl', n_opers_transformed[:, 0].shape, int_buf.shape) expr_2 = oe.contract_expression('ji,...jk,kl', - eigvecs_propagated[0].shape, intermediate.shape, + eigvecs_propagated[0].shape, (len(omega), len(n_opers), d, d), eigvecs_propagated[0].shape, optimize=[(0, 1), (0, 1)]) for g in util.progressbar_range(len(dt), show_progressbar=show_progressbar, desc='Calculating noise operators'): + if cache_intermediates: + # Assign references to the locations in the cache for the quantities + # that should be stored + sum_buf = sum_cache[g] + int_buf = _first_order_integral(omega, eigvals[g], dt[g], exp_buf, int_buf) - intermediate = expr_1(n_opers_transformed[:, g], - util.cexp(omega[:, None, None]*t[g])*int_buf, out=intermediate) + sum_buf = expr_1(n_opers_transformed[:, g], util.cexp(omega*t[g])[:, None, None]*int_buf, + out=sum_buf) - noise_operators += expr_2(eigvecs_propagated[g].conj(), intermediate, - eigvecs_propagated[g]) + noise_operators += expr_2(eigvecs_propagated[g].conj(), sum_buf, eigvecs_propagated[g], + out=sum_buf) + + if cache_intermediates: + intermediates = dict(n_opers_transformed=n_opers_transformed, + noise_operators_step=sum_cache) + return noise_operators, intermediates return noise_operators From ca66d257964787c9078a4bdd47983592a91d8563 Mon Sep 17 00:00:00 2001 From: Tobias Hangleiter Date: Tue, 2 Mar 2021 14:40:55 +0100 Subject: [PATCH 04/11] Refactor calculation of gradients Should improve performance by factor of ~2-5 in most cases. For now, there are two working solutions: one vectorized and one where one time axis is looped over, similar to the calculation of the filter function. --- filter_functions/gradient.py | 823 +++++++----- filter_functions/gradient_.py | 1855 ++++++++++++++++++++++++++++ filter_functions/pulse_sequence.py | 57 +- tests/gradient_testutil.py | 64 +- tests/test_gradient.py | 73 +- tests/test_precision.py | 2 +- 6 files changed, 2445 insertions(+), 429 deletions(-) create mode 100644 filter_functions/gradient_.py diff --git a/filter_functions/gradient.py b/filter_functions/gradient.py index ee48778..52a0faf 100644 --- a/filter_functions/gradient.py +++ b/filter_functions/gradient.py @@ -29,14 +29,18 @@ # Contact email: j.teske@fz-juelich.de # ============================================================================= """ -This module implements functions to calculate filter function and infidelity -derivatives. - -Throughout this documentation the following notation will be used: n_dt denotes -the number of time steps, n_cops the number of all control operators, n_ctrl -the number of accessible control operators (if identifiers are provided, -otherwise n_ctrl=n_cops), n_nops the number of noise operators, n_omega the -number of frequency samples and d the dimension of the system. +This module implements functions to calculate filter function and +infidelity derivatives. Currently only auto-correlated noise (i.e. no +cross-correlations) is implemented. + +Throughout this documentation the following notation will be used: + - n_dt denotes the number of time steps, + - n_cops the number of all control operators, + - n_ctrl the number of accessible control operators (if identifiers + are provided, otherwise n_ctrl=n_cops), + - n_nops the number of noise operators, + - n_omega the number of frequency samples, and + - d the dimension of the system. Functions --------- @@ -52,68 +56,96 @@ :func:`infidelity_derivative` Calculate the infidelity derivative. """ -from typing import Callable, Optional, Sequence, Union +from typing import Dict, Optional, Sequence, Tuple import numpy as np +import opt_einsum as oe from numpy import ndarray -from scipy.integrate import trapz -from filter_functions.basis import Basis -from filter_functions.types import Coefficients, Operator -from filter_functions.util import cexp +from . import numeric, superoperator, util +from .basis import Basis +from .types import Coefficients, Operator __all__ = ['liouville_derivative', 'control_matrix_at_timestep_derivative', 'calculate_derivative_of_control_matrix_from_scratch', - 'calculate_canonical_filter_function_derivative', - 'infidelity_derivative'] + 'calculate_filter_function_derivative', 'infidelity_derivative'] -def liouville_derivative( - dt: Coefficients, - propagators: ndarray, - basis: Basis, - eigvecs: ndarray, - eigvals: ndarray, - transf_control_operators: ndarray) -> ndarray: +def _derivative_integral(E: Coefficients, eigvals: Coefficients, dt: float, + out: ndarray) -> ndarray: + """ + Compute the integral appearing in the derivative of the control + matrix. Result (out) has shape (len(E), d, d, d, d). + """ + # Precompute masks and energy differences + dE = np.subtract.outer(eigvals, eigvals) + mask_dE = np.abs(dE) < 1e-7 + EdE = np.add.outer(E, dE) + mask_EdE = np.abs(EdE) < 1e-7 + EdEdE = np.add.outer(EdE, dE[~mask_dE]) + mask_EdEdE = np.abs(EdEdE) < 1e-7 + + # Case Omega_pq == 0 + tmp1 = np.divide(util.cexp(EdE*dt), EdE, where=~mask_EdE) + tmp2 = tmp1 - np.divide(1, EdE, where=~mask_EdE) + tmp2[mask_EdE] = 1j * dt + + tmp1 *= -1j * dt + tmp1 += np.divide(tmp2, EdE, where=~mask_EdE) + tmp1[mask_EdE] = dt**2 / 2 + + out[:, mask_dE] = tmp1[:, None] + + # Case Omega_pq != 0 + tmp1 = np.divide(1 - util.cexp(EdEdE*dt, where=~mask_EdEdE), EdEdE, where=~mask_EdEdE) + tmp1[mask_EdEdE] = -1j * dt + tmp1 += tmp2[..., None] + + out[:, ~mask_dE] = (tmp1 / dE[~mask_dE]).transpose(0, 3, 1, 2) + return out + + +def liouville_derivative(dt: Coefficients, propagators: ndarray, basis: Basis, eigvecs: ndarray, + eigvals: ndarray, c_opers_transformed: ndarray) -> ndarray: r""" Calculate the derivatives of the control propagators in Liouville representation. Parameters ---------- - 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}`. - propagators : array_like, shape (n_dt+1, d, d) - The propagators :math:`Q_g = P_g P_{g-1}\cdots P_0` as a (d, d) array - with *d* the dimension of the Hilbert space. - basis : Basis, shape (d**2, d, d) - The basis elements, in which the pulse control matrix will be expanded. - eigvecs : array_like, shape (n_dt, d, d) - Eigenvector matrices for each time pulse segment *g* with the first - axis counting the pulse segment, i.e. + propagators: array_like, shape (n_dt+1, d, d) + The propagators :math:`Q_g = P_g P_{g-1}\cdots P_0` as a (d, d) + array with *d* the dimension of the Hilbert space. + basis: Basis, shape (d**2, d, d) + The basis elements, in which the pulse control matrix will be + expanded. + eigvecs: array_like, shape (n_dt, d, d) + Eigenvector matrices for each time pulse segment *g* with the + first axis counting the pulse segment, i.e. ``HV == array([V_0, V_1, ...])``. - eigvals : array_like, shape (n_dt, d) - Eigenvalue vectors for each time pulse segment *g* with the first - axis counting the pulse segment, i.e. + eigvals: array_like, shape (n_dt, d) + Eigenvalue vectors for each time pulse segment *g* with the + first axis counting the pulse segment, i.e. ``HD == array([D_0, D_1, ...])``. - transf_control_operators : array_like, shape (n_dt, c_ctrl, d, d) - The control operators transformed into the eigenspace of the control - Hamiltonian. The drift operators are ignored, if identifiers for - accessible control operators are provided. + c_opers_transformed: array_like, shape (n_dt, c_ctrl, d, d) + The control operators transformed into the eigenspace of the + control Hamiltonian. The drift operators are ignored, if + identifiers for accessible control operators are provided. Returns ------- - dL: array_like, shape (n_dt, n_ctrl, n_dt, d**2, d**2) - The derivative of the control propagators in Liouville representation - :math:`\frac{\partial \mathcal{Q}_{jk}^{(g)}} + liouville_deriv: array_like, shape (n_dt, n_ctrl, n_dt, d**2, d**2) + The derivative of the control propagators in Liouville + representation :math:`\frac{\partial \mathcal{Q}_{jk}^{(g)}} {\partial u_h(t_{g^\prime})}`. - The array's indexing has shape :math:`(g,h,g^\prime,j,k)`. Notes ----- - The derivatives of the control propagators in Liouville representation are - calculated according to + The derivatives of the control propagators in Liouville + representation are calculated according to .. math:: @@ -130,46 +162,34 @@ def liouville_derivative( zero otherwise. """ - n = len(dt) - d = basis.shape[-1] - - # Allocate memory - A_mat = np.empty((d, d), dtype=complex) - partial_U = np.empty((n, transf_control_operators.shape[1], d, d), - dtype=complex) - deriv_U = np.zeros((n, n, transf_control_operators.shape[1], d, d), - dtype=complex) - - for g in (range(n)): - omega_diff = np.subtract.outer(eigvals[g], eigvals[g]) - mask = (abs(omega_diff) < 1e-10) - A_mat[mask] = dt[g] # if the integral diverges - A_mat[~mask] = (cexp(omega_diff[~mask]*dt[g]) - 1) \ - / (1j * omega_diff[~mask]) - # Calculate dU(t_g,t_{g-1})/du_h within one time step - partial_U[g] = -1j * np.einsum( - 'ia,ja,jk,hkl,kl,ml->him', propagators[g + 1], - propagators[g].conj(), eigvecs[g], transf_control_operators[g], - A_mat, eigvecs[g].conj(), optimize=['einsum_path', (3, 4), (0, 1), - (0, 3), (0, 1), (0, 1)]) - # Calculate the whole propagator derivative - for g_prime in range(g+1): # condition g' <= g-1 - deriv_U[g, g_prime] = np.einsum( - 'ij,kj,hkl,lm->him', propagators[g + 1], - propagators[g_prime + 1].conj(), partial_U[g_prime], - propagators[g_prime], optimize=['einsum_path', (0, 1), - (0, 1), (0, 1)]) - - # Now calculate derivative of Liouville representation - # Calculate traces first to save memory - sum1 = np.einsum( - 'tshba,jbc,tcd,kda->thsjk', deriv_U.conj(), basis, propagators[1:], - basis, optimize=['einsum_path', (1, 2), (0, 2), (0, 1)]) - sum2 = np.einsum( - 'tba,jbc,tshcd,kda->thsjk', propagators[1:].conj(), basis, deriv_U, - basis, optimize=['einsum_path', (0, 1), (0, 2), (0, 1)]) - liouville_deriv = sum1 + sum2 - + n, d = eigvecs.shape[:2] + # omega_i - omega_j + omega_diff = eigvals[:, :, None] - eigvals[:, None, :] + dt_broadcast = np.broadcast_to(dt[:, None, None], omega_diff.shape) + # mask = omega_diff == 0 + mask = np.broadcast_to(np.eye(d, dtype=bool), omega_diff.shape) + A_mat = np.empty(omega_diff.shape, dtype=complex) + A_mat[mask] = dt_broadcast[mask] + A_mat[~mask] = 1j*(1 - util.cexp(omega_diff[~mask]*dt_broadcast[~mask])) / omega_diff[~mask] + U_deriv = -1j * (propagators[1:] + @ propagators[:-1].conj().swapaxes(-1, -2) + @ eigvecs + @ (A_mat * c_opers_transformed.swapaxes(0, 1)) + @ eigvecs.conj().swapaxes(-1, -2)) + + # Calculate the whole propagator derivative + propagators_deriv = np.zeros((c_opers_transformed.shape[1], n-1, n, d, d), dtype=complex) + for g in range(n-1): + propagators_deriv[:, g, :g+1] = (propagators[g+1] + @ propagators[1:g+2].conj().swapaxes(-1, -2) + @ U_deriv[:, :g+1] + @ propagators[:g+1]) + + # Equivalent but usually slower contraction: 'htsba,jbc,tcd,kda->thsjk' + # Can just take 2*Re(·) when calculating x + x* + liouville_deriv = np.einsum('htsba,tjkba->thsjk', propagators_deriv.conj(), + (basis @ propagators[1:-1, None])[:, :, None] @ basis).real + liouville_deriv *= 2 return liouville_deriv @@ -179,62 +199,61 @@ def control_matrix_at_timestep_derivative( eigvals: ndarray, eigvecs: ndarray, basis: Basis, + c_opers_transformed: ndarray, n_opers: Sequence[Operator], n_coeffs: Sequence[Coefficients], - transf_control_operators: ndarray, - s_derivs: Optional[Sequence[Coefficients]] = None) -> dict: + n_coeffs_deriv: Optional[Sequence[Coefficients]] = None, + intermediates: Optional[Dict[str, ndarray]] = None +) -> Tuple[ndarray, ndarray]: r""" Calculate the control matrices and corresponding derivatives. - Calculate control matrices at each time step and the corresponding partial - derivatives of those with respect to control strength at each time step. + Calculate control matrices at each time step and the corresponding + partial derivatives of those with respect to control strength at + each time step. Parameters ---------- - omega : array_like, shape (n_omega) - Frequencies, at which the pulse control matrix is to be evaluated. - dt : array_like, shape (n_dt) + omega: array_like, shape (n_omega) + Frequencies, at which the pulse control matrix is to be + evaluated. + dt: array_like, shape (n_dt) Sequence duration, i.e. for the :math:`g`-th pulse :math:`t_g - t_{g-1}`. - eigvals : array_like, shape (n_dt, d) - Eigenvalue vectors for each time pulse segment *g* with the first - axis counting the pulse segment, i.e. + eigvals: array_like, shape (n_dt, d) + Eigenvalue vectors for each time pulse segment *g* with the + first axis counting the pulse segment, i.e. ``HD == array([D_0, D_1, ...])``. - eigvecs : array_like, shape (n_dt, d, d) - Eigenvector matrices for each time pulse segment *g* with the first - axis counting the pulse segment, i.e. + eigvecs: array_like, shape (n_dt, d, d) + Eigenvector matrices for each time pulse segment *g* with the + first axis counting the pulse segment, i.e. ``HV == array([V_0, V_1, ...])``. - basis : Basis, shape (d**2, d, d) - The basis elements, in which the pulse control matrix will be expanded. - n_opers : array_like, shape (n_nops, d, d) + basis: Basis, shape (d**2, d, d) + The basis elements, in which the pulse control matrix will be + expanded. + c_opers_transformed: array_like, shape (n_dt, n_ctrl, d, d) + The control operators transformed into the eigenspace of the + control Hamiltonian. The drift operators are ignored, if + identifiers for accessible control operators are provided. + n_opers: array_like, shape (n_nops, d, d) Noise operators :math:`B_\alpha`. - n_coeffs : array_like, shape (n_nops, n_dt) + 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. - transf_control_operators : array_like, shape (n_dt, n_ctrl, d, d) - The control operators transformed into the eigenspace of the control - Hamiltonian. The drift operators are ignored, if identifiers for - accessible control operators are provided. - s_derivs : array_like, shape (n_nops, n_ctrl, n_dt) - The derivatives of the noise susceptibilities by the control amplitudes. - Defaults to None. + n_coeffs_deriv: array_like, shape (n_nops, n_ctrl, n_dt) + The derivatives of the noise susceptibilities by the control + amplitudes. Defaults to None. + intermediates: Dict[str, ndarray], optional + Optional dictionary containing intermediate results of the + calculation of the control matrix. Returns ------- - ctrlmat_data : dict {'R_g': R_g, 'dR_g': gradient} - * **R_g** *(array_like, shape (n_dt, n_nops, d**2, n_omega))* - The control matrix at each time step - :math:`\mathcal{R}_{\alpha j}^{(g)}(\omega)` is identified with R_g. - The array's indexing has shape :math:`(g,\alpha,j,\omega)`. - - * **dR_g** *(array_like, shape (n_dt, n_nops, d**2, n_ctrl, n_omega))* - The corresponding derivative with respect to the control strength - :math:`\frac{\partial\mathcal{R}_{\alpha j}^{(g)}(\omega)} - {\partial u_h(t_{g^\prime})}` is identified with dR_g. The array's - indexing has shape :math:`(g^\prime,\alpha,j,h,\omega)`. Here, only one - time axis is needed, since the derivative is zero for - :math:`g\neq g^\prime`. - + ctrlmat_g: ndarray, shape (n_dt, n_nops, d**2, n_omega) + The individual control matrices of all time steps + ctrlmat_g_deriv: ndarray, shape (n_dt, n_nops, d**2, n_ctrl, n_omega) + The corresponding derivative with respect to the control + strength :math:`\frac{\partial\mathcal{B}_{\alpha j}^{(g)}(\omega)} Notes ----- @@ -242,7 +261,7 @@ def control_matrix_at_timestep_derivative( .. math:: - \mathcal{R}_{\alpha j}^{(g)}(\omega) = s_\alpha^{(g)}\mathrm{tr} + \mathcal{B}_{\alpha j}^{(g)}(\omega) = s_\alpha^{(g)}\mathrm{tr} \left([\bar{B}_\alpha \circ I_1^{(g)}(\omega)] \bar{C}_j \right), where @@ -253,12 +272,12 @@ def control_matrix_at_timestep_derivative( - \omega_m^{(g)}) \Delta t_g) - 1} {\mathrm{i}(\omega + \omega_n^{(g)} - \omega_m^{(g)})} - The derivative of the control matrix with respect to the control strength - at different time steps is calculated according to + The derivative of the control matrix with respect to the control + strength at different time steps is calculated according to .. math:: - \frac{\partial \mathcal{R}_{\alpha j}^{(g)}(\omega)} + \frac{\partial \mathcal{B}_{\alpha j}^{(g)}(\omega)} {\partial u_h(t_{g^\prime})} = \mathrm{i}\delta_{gg^\prime} s_\alpha^{(g)} \mathrm{tr} \left( \bar{B}_{\alpha} \cdot \mathbb{M} \right) @@ -266,11 +285,11 @@ def control_matrix_at_timestep_derivative( \left( (\overline{B}_{\alpha} \circ I_1^{(g)}{}(\omega)) \cdot \overline{C}_{j}) \right). - We assume that the noise susceptibility :math:`s` only depends locally - on the time i.e. :math:`\partial_{u(t_g)} s(t_{g^\prime}) + We assume that the noise susceptibility :math:`s` only depends + locally on the time i.e. :math:`\partial_{u(t_g)} s(t_{g^\prime}) = \delta_{gg^\prime} \partial_{u(t_g)} s(t_g)` - If denoting :math:`\Delta\omega_{ij} = \omega_i^{(g)} - \omega_j^{(g)}` - the integral part is encapsulated in + If denoting :math:`\Delta\omega_{ij} = \omega_i^{(g)} - + \omega_j^{(g)}` the integral part is encapsulated in .. math:: @@ -297,84 +316,149 @@ def control_matrix_at_timestep_derivative( """ d = eigvecs.shape[-1] n_dt = len(dt) - E = omega - - # Precompute some transformation into eigenspace of control Hamiltonian - path = ['einsum_path', (0, 1), (0, 1)] - VBV = np.einsum('gji,ajk,gkl->gail', eigvecs.conj(), n_opers, eigvecs, - optimize=path) - VCV = np.einsum('tnm,jnk,tkl->tjml', eigvecs.conj(), basis, eigvecs, - optimize=path) - - # Allocate memory - R_g = np.empty((n_dt, len(n_opers), len(basis), len(E)), dtype=complex) - R_g_d_s = np.empty( - (n_dt, len(n_opers), len(basis), len(transf_control_operators[0]), - len(E)), dtype=complex) - # For calculating dR_g: d,d = p,q, d,d = m,n - integral_deriv = np.empty((n_dt, len(E), d, d, d, d), dtype=complex) - + n_ctrl = len(c_opers_transformed[0]) + n_nops = len(n_opers) + n_omega = len(omega) + + # Cannot grab this from cached intermediates because those include the propagators + basis_transformed = numeric._transform_by_unitary(eigvecs[:, None], basis[None], + out=np.empty((n_dt, d**2, d, d), complex)) + if not intermediates: + n_opers_transformed = numeric._transform_hamiltonian(eigvecs, n_opers, + n_coeffs).swapaxes(0, 1) + exp_buf, integral = np.empty((2, n_omega, d, d), dtype=complex) + else: + n_opers_transformed = intermediates['n_opers_transformed'].swapaxes(0, 1) + first_order_integral = intermediates['first_order_integral'] + + deriv_integral = np.empty((n_dt, n_omega, d, d, d, d), dtype=complex) + ctrlmat_g = np.empty((n_dt, n_nops, d**2, n_omega), dtype=complex) + ctrlmat_g_deriv_s = np.empty((n_dt, n_nops, d**2, n_ctrl, n_omega), dtype=complex) + # Expression for calculating the control matrix during single time step + expr = oe.contract_expression('bcd,adc,hdc->abh', basis.shape, basis.shape, + (n_omega, d, d), optimize=[(0, 1), (0, 1)]) for g in range(n_dt): - # Any combination of \omega_m-\omega_n (R_g), \omega_p-\omega_q (dR_g) - dE = np.subtract.outer(eigvals[g], eigvals[g]) - # Any combination of \omega+\omega_m-\omega_n (R_g) or - # \omega-\omega_m+\omega_n (dR_g) - EdE = np.add.outer(E, dE) - - # 1) Calculation of the control matrix R_g at each time step - integral_Rg = np.empty((len(E), d, d), dtype=complex) - # Mask the integral to avoid convergence problems - mask_Rg = np.abs(EdE*dt[g]) <= 1e-7 - integral_Rg[mask_Rg] = dt[g] - integral_Rg[~mask_Rg] = (cexp(EdE[~mask_Rg]*dt[g]) - 1) \ - / (1j*(EdE[~mask_Rg])) - - R_g[g] = np.einsum('a,bcd,adc,hdc->abh', n_coeffs[:, g], VCV[g], - VBV[g], integral_Rg, - optimize=['einsum_path', (0, 2), (0, 2), (0, 1)]) - - # Add the dependency of the susceptibilities - # s_derivs has shape (n_nops, n_ctrl, n_dt) - # VCV has shape (n_dt, d**2, d, d) - # VBV has shape (n_dt, n_nops, d, d) - # integral_Rg has shape (n_omega, d, d) - # R_g_d_s shall be of shape (n_dt, n_nops, d**2, n_ctrl, n_omega) - if s_derivs is not None: - R_g_d_s[g] = np.einsum( - 'ae,bcd,adc,hdc->abeh', s_derivs[:, :, g], VCV[g], VBV[g], - integral_Rg, optimize=['einsum_path', (0, 2), (0, 2), (0, 1)]) - - # 2) Calculation of derivatives of control matrices at each time step - mask_deriv = (abs(dE) < 1e-15) # case: \omega_p-\omega_q = 0 - # calculation if omega_diff = 0 - n_case = sum(sum(mask_deriv)) - a = dt[g]*cexp(EdE*dt[g]) / (1j * EdE) \ - + (cexp(EdE*dt[g]) - 1) / (EdE)**2 - integral_deriv[g, :, mask_deriv] = np.concatenate(([[a]*n_case]), - axis=0) - # calculation if omega_diff != 0 - b1 = - (cexp(np.add.outer(EdE, dE[~mask_deriv])*dt[g]) - 1) \ - / (np.add.outer(EdE, dE[~mask_deriv])) / dE[~mask_deriv] - b2 = + np.divide.outer(((cexp(EdE*dt[g]) - 1) / EdE), dE[~mask_deriv]) - integral_deriv[g, :, ~mask_deriv] = (b1 + b2).transpose(3, 0, 1, 2) - - # Computation of the derivative/ gradient - I_circ_H = np.einsum('toijnm,thij->tohijnm', integral_deriv, - transf_control_operators) - M_mat = (np.einsum('tjmk,tohknnm->tojhmn', VCV, I_circ_H) - - np.einsum('tohmknm,tjkn->tojhmn', I_circ_H, VCV)) - gradient = 1j * np.einsum('at,tamn,tojhnm->tajho', n_coeffs, VBV, M_mat, - optimize=['einsum_path', (1, 2), (0, 1)]) - if s_derivs is not None: - gradient += R_g_d_s - ctrlmat_data = {'R_g': R_g, 'dR_g': gradient} - return ctrlmat_data + deriv_integral[g] = _derivative_integral(omega, eigvals[g], dt[g], deriv_integral[g]) + if not intermediates: + integral = numeric._first_order_integral(omega, eigvals[g], dt[g], exp_buf, integral) + else: + integral = first_order_integral[g] + + ctrlmat_g[g] = expr(basis_transformed[g], n_opers_transformed[g], integral, + out=ctrlmat_g[g]) + + if n_coeffs_deriv is not None: + # equivalent contraction: 'ah,a,ako->akho', but this faster + ctrlmat_g_deriv_s[g] = ( + (n_coeffs_deriv[..., g] / n_coeffs[:, g, None])[:, None, :, None] + * ctrlmat_g[g, :, :, None] + ) + + # Basically a tensor product + # K = np.einsum('taij,thkl->tahikjl', VBV, VHV) + # L = K.transpose(0, 1, 2, 4, 3, 6, 5) + K = util.tensor(n_opers_transformed[:, :, None], c_opers_transformed[:, None]) + L = util.tensor_transpose(K, (1, 0), [[d, d], [d, d]]) + k = np.diagonal(K.reshape(n_dt, n_nops, n_ctrl, d, d, d, d), 0, -2, -3) + l = np.diagonal(L.reshape(n_dt, n_nops, n_ctrl, d, d, d, d), 0, -2, -3) + i1 = np.diagonal(deriv_integral, 0, -2, -3) + i2 = np.diagonal(deriv_integral, 0, -1, -4) + # Reshaping in F-major is faster than not (~ factor 2-4) + M = np.einsum( + 'tahpm,topm->tahop', + l.reshape(n_dt, n_nops, n_ctrl, d**2, d, order='F'), + i1.reshape(n_dt, n_omega, d**2, d, order='F') + ).reshape(n_dt, n_nops, n_ctrl, n_omega, d, d, order='F') + if d == 2: + # Life a bit simpler + mask = np.eye(d, dtype=bool) + M[..., mask] -= M[..., mask][..., ::-1] + M[..., ~mask] *= 2 + else: + M -= np.einsum( + 'tahpn,topn->tahop', + k.swapaxes(-2, -3).reshape(n_dt, n_nops, n_ctrl, d**2, d, order='F'), + i2.reshape(n_dt, n_omega, d**2, d, order='F') + ).reshape(n_dt, n_nops, n_ctrl, n_omega, d, d, order='F').swapaxes(-1, -2) + + # Expand in basis transformed to eigenspace + ctrlmat_g_deriv = oe.contract('tjnk,tahokn->tajho', 1j*basis_transformed, M) + + if n_coeffs_deriv is not None: + ctrlmat_g_deriv += ctrlmat_g_deriv_s + + return ctrlmat_g, ctrlmat_g_deriv + + +def control_matrix_at_timestep_derivative_loop( + omega: Coefficients, + dt: Coefficients, + eigvals: ndarray, + eigvecs: ndarray, + basis_transformed, + c_opers_transformed: ndarray, + n_opers_transformed, + n_coeffs: Sequence[Coefficients], + n_coeffs_deriv: Sequence[Coefficients], + ctrlmat_step, + phase_factor, + deriv_integral, + integral, + ctrlmat_expr +) -> Tuple[ndarray, ndarray]: + """".""" + d = len(eigvecs) + d2 = d**2 + n_ctrl = len(c_opers_transformed) + n_nops = len(n_opers_transformed) + n_omega = len(omega) + + deriv_integral = _derivative_integral(omega, eigvals, dt, out=deriv_integral) + ctrlmat_step = ctrlmat_expr(phase_factor, basis_transformed, n_opers_transformed, integral, + out=ctrlmat_step) + + # Basically a tensor product + # K = np.einsum('taij,thkl->tahikjl', VBV, VHV) + # L = K.transpose(0, 1, 2, 4, 3, 6, 5) + K = util.tensor(n_opers_transformed[:, None], c_opers_transformed[None]) + L = util.tensor_transpose(K, (1, 0), [[d, d], [d, d]]) + k = np.diagonal(K.reshape(n_nops, n_ctrl, d, d, d, d), 0, -2, -3) + l = np.diagonal(L.reshape(n_nops, n_ctrl, d, d, d, d), 0, -2, -3) + i1 = np.diagonal(deriv_integral, 0, -2, -3) + i2 = np.diagonal(deriv_integral, 0, -1, -4) + # reshaping in F-major is significantly faster than C-major (~ factor 2-4) + M = np.einsum( + 'ahpm,opm->ahop', + l.reshape(n_nops, n_ctrl, d2, d, order='F'), + i1.reshape(n_omega, d2, d, order='F') + ).reshape(n_nops, n_ctrl, n_omega, d, d, order='F') + if d == 2: + # Life a bit simpler + mask = np.eye(d, dtype=bool) + M[..., mask] -= M[..., mask][..., ::-1] + M[..., ~mask] *= 2 + else: + M -= np.einsum( + 'ahpn,opn->ahop', + k.swapaxes(-2, -3).reshape(n_nops, n_ctrl, d2, d, order='F'), + i2.reshape(n_omega, d2, d, order='F') + ).reshape(n_nops, n_ctrl, n_omega, d, d, order='F').swapaxes(-1, -2) + + # Expand in basis transformed to eigenspace + ctrlmat_step_deriv = oe.contract('o,jnk,ahokn->ajho', phase_factor, 1j*basis_transformed, M, + optimize=[(1, 2), (0, 1)]) + + if n_coeffs_deriv is not None: + # equivalent contraction: 'ah,a,ako->akho', but this faster + ctrlmat_step_deriv += ((n_coeffs_deriv / n_coeffs[:, None])[:, None, :, None] + * ctrlmat_step[:, :, None]) + + return ctrlmat_step, ctrlmat_step_deriv def calculate_derivative_of_control_matrix_from_scratch( omega: Coefficients, propagators: ndarray, - Q_Liou: ndarray, eigvals: ndarray, eigvecs: ndarray, basis: Basis, @@ -385,64 +469,70 @@ def calculate_derivative_of_control_matrix_from_scratch( c_opers: Sequence[Operator], all_identifiers: Sequence[str], control_identifiers: Optional[Sequence[str]] = None, - s_derivs: Optional[Sequence[Coefficients]] = None) -> ndarray: + n_coeffs_deriv: Optional[Sequence[Coefficients]] = None, + intermediates: Optional[Dict[str, ndarray]] = None +) -> ndarray: r""" Calculate the derivative of the control matrix from scratch. Parameters ---------- - omega : array_like, shape (n_omega,) - Frequencies, at which the pulse control matrix is to be evaluated. - propagators : array_like, shape (n_dt+1, d, d) - The propagators :math:`Q_g = P_g P_{g-1}\cdots P_0` as a (d, d) array - with *d* the dimension of the Hilbert space. - Q_Liou : ndarray, shape (n_dt+1, d**2, d**2) - The Liouville representation of the cumulative control propagators - U_c(t_g,0). - eigvals : array_like, shape (n_dt, d) - Eigenvalue vectors for each time pulse segment *g* with the first - axis counting the pulse segment, i.e. + omega: array_like, shape (n_omega,) + Frequencies, at which the pulse control matrix is to be + evaluated. + propagators: array_like, shape (n_dt+1, d, d) + The propagators :math:`Q_g = P_g P_{g-1}\cdots P_0` as a (d, d) + array with *d* the dimension of the Hilbert space. + eigvals: array_like, shape (n_dt, d) + Eigenvalue vectors for each time pulse segment *g* with the + first axis counting the pulse segment, i.e. ``HD == array([D_0, D_1, ...])``. - eigvecs : array_like, shape (n_dt, d, d) - Eigenvector matrices for each time pulse segment *g* with the first - axis counting the pulse segment, i.e. + eigvecs: array_like, shape (n_dt, d, d) + Eigenvector matrices for each time pulse segment *g* with the + first axis counting the pulse segment, i.e. ``HV == array([V_0, V_1, ...])``. - basis : Basis, shape (d**2, d, d) - The basis elements, in which the pulse control matrix will be expanded. - t : array_like, shape (n_dt+1), optional + basis: Basis, shape (d**2, d, d) + The basis elements, in which the pulse control matrix will be + expanded. + t: array_like, shape (n_dt+1), optional The absolute times of the different segments. - 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}`. - n_opers : array_like, shape (n_nops, d, d) + n_opers: array_like, shape (n_nops, d, d) Noise operators :math:`B_\alpha`. - n_coeffs : array_like, shape (n_nops, n_dt) + 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. - c_opers : array_like, shape (n_cops, d, d) + c_opers: array_like, shape (n_cops, d, d) Control operators :math:`H_k`. - all_identifiers : array_like, shape (n_cops) + all_identifiers: array_like, shape (n_cops) Identifiers of all control operators. - control_identifiers : Sequence[str], shape (n_ctrl), Optional - Sequence of strings with the control identifiers to distinguish between - accessible control and drift Hamiltonian. The default is None. - s_derivs : array_like, shape (n_nops, n_ctrl, n_dt) - The derivatives of the noise susceptibilities by the control amplitudes. - Defaults to None. + control_identifiers: Sequence[str], shape (n_ctrl), Optional + Sequence of strings with the control identifiers to distinguish + between accessible control and drift Hamiltonian. The default is + None. + n_coeffs_deriv: array_like, shape (n_nops, n_ctrl, n_dt) + The derivatives of the noise susceptibilities by the control + amplitudes. Defaults to None. + intermediates: Dict[str, ndarray], optional + Optional dictionary containing intermediate results of the + calculation of the control matrix. Raises ------ ValueError - If the given identifiers *control_identifier* are not subset of the - total identifiers *all_identifiers* of all control operators. + If the given identifiers *control_identifier* are not subset of + the total identifiers *all_identifiers* of all control + operators. Returns ------- - dR : array_like, shape (n_nops, d**2, n_dt, n_ctrl, n_omega) - Partial derivatives of the total control matrix with respect to each - control direction - :math:`\frac{\partial R_{\alpha k}(\omega)}{\partial u_h(t_{g'})}`. - The array's indexing has shape :math:`(\alpha,k,g^\prime,h,\omega)`. + ctrlmat_deriv: ndarray, shape (n_ctrl, n_omega, n_dt, n_nops, d**2) + Partial derivatives of the total control matrix with respect to + each control direction + :math:`\frac{\partial\mathcal{B}_{\alpha k}(\omega)}{\partial + u_h(t_{g'})}`. Notes ----- @@ -450,11 +540,11 @@ def calculate_derivative_of_control_matrix_from_scratch( .. math :: - \frac{\partial R_{\alpha k}(\omega)}{\partial u_h(t_{g'})} = + \frac{\partial\mathcal{B}_{\alpha k}(\omega)}{\partial u_h(t_{g'})} = \sum_{g=1}^G \mathrm{e}^{\mathrm{i}\omega t_{g-1}}\cdot\left(\sum_j - \left[\frac{\partial R_{\alpha j}^{(g)}(\omega)} + \left[\frac{\partial\mathcal{B}_{\alpha j}^{(g)}(\omega)} {\partial u_h(t_{g'})} \cdot \mathcal{Q}_{jk}^{(g-1)} - + R_{\alpha j}^{(g)}(\omega) + + \mathcal{B}_{\alpha j}^{(g)}(\omega) \cdot\frac{\partial \mathcal{Q}_{jk}^{(g-1)}} {\partial u_h(t_{g'})} \right] \right) @@ -465,71 +555,133 @@ def calculate_derivative_of_control_matrix_from_scratch( """ # Distinction between control and drift operators and only calculate the # derivatives in control direction + try: + idx = util.get_indices_from_identifiers(all_identifiers, control_identifiers) + except ValueError as err: + raise ValueError('Given control identifiers have to be a subset of (drift+control) ' + + 'Hamiltonian!') from err + + c_opers_transformed = numeric._transform_hamiltonian(eigvecs, c_opers[idx]).swapaxes(0, 1) + propagators_liouville = superoperator.liouville_representation(propagators[:-1], basis) + propagators_liouville_deriv = liouville_derivative(dt, propagators, basis, eigvecs, eigvals, + c_opers_transformed) + + ctrlmat_g, ctrlmat_g_deriv = control_matrix_at_timestep_derivative( + omega, dt, eigvals, eigvecs, basis, c_opers_transformed, n_opers, n_coeffs, n_coeffs_deriv, + intermediates + ) - path = ['einsum_path', (0, 1), (0, 1)] - if (control_identifiers is None): - VHV = np.einsum('tji,hjk,tkl->thil', eigvecs.conj(), c_opers, eigvecs, - optimize=path) - elif (set(control_identifiers) <= set(all_identifiers)): - dict_id_oper = dict(zip(all_identifiers, c_opers)) - control = [dict_id_oper[element] for element in control_identifiers] - VHV = np.einsum('tji,hjk,tkl->thil', eigvecs.conj(), control, eigvecs, - optimize=path) + if intermediates: + phase_factors = intermediates['phase_factors'].T else: - raise ValueError('Given control identifiers have to be a \ - subset of (drift+control) Hamiltonian!') - # Get derivative of Liouville, control matrix at each time step, derivative - # derivative of control matrix at each time step - dL = liouville_derivative( - dt=dt, propagators=propagators, basis=basis, eigvecs=eigvecs, - eigvals=eigvals, transf_control_operators=VHV) - ctrlmat_data = control_matrix_at_timestep_derivative( - omega=omega, - dt=dt, - eigvals=eigvals, - eigvecs=eigvecs, - basis=basis, - n_opers=n_opers, - n_coeffs=n_coeffs, - transf_control_operators=VHV, - s_derivs=s_derivs - ) - ctrlmat_g = ctrlmat_data['R_g'] - ctrlmat_g_deriv = ctrlmat_data['dR_g'] + phase_factors = util.cexp(np.multiply.outer(omega, t[:-1])) - # Calculate the derivative of the total control matrix - exp_factor = cexp(np.multiply.outer(t, omega)) - summand1 = np.einsum('to,tajho,tjk->aktho', exp_factor[:-1], - ctrlmat_g_deriv, Q_Liou[:-1], - optimize=['einsum_path', (1, 2), (0, 1)]) - summand2 = np.einsum('to,tajo,thsjk->aksho', exp_factor[1:-1], - ctrlmat_g[1:], dL[:-1], - optimize=['einsum_path', (0, 1), (0, 1)]) + # Equivalent (but slower) einsum contraction for first term: + # ctrlmat_deriv = np.einsum('ot,tajho,tjk->hotak', + # phase_factors, ctrlmat_g_deriv, propagators_liouville, + # optimize=['einsum_path', (0, 2), (0, 1)]) + ctrlmat_deriv = (ctrlmat_g_deriv.transpose(3, 4, 0, 1, 2) + @ (phase_factors[..., None, None] * propagators_liouville)) + ctrlmat_deriv += np.einsum('ot,tajo,thsjk->hosak', + phase_factors[:, 1:], ctrlmat_g[1:], propagators_liouville_deriv, + optimize=['einsum_path', (0, 1), (0, 1)]) - dR = summand1 + summand2 - return dR + return ctrlmat_deriv -def calculate_canonical_filter_function_derivative( - R: ndarray, - deriv_R: ndarray) -> ndarray: +def calculate_derivative_of_control_matrix_from_scratch_loop( + omega: Coefficients, + propagators: ndarray, + eigvals: ndarray, + eigvecs: ndarray, + basis: Basis, + t: Coefficients, + dt: Coefficients, + n_opers: Sequence[Operator], + n_coeffs: Sequence[Coefficients], + c_opers: Sequence[Operator], + all_identifiers: Sequence[str], + control_identifiers: Optional[Sequence[str]] = None, + n_coeffs_deriv: Optional[Sequence[Coefficients]] = None, + intermediates: Dict[str, ndarray] = None +) -> ndarray: + r""".""" + # Distinction between control and drift operators and only + # calculate the derivatives in control direction + try: + idx = util.get_indices_from_identifiers(all_identifiers, control_identifiers) + except ValueError as err: + raise ValueError('Given control identifiers have to be a subset of (drift+control) ' + + 'Hamiltonian!') from err + + d = eigvecs.shape[-1] + n_dt = len(dt) + n_ctrl = len(idx) + n_nops = len(n_opers) + n_omega = len(omega) + + basis_transformed = numeric._transform_by_unitary(eigvecs[:, None], basis[None], + out=np.empty((n_dt, d**2, d, d), complex)) + c_opers_transformed = numeric._transform_hamiltonian(eigvecs, c_opers[idx]).swapaxes(0, 1) + if intermediates is None: + n_opers_transformed = numeric._transform_hamiltonian(eigvecs, n_opers, + n_coeffs).swapaxes(0, 1) + exp_buf, integral = np.empty((2, n_omega, d, d), dtype=complex) + else: + n_opers_transformed = intermediates['n_opers_transformed'].swapaxes(0, 1) + + propagators_liouville = superoperator.liouville_representation(propagators[:-1], basis) + propagators_liouville_deriv = liouville_derivative(dt, propagators, basis, eigvecs, eigvals, + c_opers_transformed) + + deriv_integral = np.empty((n_omega, d, d, d, d), dtype=complex) + ctrlmat_deriv = np.empty((n_ctrl, n_omega, n_dt, n_nops, d**2), dtype=complex) + ctrlmat_step = np.empty((n_dt, n_nops, d**2, n_omega), dtype=complex) + ctrlmat_expr = oe.contract_expression('o,icd,adc,odc->aio', (len(omega),), basis.shape, + n_opers.shape, (len(omega), d, d), + optimize=[(0, 3), (0, 1), (0, 1)]) + for g in range(n_dt): + if intermediates is None: + integral = numeric._first_order_integral(omega, eigvals[g], dt[g], exp_buf, integral) + else: + # Grab from cache + integral = intermediates['first_order_integral'][g] + + phase_factor = util.cexp(omega * t[g]) + n_coeff_deriv = n_coeffs_deriv if n_coeffs_deriv is None else n_coeffs_deriv[:, :, g] + + ctrlmat_step[g], ctrlmat_step_deriv = control_matrix_at_timestep_derivative_loop( + omega, dt[g], eigvals[g], eigvecs[g], basis_transformed[g], c_opers_transformed[g], + n_opers_transformed[g], n_coeffs[:, g], n_coeff_deriv, ctrlmat_step[g], + phase_factor, deriv_integral, integral, ctrlmat_expr + ) + ctrlmat_deriv[:, :, g] = (ctrlmat_step_deriv.transpose(2, 3, 0, 1) + @ propagators_liouville[g]) + + # opt_einsum a lot faster here + ctrlmat_deriv += oe.contract('tajo,thsjk->hosak', + ctrlmat_step[1:], propagators_liouville_deriv) + + return ctrlmat_deriv + + +def calculate_filter_function_derivative(ctrlmat: ndarray, ctrlmat_deriv: ndarray) -> ndarray: r""" Compute the filter function derivative from the control matrix. Parameters ---------- - R : array_like, shape (n_nops, d**2, n_omega) + ctrlmat: array_like, shape (n_nops, d**2, n_omega) The control matrix. - deriv_R: array_like, shape (n_nops, d**2, n_t, n_ctrl, n_omega) + ctrlmat_deriv: array_like, shape (n_nops, d**2, n_t, n_ctrl, n_omega) The derivative of the control matrix. Returns ------- - deriv_filter_function : ndarray, shape (n_nops, n_dt, n_ctrl, n_omega) + filter_function_deriv: ndarray, shape (n_nops, n_dt, n_ctrl, n_omega) The regular filter functions' derivatives for variation in each control contribution :math:`\frac{\partial F_\alpha(\omega)}{\partial u_h(t_{g'})}`. - The array's indexing has shape :math:`(\alpha,g^\prime,h,\omega)`. Notes ----- @@ -538,60 +690,61 @@ def calculate_canonical_filter_function_derivative( .. math :: \frac{\partial F_\alpha(\omega)}{\partial u_h(t_{g'})} - = 2 \mathrm{Re} \left( \sum_k R_{\alpha k}^\ast(\omega) - \frac{\partial R_{\alpha k}(\omega)} + = 2 \mathrm{Re}\left(\sum_k + \mathcal{B}_{\alpha k}^\ast(\omega) + \frac{\partial\mathcal{B}_{\alpha k}(\omega)} {\partial u_h(t_{g'})} \right) """ - summe = np.einsum('ako,aktho->atho', R.conj(), deriv_R) - return 2*summe.real + return 2*np.einsum('ako,hotak->atho', ctrlmat.conj(), ctrlmat_deriv).real def infidelity_derivative( pulse: 'PulseSequence', - S: Union[Coefficients, Callable], + spectrum: Coefficients, omega: Coefficients, control_identifiers: Optional[Sequence[str]] = None, - s_derivs: Optional[Sequence[Coefficients]] = None + n_coeffs_deriv: Optional[Sequence[Coefficients]] = None ) -> ndarray: - r""" - Calculate the infidelity derivative. - - Calculate the entanglement infidelity derivative of the ``PulseSequence`` - *pulse*. + r"""Calculate the entanglement infidelity derivative of the + ``PulseSequence`` *pulse*. Parameters ---------- - pulse : PulseSequence - The ``PulseSequence`` instance, for which to calculate the infidelity - for. - S : array_like + pulse: PulseSequence + The ``PulseSequence`` instance for which to calculate the + infidelity. + spectrum: array_like, shape ([[n_nops,] n_nops,] omega) The two-sided noise power spectral density in units of inverse - frequencies as an array of shape (n_omega,) or (n_nops, n_omega). In - the first case, the same spectrum is taken for all noise operators, in - the second, it is assumed that there are no correlations between - different noise sources and thus there is one spectrum for each noise - operator. - omega : array_like, shape (n_omega) + frequencies as an array of shape (n_omega,), (n_nops, n_omega), + or (n_nops, n_nops, n_omega). In the first case, the same + spectrum is taken for all noise operators, in the second, it is + assumed that there are no correlations between different noise + sources and thus there is one spectrum for each noise operator. + In the third and most general case, there may be a spectrum for + each pair of noise operators corresponding to the correlations + between them. n_nops is the number of noise operators considered + and should be equal to ``len(n_oper_identifiers)``. + omega: array_like, shape (n_omega,) The frequencies at which the integration is to be carried out. - control_identifiers : Sequence[str], shape (n_ctrl) - Sequence of strings with the control identifiern to distinguish between - accessible control and drift Hamiltonian. - s_derivs : array_like, shape (n_nops, n_ctrl, n_dt) - The derivatives of the noise susceptibilities by the control amplitudes. - Defaults to None. + control_identifiers: Sequence[str], shape (n_ctrl,) + Sequence of strings with the control identifiern to distinguish + between accessible control and drift Hamiltonian. + n_coeffs_deriv: array_like, shape (n_nops, n_ctrl, n_dt) + The derivatives of the noise susceptibilities by the control + amplitudes. Defaults to None. Raises ------ ValueError - If the provided noise spectral density does not fit expected shape. + If the provided noise spectral density does not fit expected + shape. Returns ------- - deriv_infid : ndarray, shape (n_nops, n_dt, n_ctrl) - Array with the derivative of the infidelity for each noise source taken - for each control direction at each time step - :math:`\frac{\partial I_e}{\partial u_h(t_{g'})}`. The array's indexing - has shape :math:`(\alpha,g^\prime,h)`. + infid_deriv: ndarray, shape (n_nops, n_dt, n_ctrl) + Array with the derivative of the infidelity for each noise + source taken for each control direction at each time step + :math:`\frac{\partial I_e}{\partial u_h(t_{g'})}`. Notes ----- @@ -603,16 +756,15 @@ def infidelity_derivative( \int_{-\infty}^\infty \frac{d\omega}{2\pi} S_\alpha(\omega) - \frac{\partial F_\alpha (\omega)} + \frac{\partial F_\alpha(\omega)} {\partial u_h(t_{g'})} with :math:`S_{\alpha}(\omega)` the noise spectral density and :math:`F_{\alpha}(\omega)` the canonical filter function for noise source :math:`\alpha`. - To convert to the average gate infidelity, use the - following relation given by Horodecki et al. [Hor99]_ and - Nielsen [Nie02]_: + To convert to the average gate infidelity, use the following + relation given by Horodecki et al. [Hor99]_ and Nielsen [Nie02]_: .. math:: @@ -633,22 +785,11 @@ def infidelity_derivative( Section A: General, Atomic and Solid State Physics, 303(4), 249–252. https://doi.org/10.1016/S0375-9601(02)01272-0 """ - n_ops = len(pulse.n_opers) - S = np.asarray(S) - if(S.shape[0] == n_ops): - S_all = S - elif(S.shape[0] == 1): - S_all = np.array([S[0]]*n_ops) - elif(S.shape[0] == len(omega)): - S_all = np.array([S]*n_ops) - else: - raise ValueError('Not fitting shape of S.') - - deriv_F = pulse.get_filter_function_derivative( - omega=omega, contorl_identifier=control_identifiers, s_derivs=s_derivs) - d = pulse.d - integrand = np.einsum('ao,atho->atho', S_all, deriv_F) + spectrum = numeric._parse_spectrum(spectrum, omega, range(len(pulse.n_opers))) + filter_function_deriv = pulse.get_filter_function_derivative(omega, control_identifiers, + n_coeffs_deriv) - deriv_infid = trapz(integrand, omega) / (2*np.pi*d) + integrand = np.einsum('ao,atho->atho', spectrum, filter_function_deriv) + infid_deriv = util.integrate(integrand, omega) / (2*np.pi*pulse.d) - return deriv_infid + return infid_deriv diff --git a/filter_functions/gradient_.py b/filter_functions/gradient_.py new file mode 100644 index 0000000..4ef2dfa --- /dev/null +++ b/filter_functions/gradient_.py @@ -0,0 +1,1855 @@ +# -*- coding: utf-8 -*- +# ============================================================================= +# filter_functions +# Copyright (C) 2019 Quantum Technology Group, RWTH Aachen University +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# This module is an extension of the filter_functions package written by +# Tobias Hangleiter. Its implementation was center of the Bachelor thesis +# 'Filter Function Derivative for Quantum Optimal Control' by Isabel Le. +# +# Contact email: isabel.le@rwth-aachen.de +# +# The code has been extended by Julian Teske such that the noise +# susceptibilities or noise coefficients can depend on the control +# amplitudes as well. +# +# Contact email: j.teske@fz-juelich.de +# ============================================================================= +""" +This module implements functions to calculate filter function and +infidelity derivatives. + +Throughout this documentation the following notation will be used: + - n_dt denotes the number of time steps, + - n_cops the number of all control operators, + - n_ctrl the number of accessible control operators (if identifiers + are provided, otherwise n_ctrl=n_cops), + - n_nops the number of noise operators, + - n_omega the number of frequency samples, and + - d the dimension of the system. + +Functions +--------- +:func:`liouville_derivative` + Calculate the derivatives of the control propagators in Liouville + representation. +:func:`control_matrix_at_timestep_derivative` + Calculate the control matrices and corresponding derivatives. +:func:`calculate_derivative_of_control_matrix_from_scratch` + Calculate the derivative of the control matrix from scratch. +:func:`calculate_canonical_filter_function_derivative` + Compute the filter function derivative from the control matrix. +:func:`infidelity_derivative` + Calculate the infidelity derivative. +""" +from typing import Callable, Dict, Optional, Sequence, Tuple, Union + +import numpy as np +import opt_einsum as oe +from numpy import ndarray +from scipy.integrate import trapz + +from . import numeric, superoperator, util +from .basis import Basis +from .types import Coefficients, Operator + +__all__ = ['liouville_derivative', 'control_matrix_at_timestep_derivative', + 'calculate_derivative_of_control_matrix_from_scratch', + 'calculate_filter_function_derivative', 'infidelity_derivative'] + + +def _derivative_integral_(E, eigvals, dt, deriv_integral): + """Calculate I_nmpq, but return I_pqnm""" + # Any combination of \omega_p-\omega_q (dR_g) + dE = np.subtract.outer(eigvals, eigvals) + mask_dE = np.abs(dE) < 1e-7 + # Any combination of \omega+\omega_n-\omega_m (dR_g) + EdE = np.add.outer(E, dE) + mask_EdE = np.abs(EdE) < 1e-7 + # Any combination of \omega+\omega_n-\omega_m+\omega_p-\omega_q (dR_g) + EdEdE = np.add.outer(EdE, dE[~mask_dE]) + + # calculation if omega_diff = 0 + tmp1 = np.divide(util.cexp(EdE*dt), EdE, where=~mask_EdE) + tmp2 = tmp1 - np.divide(1, EdE, where=~mask_EdE) + + tmp1 *= -1j * dt + tmp1 += np.divide(tmp2, EdE, where=~mask_EdE) + tmp1[mask_EdE] = dt**2 / 2 + + # tmp1 = util.cexp(EdE*dt) / EdE + # tmp2 = tmp1 - 1 / EdE + # tmp1 = -1j * dt * tmp1 + tmp2 / EdE + + deriv_integral[:, mask_dE] = tmp1[:, None] + + # # calculation if omega_diff != 0 + # tmp1 = np.divide(1 - util.cexp(EdEdE*dt), EdEdE, where=~mask_ndE_nEdE) + # tmp1[mask_ndE_EdE] = -1j * dt + # tmp1 += tmp2[:, None, None] + # tmp1 = np.divide(tmp1, dE, where=~mask_dE, out=tmp1) + # deriv_integral[:, ~mask_dE] = tmp1[:, ~mask_dE] + # calculation if omega_diff != 0 + + # pq on last axis in EdEdE, pq on 1,2nd axes in deriv_integral + deriv_integral[:, ~mask_dE] = ((1 - util.cexp(EdEdE*dt)) / EdEdE / dE[~mask_dE]).transpose(0, 3, 1, 2) + deriv_integral[:, ~mask_dE] += np.divide.outer(tmp2, dE[~mask_dE]).transpose(0, 3, 1, 2) + return deriv_integral + + +def _derivative_integral(E, eigvals, dt, out): + dE = np.subtract.outer(eigvals, eigvals) + mask_dE = np.abs(dE) < 1e-7 + EdE = np.add.outer(E, dE) + mask_EdE = np.abs(EdE) < 1e-7 + EdEdE = np.add.outer(EdE, dE[~mask_dE]) + mask_EdEdE = np.abs(EdEdE) < 1e-7 + + # Omega_pq == 0 + tmp1 = np.divide(util.cexp(EdE*dt), EdE, where=~mask_EdE) + tmp2 = tmp1 - np.divide(1, EdE, where=~mask_EdE) + tmp2[mask_EdE] = 1j * dt + + tmp1 *= -1j * dt + tmp1 += np.divide(tmp2, EdE, where=~mask_EdE) + tmp1[mask_EdE] = dt**2 / 2 + + out[:, mask_dE] = tmp1[:, None] + + # Omega_pq != 0 + tmp1 = np.divide(1 - util.cexp(EdEdE*dt), EdEdE, where=~mask_EdEdE) + tmp1[mask_EdEdE] = -1j * dt + tmp1 += tmp2[..., None] + + out[:, ~mask_dE] = (tmp1 / dE[~mask_dE]).transpose(0, 3, 1, 2) + + return out + + +def liouville_derivative( + dt: Coefficients, + propagators: ndarray, + basis: Basis, + eigvecs: ndarray, + eigvals: ndarray, + transf_control_operators: ndarray) -> ndarray: + r""" + Calculate the derivatives of the control propagators in Liouville + representation. + + Parameters + ---------- + dt : array_like, shape (n_dt) + Sequence duration, i.e. for the :math:`g`-th pulse + :math:`t_g - t_{g-1}`. + propagators : array_like, shape (n_dt+1, d, d) + The propagators :math:`Q_g = P_g P_{g-1}\cdots P_0` as a (d, d) array + with *d* the dimension of the Hilbert space. + basis : Basis, shape (d**2, d, d) + The basis elements, in which the pulse control matrix will be expanded. + eigvecs : array_like, shape (n_dt, d, d) + Eigenvector matrices for each time pulse segment *g* with the first + axis counting the pulse segment, i.e. + ``HV == array([V_0, V_1, ...])``. + eigvals : array_like, shape (n_dt, d) + Eigenvalue vectors for each time pulse segment *g* with the first + axis counting the pulse segment, i.e. + ``HD == array([D_0, D_1, ...])``. + transf_control_operators : array_like, shape (n_dt, c_ctrl, d, d) + The control operators transformed into the eigenspace of the control + Hamiltonian. The drift operators are ignored, if identifiers for + accessible control operators are provided. + + Returns + ------- + dL: array_like, shape (n_dt, n_ctrl, n_dt, d**2, d**2) + The derivative of the control propagators in Liouville representation + :math:`\frac{\partial \mathcal{Q}_{jk}^{(g)}} + {\partial u_h(t_{g^\prime})}`. + The array's indexing has shape :math:`(g,h,g^\prime,j,k)`. + + Notes + ----- + The derivatives of the control propagators in Liouville representation are + calculated according to + + .. math:: + + \frac{\partial\mathcal{Q}_{jk}^{(g-1)}}{\partial u_h(t_{g^\prime})} &= + \Theta_{g-1}(g^\prime) \mathrm{tr}\Big( + \frac{\partial U_c^\dagger(t_{g-1},0)}{\partial u_h(t_{g^\prime})} + C_j U_c(t_{g-1},0) C_k\Big)\\ + &+\Theta_{g-1}(g^\prime)\mathrm{tr}\Big(U_c^\dagger(t_{g-1},0)C_j + \frac{\partial U_c(t_{g-1},0)} + {\partial u_h(t_{g^\prime})}C_k + \Big), + + where :math:`\Theta_{g-1}(g^\prime)` being 1 if :math:`g^\prime < g-1` and + zero otherwise. + + """ + n = len(dt) + d = basis.shape[-1] + + # Allocate memory + A_mat = np.empty((d, d), dtype=complex) + partial_U = np.empty((n, transf_control_operators.shape[1], d, d), + dtype=complex) + deriv_U = np.zeros((n, n, transf_control_operators.shape[1], d, d), + dtype=complex) + + for g in (range(n)): + omega_diff = np.subtract.outer(eigvals[g], eigvals[g]) + mask = (abs(omega_diff) < 1e-10) + A_mat[mask] = dt[g] # if the integral diverges + A_mat[~mask] = (util.cexp(omega_diff[~mask]*dt[g]) - 1) \ + / (1j * omega_diff[~mask]) + # Calculate dU(t_g,t_{g-1})/du_h within one time step + partial_U[g] = -1j * np.einsum( + 'ia,ja,jk,hkl,kl,ml->him', propagators[g + 1], + propagators[g].conj(), eigvecs[g], transf_control_operators[g], + A_mat, eigvecs[g].conj(), optimize=['einsum_path', (3, 4), (0, 1), + (0, 3), (0, 1), (0, 1)]) + # Calculate the whole propagator derivative + for g_prime in range(g+1): # condition g' <= g-1 + deriv_U[g, g_prime] = np.einsum( + 'ij,kj,hkl,lm->him', propagators[g + 1], + propagators[g_prime + 1].conj(), partial_U[g_prime], + propagators[g_prime], optimize=['einsum_path', (0, 1), + (0, 1), (0, 1)]) + + # Now calculate derivative of Liouville representation + # Calculate traces first to save memory + sum1 = np.einsum( + 'tshba,jbc,tcd,kda->thsjk', deriv_U.conj(), basis, propagators[1:], + basis, optimize=['einsum_path', (1, 2), (0, 2), (0, 1)]) + sum2 = np.einsum( + 'tba,jbc,tshcd,kda->thsjk', propagators[1:].conj(), basis, deriv_U, + basis, optimize=['einsum_path', (0, 1), (0, 2), (0, 1)]) + liouville_deriv = sum1 + sum2 + + return liouville_deriv + + +def liouville_derivative_faster( + dt: Coefficients, + propagators: ndarray, + basis: Basis, + eigvecs: ndarray, + eigvals: ndarray, + transf_control_operators: ndarray) -> ndarray: + """.""" + n = len(dt) + d = basis.shape[-1] + + # Allocate memory + A_mat = np.empty((d, d), dtype=complex) + omega_diff = np.empty((d, d), dtype=float) + partial_U = np.empty((n, transf_control_operators.shape[1], d, d), + dtype=complex) + deriv_U = np.zeros((n, n, transf_control_operators.shape[1], d, d), + dtype=complex) + + for g in (range(n)): + omega_diff = np.subtract.outer(eigvals[g], eigvals[g], out=omega_diff) + mask = omega_diff == 0 + A_mat[mask] = dt[g] # if the integral diverges + A_mat[~mask] = (util.cexp(omega_diff[~mask]*dt[g]) - 1) / (1j * omega_diff[~mask]) + # Calculate dU(t_g,t_{g-1})/du_h within one time step + partial_U[g] = -1j * (propagators[g+1] @ propagators[g].conj().T @ eigvecs[g] @ + (transf_control_operators[g] * A_mat) @ eigvecs[g].conj().T) + # # Calculate the whole propagator derivative + deriv_U[g, :g+1] = (propagators[g+1] @ propagators[1:g+2].conj().swapaxes(-1, -2) @ + partial_U[:g+1].swapaxes(0, 1) @ propagators[:g+1]).swapaxes(0, 1) + + # Now calculate derivative of Liouville representation + # Calculate traces first to save memory + liouville_deriv = np.einsum( + 'tshba,jbc,tcd,kda->thsjk', deriv_U.conj(), basis, propagators[1:], + basis, optimize=['einsum_path', (1, 2), (0, 2), (0, 1)]) + + return 2*liouville_deriv.real + + +def liouville_derivative_vectorized( + dt: Coefficients, + propagators: ndarray, + basis: Basis, + eigvecs: ndarray, + eigvals: ndarray, + VHV: ndarray +) -> ndarray: + """".""" + n, d = eigvecs.shape[:2] + # omega_i - omega_j + omega_diff = eigvals[:, :, None] - eigvals[:, None] + dt_broadcast = np.broadcast_to(dt[:, None, None], omega_diff.shape) + # mask = omega_diff == 0 + mask = np.broadcast_to(np.eye(d, dtype=bool), omega_diff.shape) + A_mat = np.empty(omega_diff.shape, dtype=complex) + A_mat[mask] = dt_broadcast[mask] + # A_mat[~mask] = (util.cexp((omega_diff[~mask]*dt_broadcast[~mask])) - 1) / (1j * omega_diff[~mask]) + A_mat[~mask] = np.expm1(1j*(omega_diff[~mask]*dt_broadcast[~mask])) / (1j * omega_diff[~mask]) + U_deriv = -1j * (propagators[1:] + @ propagators[:-1].conj().swapaxes(-1, -2) + @ eigvecs + @ (A_mat * VHV.swapaxes(0, 1)) + @ eigvecs.conj().swapaxes(-1, -2)) + + # Calculate the whole propagator derivative + propagators_deriv = np.zeros((VHV.shape[1], n-1, n, d, d), dtype=complex) + for g in range(n-1): + propagators_deriv[:, g, :g+1] = (propagators[g+1] + @ propagators[1:g+2].conj().swapaxes(-1, -2) + @ U_deriv[:, :g+1] + @ propagators[:g+1]) + + # liouville_deriv = np.einsum('htsba,tjkba->thsjk', propagators_deriv.conj(), + # (basis @ propagators[1:-1, None])[:, :, None] @ basis).real + liouville_deriv = oe.contract('htsba,jbc,tcd,kda->thsjk', + propagators_deriv.conj(), basis, propagators[1:-1], basis, + optimize=[(0, 2), (0, 2), (0, 1)]).real + liouville_deriv *= 2 + + return liouville_deriv + + +def liouville_derivative_vectorized_completely( + dt: Coefficients, + propagators: ndarray, + basis: Basis, + eigvecs: ndarray, + eigvals: ndarray, + VHV: ndarray +) -> ndarray: + n, d = eigvecs.shape[:2] + # omega_i - omega_j + omega_diff = eigvals[:, :, None] - eigvals[:, None] + dt_broadcast = np.broadcast_to(dt[:, None, None], omega_diff.shape) + # mask = omega_diff == 0 + mask = np.broadcast_to(np.eye(d, dtype=bool), omega_diff.shape) + A_mat = np.empty(omega_diff.shape, dtype=complex) + A_mat[mask] = dt_broadcast[mask] + A_mat[~mask] = np.expm1(1j*(omega_diff[~mask]*dt_broadcast[~mask])) / (1j * omega_diff[~mask]) + U_deriv = -1j * (propagators[1:] + @ propagators[:-1].conj().swapaxes(-1, -2) + @ eigvecs + @ (A_mat * VHV.swapaxes(0, 1)) + @ eigvecs.conj().swapaxes(-1, -2)) + + # opt_einsum performs some magic with intermediate terms here, faster than explicit loop + propagators_deriv = oe.contract('tab,scb,hscd,sde->htsae', propagators[1:-1], + propagators[1:].conj(), U_deriv, propagators[:-1], + optimize=[(1, 2), (1, 2), (0, 1)]) + # Derivative is zero for times in the future of the control step + propagators_deriv[:, ~np.tri(n - 1, n, dtype=bool)] = 0 + + liouville_deriv = oe.contract('htsba,jbc,tcd,kda->thsjk', + propagators_deriv.conj(), basis, propagators[1:-1], basis, + optimize=[(0, 2), (0, 2), (0, 1)]).real + # oe.contract('xtba,jbc,tcd,kda->xtjk', propagators_deriv.transpose(0,2,1,3,4).reshape(5*100,99,4,4, order='F').conj(), basis, propagators[1:-1], basis, optimize=p).real.reshape(5,100,99,16,16,order='F') + liouville_deriv *= 2 + + return liouville_deriv + + +def liouville_derivative_vectorized_loop( + dt: Coefficients, + propagators: ndarray, + basis: Basis, + eigvecs: ndarray, + eigvals: ndarray, + c_opers_transformed: ndarray +) -> ndarray: + """".""" + n, d = eigvecs.shape[:2] + # omega_i - omega_j + omega_diff = eigvals[:, :, None] - eigvals[:, None] + dt_broadcast = np.broadcast_to(dt[:, None, None], omega_diff.shape) + # mask = omega_diff == 0 + mask = np.broadcast_to(np.eye(d, dtype=bool), omega_diff.shape) + integral = np.empty(omega_diff.shape, dtype=complex) + integral[mask] = dt_broadcast[mask] + # integral[~mask] = (util.cexp((omega_diff[~mask]*dt_broadcast[~mask])) - 1) / (1j * omega_diff[~mask]) + integral[~mask] = (np.expm1(1j*(omega_diff[~mask]*dt_broadcast[~mask])) + / (1j * omega_diff[~mask])) + + return -1j * (propagators[1:] + @ propagators[:-1].conj().swapaxes(-1, -2) + @ eigvecs + @ (integral * c_opers_transformed.swapaxes(0, 1)) + @ eigvecs.conj().swapaxes(-1, -2)) + + +def liouville_derivative_vectorized_completely_loop(propagator_deriv, basis, propagator, expr, + out): + liouville_deriv = expr(propagator_deriv.conj(), basis, propagator, basis, out=out) + liouville_deriv.real *= 2 + return liouville_deriv + + +def propagators_derivative_vectorized_completely( + dt: Coefficients, + propagators: ndarray, + basis: Basis, + eigvecs: ndarray, + eigvals: ndarray, + c_opers_transformed: ndarray +) -> ndarray: + + n, d = eigvecs.shape[:2] + # omega_i - omega_j + omega_diff = eigvals[:, :, None] - eigvals[:, None] + dt_broadcast = np.broadcast_to(dt[:, None, None], omega_diff.shape) + # mask = omega_diff == 0 + mask = np.broadcast_to(np.eye(d, dtype=bool), omega_diff.shape) + A_mat = np.empty(omega_diff.shape, dtype=complex) + A_mat[mask] = dt_broadcast[mask] + A_mat[~mask] = np.expm1(1j*(omega_diff[~mask]*dt_broadcast[~mask])) / (1j * omega_diff[~mask]) + U_deriv = -1j * (propagators[1:] + @ propagators[:-1].conj().swapaxes(-1, -2) + @ eigvecs + @ (A_mat * c_opers_transformed.swapaxes(0, 1)) + @ eigvecs.conj().swapaxes(-1, -2)) + + # opt_einsum performs some magic with intermediate terms here, faster than explicit loop + propagators_deriv = oe.contract('tab,scb,hscd,sde->htsae', propagators[1:-1], + propagators[1:].conj(), U_deriv, propagators[:-1], + optimize=[(1, 2), (1, 2), (0, 1)]) + # Derivative is zero for times in the future of the control step + propagators_deriv[:, ~np.tri(n - 1, n, dtype=bool)] = 0 + return propagators_deriv + + +def propagator_derivative( + dt: Coefficients, + propagators: ndarray, + eigvecs: ndarray, + eigvals: ndarray, + transf_control_operators: ndarray) -> ndarray: + """".""" + n, d = eigvecs.shape[:2] + # omega_i - omega_j + omega_diff = eigvals[:, :, None] - eigvals[:, None] + dt_broadcast = np.broadcast_to(dt[:, None, None], omega_diff.shape) + # mask = omega_diff == 0 + mask = np.broadcast_to(np.eye(d, dtype=bool), omega_diff.shape) + A_mat = np.empty(omega_diff.shape, dtype=complex) + A_mat[mask] = dt_broadcast[mask] + # A_mat[~mask] = (util.cexp((omega_diff[~mask]*dt_broadcast[~mask])) - 1) / (1j * omega_diff[~mask]) + A_mat[~mask] = np.expm1(1j*(omega_diff[~mask]*dt_broadcast[~mask])) / (1j * omega_diff[~mask]) + partial_U = -1j * (propagators[1:] @ propagators[:-1].conj().swapaxes(-1, -2) + @ eigvecs + @ (transf_control_operators.swapaxes(0, 1) * A_mat) + @ eigvecs.conj().swapaxes(-1, -2)) + + # Calculate the whole propagator derivative + derivative = np.zeros((transf_control_operators.shape[1], n, n, d, d), dtype=complex) + for g in range(n): + derivative[:, g, :g+1] = (propagators[g+1] @ propagators[1:g+2].conj().swapaxes(-1, -2) @ + partial_U[:, :g+1] @ propagators[:g+1]) + + return derivative + + +def propagator_derivative_factor( + dt: Coefficients, + propagators: ndarray, + eigvecs: ndarray, + eigvals: ndarray, + VHV: ndarray) -> ndarray: + """".""" + n, d = eigvecs.shape[:2] + # omega_i - omega_j + omega_diff = eigvals[:, :, None] - eigvals[:, None] + dt_broadcast = np.broadcast_to(dt[:, None, None], omega_diff.shape) + # mask = omega_diff == 0 + mask = np.broadcast_to(np.eye(d, dtype=bool), omega_diff.shape) + A_mat = np.empty(omega_diff.shape, dtype=complex) + A_mat[mask] = dt_broadcast[mask] + # A_mat[~mask] = (util.cexp((omega_diff[~mask]*dt_broadcast[~mask])) - 1) / (1j * omega_diff[~mask]) + A_mat[~mask] = np.expm1(1j*(omega_diff[~mask]*dt_broadcast[~mask])) / (1j * omega_diff[~mask]) + partial_U = -1j * (propagators[1:] @ propagators[:-1].conj().swapaxes(-1, -2) + @ eigvecs + @ (VHV.swapaxes(0, 1) * A_mat) + @ eigvecs.conj().swapaxes(-1, -2)) + + # Calculate the whole propagator derivative + derivative = np.zeros((VHV.shape[1], n, n, d, d), dtype=complex) + for g in range(n): + derivative[:, g, :g+1] = (propagators[1:g+2].conj().swapaxes(-1, -2) + @ partial_U[:, :g+1] + @ propagators[:g+1]) + + return derivative + + +def control_matrix_at_timestep_derivative( + omega: Coefficients, + dt: Coefficients, + eigvals: ndarray, + eigvecs: ndarray, + basis: Basis, + n_opers: Sequence[Operator], + n_coeffs: Sequence[Coefficients], + transf_control_operators: ndarray, + s_derivs: Optional[Sequence[Coefficients]] = None) -> dict: + r""" + Calculate the control matrices and corresponding derivatives. + + Calculate control matrices at each time step and the corresponding partial + derivatives of those with respect to control strength at each time step. + + Parameters + ---------- + omega : array_like, shape (n_omega) + Frequencies, at which the pulse control matrix is to be evaluated. + dt : array_like, shape (n_dt) + Sequence duration, i.e. for the :math:`g`-th pulse + :math:`t_g - t_{g-1}`. + eigvals : array_like, shape (n_dt, d) + Eigenvalue vectors for each time pulse segment *g* with the first + axis counting the pulse segment, i.e. + ``HD == array([D_0, D_1, ...])``. + eigvecs : array_like, shape (n_dt, d, d) + Eigenvector matrices for each time pulse segment *g* with the first + axis counting the pulse segment, i.e. + ``HV == array([V_0, V_1, ...])``. + basis : Basis, shape (d**2, d, d) + The basis elements, in which the pulse control matrix will be expanded. + n_opers : array_like, shape (n_nops, d, d) + Noise operators :math:`B_\alpha`. + 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. + transf_control_operators : array_like, shape (n_dt, n_ctrl, d, d) + The control operators transformed into the eigenspace of the control + Hamiltonian. The drift operators are ignored, if identifiers for + accessible control operators are provided. + s_derivs : array_like, shape (n_nops, n_ctrl, n_dt) + The derivatives of the noise susceptibilities by the control amplitudes. + Defaults to None. + + Returns + ------- + ctrlmat_data : dict {'R_g': R_g, 'dR_g': gradient} + * **R_g** *(array_like, shape (n_dt, n_nops, d**2, n_omega))* + The control matrix at each time step + :math:`\mathcal{R}_{\alpha j}^{(g)}(\omega)` is identified with R_g. + The array's indexing has shape :math:`(g,\alpha,j,\omega)`. + + * **dR_g** *(array_like, shape (n_dt, n_nops, d**2, n_ctrl, n_omega))* + The corresponding derivative with respect to the control strength + :math:`\frac{\partial\mathcal{R}_{\alpha j}^{(g)}(\omega)} + {\partial u_h(t_{g^\prime})}` is identified with dR_g. The array's + indexing has shape :math:`(g^\prime,\alpha,j,h,\omega)`. Here, only one + time axis is needed, since the derivative is zero for + :math:`g\neq g^\prime`. + + + Notes + ----- + The control matrix at each time step is evaluated according to + + .. math:: + + \mathcal{R}_{\alpha j}^{(g)}(\omega) = s_\alpha^{(g)}\mathrm{tr} + \left([\bar{B}_\alpha \circ I_1^{(g)}(\omega)] \bar{C}_j \right), + + where + + .. math:: + + I_{1,nm}^{(g)}(\omega) = \frac{\exp(\mathrm{i}(\omega + \omega_n^{(g)} + - \omega_m^{(g)}) \Delta t_g) - 1} + {\mathrm{i}(\omega + \omega_n^{(g)} - \omega_m^{(g)})} + + The derivative of the control matrix with respect to the control strength + at different time steps is calculated according to + + .. math:: + + \frac{\partial \mathcal{R}_{\alpha j}^{(g)}(\omega)} + {\partial u_h(t_{g^\prime})} = + \mathrm{i}\delta_{gg^\prime} s_\alpha^{(g)} \mathrm{tr} + \left( \bar{B}_{\alpha} \cdot \mathbb{M} \right) + + \frac{\partial s_\alpha^{(g)}}{u_h (t_{g^\prime})} \text{tr} + \left( (\overline{B}_{\alpha} \circ I_1^{(g)}{}(\omega)) \cdot + \overline{C}_{j}) \right). + + We assume that the noise susceptibility :math:`s` only depends locally + on the time i.e. :math:`\partial_{u(t_g)} s(t_{g^\prime}) + = \delta_{gg^\prime} \partial_{u(t_g)} s(t_g)` + If denoting :math:`\Delta\omega_{ij} = \omega_i^{(g)} - \omega_j^{(g)}` + the integral part is encapsulated in + + .. math:: + + \mathbb{M}_{mn} = \left[ \bar{C}_j, \mathbb{I}^{(mn)} + \circ \bar{H}_h \right]_{mn}, + + with + + .. math:: + + \mathbb{I}_{pq}^{(mn)} &= \delta_{pq} \left( + \frac{\Delta t_g \cdot + \exp(\mathrm{i}(\omega + \Delta\omega_{nm})\Delta t_g)} + {\mathrm{i}(\omega + \Delta\omega_{nm})} + + \frac{\exp(\mathrm{i}(\omega + \Delta\omega_{nm})\Delta t_g) - 1} + {(\omega + \Delta\omega_{nm})^2}\right)\\ + &+ \frac{1-\delta_{pq}}{\mathrm{i}\Delta\omega_{pq}} \cdot + \frac{\exp(\mathrm{i}(\omega + \Delta\omega_{nm} + + \Delta\omega_{pq})\Delta t_g) - 1} + {\mathrm{i}(\omega + \Delta\omega_{nm} + \Delta\omega_{pq})}\\ + &- \frac{1-\delta_{pq}}{\mathrm{i}\Delta\omega_{pq}} \cdot + \frac{\exp(\mathrm{i}(\omega + \Delta\omega_{nm})\Delta t_g) - 1} + {\mathrm{i}(\omega + \Delta\omega_{nm})} + """ + d = eigvecs.shape[-1] + n_dt = len(dt) + E = omega + + # Precompute some transformation into eigenspace of control Hamiltonian + path = ['einsum_path', (0, 1), (0, 1)] + VBV = np.einsum('gji,ajk,gkl->gail', eigvecs.conj(), n_opers, eigvecs, + optimize=path) + VCV = np.einsum('tnm,jnk,tkl->tjml', eigvecs.conj(), basis, eigvecs, + optimize=path) + + # Allocate memory + R_g = np.empty((n_dt, len(n_opers), len(basis), len(E)), dtype=complex) + R_g_d_s = np.empty( + (n_dt, len(n_opers), len(basis), len(transf_control_operators[0]), + len(E)), dtype=complex) + # For calculating dR_g: d,d = p,q, d,d = m,n + integral_deriv = np.empty((n_dt, len(E), d, d, d, d), dtype=complex) + + for g in range(n_dt): + # Any combination of \omega_m-\omega_n (R_g), \omega_p-\omega_q (dR_g) + dE = np.subtract.outer(eigvals[g], eigvals[g]) + # Any combination of \omega+\omega_m-\omega_n (R_g) or + # \omega-\omega_m+\omega_n (dR_g) + EdE = np.add.outer(E, dE) + + # 1) Calculation of the control matrix R_g at each time step + integral_Rg = np.empty((len(E), d, d), dtype=complex) + # Mask the integral to avoid convergence problems + mask_Rg = np.abs(EdE*dt[g]) <= 1e-7 + integral_Rg[mask_Rg] = dt[g] + integral_Rg[~mask_Rg] = (util.cexp(EdE[~mask_Rg]*dt[g]) - 1) \ + / (1j*(EdE[~mask_Rg])) + + R_g[g] = np.einsum('a,bcd,adc,hdc->abh', n_coeffs[:, g], VCV[g], + VBV[g], integral_Rg, + optimize=['einsum_path', (0, 2), (0, 2), (0, 1)]) + + # Add the dependency of the susceptibilities + # s_derivs has shape (n_nops, n_ctrl, n_dt) + # VCV has shape (n_dt, d**2, d, d) + # VBV has shape (n_dt, n_nops, d, d) + # integral_Rg has shape (n_omega, d, d) + # R_g_d_s shall be of shape (n_dt, n_nops, d**2, n_ctrl, n_omega) + if s_derivs is not None: + R_g_d_s[g] = np.einsum( + 'ae,bcd,adc,hdc->abeh', s_derivs[:, :, g], VCV[g], VBV[g], + integral_Rg, optimize=['einsum_path', (0, 2), (0, 2), (0, 1)]) + + # 2) Calculation of derivatives of control matrices at each time step + mask_deriv = (abs(dE) < 1e-15) # case: \omega_p-\omega_q = 0 + # calculation if omega_diff = 0 + n_case = sum(sum(mask_deriv)) + a = dt[g]*util.cexp(EdE*dt[g]) / (1j * EdE) \ + + (util.cexp(EdE*dt[g]) - 1) / (EdE)**2 + integral_deriv[g, :, mask_deriv] = np.concatenate(([[a]*n_case]), + axis=0) + # calculation if omega_diff != 0 + b1 = - (util.cexp(np.add.outer(EdE, dE[~mask_deriv])*dt[g]) - 1) \ + / (np.add.outer(EdE, dE[~mask_deriv])) / dE[~mask_deriv] + b2 = + np.divide.outer(((util.cexp(EdE*dt[g]) - 1) / EdE), dE[~mask_deriv]) + integral_deriv[g, :, ~mask_deriv] = (b1 + b2).transpose(3, 0, 1, 2) + + # Computation of the derivative/ gradient + I_circ_H = np.einsum('toijnm,thij->tohijnm', integral_deriv, + transf_control_operators) + M_mat = (np.einsum('tjmk,tohknnm->tojhmn', VCV, I_circ_H) + - np.einsum('tohmknm,tjkn->tojhmn', I_circ_H, VCV)) + gradient = 1j * np.einsum('at,tamn,tojhnm->tajho', n_coeffs, VBV, M_mat, + optimize=['einsum_path', (1, 2), (0, 1)]) + if s_derivs is not None: + gradient += R_g_d_s + ctrlmat_data = {'R_g': R_g, 'dR_g': gradient} + return ctrlmat_data + + +def control_matrix_at_timestep_derivative_refactored( + omega: Coefficients, + dt: Coefficients, + eigvals: ndarray, + eigvecs: ndarray, + basis: Basis, + c_opers_transformed: ndarray, + n_opers: Sequence[Operator], + n_coeffs: Sequence[Coefficients], + n_coeffs_deriv: Optional[Sequence[Coefficients]] = None, + intermediates: Optional[Dict[str, ndarray]] = None +) -> Tuple[ndarray, ndarray]: + """".""" + d = eigvecs.shape[-1] + n_dt = len(dt) + n_ctrl = len(c_opers_transformed[0]) + n_nops = len(n_opers) + E = omega + + # Precompute some transformation into eigenspace of control Hamiltonian + path = ['einsum_path', (0, 1), (0, 1)] + VBV = np.einsum('gji,ajk,gkl->gail', eigvecs.conj(), n_opers, eigvecs, optimize=path) + VCV = np.einsum('tnm,jnk,tkl->tjml', eigvecs.conj(), basis, eigvecs, optimize=path) + + # Allocate memory + R_g = np.empty((n_dt, n_nops, d**2, len(E)), dtype=complex) + buffer = np.empty((n_nops, d**2, len(E)), dtype=complex) + R_g_d_s = np.empty((n_dt, n_nops, d**2, n_ctrl, len(E)), dtype=complex) + # For calculating dR_g: d,d = p,q, d,d = m,n + integral_deriv = np.empty((n_dt, len(E), d, d, d, d), dtype=complex) + exp_buf, integral_Rg = np.empty((2, len(E), d, d), dtype=complex) + + expr = oe.contract_expression('bcd,adc,hdc->abh', VCV[0].shape, VBV[0].shape, + integral_Rg.shape, optimize=[(0, 1), (0, 1)]) + for g in range(n_dt): + integral_Rg = numeric._first_order_integral(E, eigvals[g], dt[g], exp_buf, integral_Rg) + + buffer = expr(VCV[g], VBV[g], integral_Rg, out=buffer) + R_g[g] = n_coeffs[:, g, None, None] * buffer + + # Add the dependency of the susceptibilities + # n_coeffs_deriv has shape (n_nops, n_ctrl, n_dt) + # VCV has shape (n_dt, d**2, d, d) + # VBV has shape (n_dt, n_nops, d, d) + # integral_Rg has shape (n_omega, d, d) + # R_g_d_s shall be of shape (n_dt, n_nops, d**2, n_ctrl, n_omega) + if n_coeffs_deriv is not None: + R_g_d_s[g] = n_coeffs_deriv[:, None, :, g, None] * buffer[:, :, None] + + integral_deriv[g] = _derivative_integral(E, eigvals[g], dt[g], integral_deriv[g]) + + # Computation of the derivative/ gradient + # I_circ_H = np.einsum('toijnm,thij->tohijnm', integral_deriv, c_opers_transformed) + I_circ_H = integral_deriv[:, :, None] * c_opers_transformed[:, None, ..., None, None] + if d == 2: + # Life a bit simpler + mask = np.eye(d, dtype=bool) + M_mat = np.zeros((n_dt, len(E), d**2, n_ctrl, d, d), dtype=complex) + M_mat[:, :, 1:] = np.einsum('tjmk,tohknnm->tojhmn', VCV[:, 1:], I_circ_H) + M_mat[:, :, 1:, :, mask] -= M_mat[:, :, 1:, :, mask][..., ::-1] + M_mat[:, :, 1:, :, ~mask] *= 2 + # M_mat -= M_mat.conj().swapaxes(-1, -2) + else: + M_mat = np.einsum('tjmk,tohknnm->tojhmn', VCV, I_circ_H) + M_mat -= np.einsum('tjkn,tohmknm->tojhmn', VCV, I_circ_H) + + ctrlmat_g_deriv = 1j*np.einsum('tamn,tojhnm->tajho', n_coeffs.T[..., None, None] * VBV, M_mat) + + if n_coeffs_deriv is not None: + ctrlmat_g_deriv += R_g_d_s + + return R_g, ctrlmat_g_deriv + + +def control_matrix_at_timestep_derivative_faster( + omega: Coefficients, + dt: Coefficients, + eigvals: ndarray, + eigvecs: ndarray, + basis: Basis, + c_opers_transformed: ndarray, + n_opers: Sequence[Operator], + n_coeffs: Sequence[Coefficients], + n_coeffs_deriv: Optional[Sequence[Coefficients]] = None +) -> Tuple[ndarray, ndarray]: + """".""" + d = eigvecs.shape[-1] + d2 = d**2 + n_dt = len(dt) + n_ctrl = len(c_opers_transformed[0]) + n_nops = len(n_opers) + n_omega = len(omega) + E = omega + + # Precompute some transformation into eigenspace of control Hamiltonian + n_opers_transformed = numeric._transform_hamiltonian(eigvecs, n_opers, n_coeffs).swapaxes(0, 1) + basis_transformed = numeric._transform_by_unitary(eigvecs[:, None], basis[None], + out=np.empty((n_dt, d2, d, d), complex)) + + # Allocate memory + ctrlmat_g = np.empty((n_dt, n_nops, d2, len(E)), dtype=complex) + ctrlmat_g_deriv_s = np.empty((n_dt, n_nops, d2, n_ctrl, len(E)), dtype=complex) + # For calculating dR_g: d,d = p,q, d,d = m,n + deriv_integral = np.empty((n_dt, len(E), d, d, d, d), dtype=complex) + exp_buf, integral = np.empty((2, len(E), d, d), dtype=complex) + + expr = oe.contract_expression('bcd,adc,hdc->abh', basis.shape, basis.shape, + integral.shape, optimize=[(0, 1), (0, 1)]) + for g in range(n_dt): + integral = numeric._first_order_integral(E, eigvals[g], dt[g], exp_buf, integral) + deriv_integral[g] = _derivative_integral(E, eigvals[g], dt[g], deriv_integral[g]) + + ctrlmat_g[g] = expr(basis_transformed[g], n_opers_transformed[g], integral, + out=ctrlmat_g[g]) + if n_coeffs_deriv is not None: + # equivalent contraction: 'ah,a,ako->akho', but this faster + ctrlmat_g_deriv_s[g] = ( + (n_coeffs_deriv[..., g] / n_coeffs[:, g, None])[:, None, :, None] + * ctrlmat_g[g, :, :, None] + ) + + # Basically a tensor product + # K = np.einsum('taij,thkl->tahikjl', VBV, VHV) + # L = K.transpose(0, 1, 2, 4, 3, 6, 5) + K = util.tensor(n_opers_transformed[:, :, None], c_opers_transformed[:, None]) + L = util.tensor_transpose(K, (1, 0), [[d, d], [d, d]]) + k = np.diagonal(K.reshape(n_dt, n_nops, n_ctrl, d, d, d, d), 0, -2, -3) + l = np.diagonal(L.reshape(n_dt, n_nops, n_ctrl, d, d, d, d), 0, -2, -3) + i1 = np.diagonal(deriv_integral, 0, -2, -3) + i2 = np.diagonal(deriv_integral, 0, -1, -4) + # Let magic ensue. Reshaping in F-major is faster than not (~ factor 2-4) + M = np.einsum( + 'tahpm,topm->tahop', + l.reshape(n_dt, n_nops, n_ctrl, d2, d, order='F'), + i1.reshape(n_dt, n_omega, d2, d, order='F') + ).reshape(n_dt, n_nops, n_ctrl, n_omega, d, d, order='F') + if d == 2: + # Life a bit simpler + mask = np.eye(d, dtype=bool) + M[..., mask] -= M[..., mask][..., ::-1] + M[..., ~mask] *= 2 + else: + M -= np.einsum( + 'tahpn,topn->tahop', + k.swapaxes(-2, -3).reshape(n_dt, n_nops, n_ctrl, d2, d, order='F'), + i2.reshape(n_dt, n_omega, d2, d, order='F') + ).reshape(n_dt, n_nops, n_ctrl, n_omega, d, d, order='F').swapaxes(-1, -2) + + # Expand in basis transformed to eigenspace + ctrlmat_g_deriv = np.einsum('tjnk,tahokn->tajho', 1j*basis_transformed, M) + + if n_coeffs_deriv is not None: + ctrlmat_g_deriv += ctrlmat_g_deriv_s + + return ctrlmat_g, ctrlmat_g_deriv + + +def _control_matrix_at_timestep_derivative_loop_1( + omega: Coefficients, + dt: Coefficients, + eigvals: ndarray, + eigvecs: ndarray, + basis_transformed, + c_opers_transformed: ndarray, + n_opers_transformed, + n_coeffs: Sequence[Coefficients], + n_coeffs_deriv: Sequence[Coefficients], + ctrlmat_step, + deriv_integral, + exp_buf, + integral, + expr +) -> Tuple[ndarray, ndarray]: + """".""" + d = len(eigvecs) + d2 = d**2 + n_ctrl = len(c_opers_transformed) + n_nops = len(n_opers_transformed) + n_omega = len(omega) + + integral = numeric._first_order_integral(omega, eigvals, dt, exp_buf, integral) + deriv_integral = _derivative_integral(omega, eigvals, dt, deriv_integral) + ctrlmat_step = expr(basis_transformed, n_opers_transformed, integral, out=ctrlmat_step) + if n_coeffs_deriv is not None: + # equivalent contraction: 'ah,a,ako->akho', but this faster + ctrlmat_step_deriv_s = ((n_coeffs_deriv / n_coeffs[:, None])[:, None, :, None] + * ctrlmat_step[:, :, None]) + + # Basically a tensor product + # K = np.einsum('taij,thkl->tahikjl', VBV, VHV) + # L = K.transpose(0, 1, 2, 4, 3, 6, 5) + K = util.tensor(n_opers_transformed[:, None], c_opers_transformed[None]) + L = util.tensor_transpose(K, (1, 0), [[d, d], [d, d]]) + k = np.diagonal(K.reshape(n_nops, n_ctrl, d, d, d, d), 0, -2, -3) + l = np.diagonal(L.reshape(n_nops, n_ctrl, d, d, d, d), 0, -2, -3) + i1 = np.diagonal(deriv_integral, 0, -2, -3) + i2 = np.diagonal(deriv_integral, 0, -1, -4) + # Let magic ensue. Reshaping in F-major is faster than not (~ factor 2-4) + M = np.einsum( + 'ahpm,opm->ahop', + l.reshape(n_nops, n_ctrl, d2, d, order='F'), + i1.reshape(n_omega, d2, d, order='F') + ).reshape(n_nops, n_ctrl, n_omega, d, d, order='F') + if d == 2: + # Life a bit simpler + mask = np.eye(d, dtype=bool) + M[..., mask] -= M[..., mask][..., ::-1] + M[..., ~mask] *= 2 + else: + M -= np.einsum( + 'ahpn,opn->ahop', + k.swapaxes(-2, -3).reshape(n_nops, n_ctrl, d2, d, order='F'), + i2.reshape(n_omega, d2, d, order='F') + ).reshape(n_nops, n_ctrl, n_omega, d, d, order='F').swapaxes(-1, -2) + + # Expand in basis transformed to eigenspace + ctrlmat_step_deriv = np.einsum('jnk,ahokn->ajho', 1j*basis_transformed, M) + + if n_coeffs_deriv is not None: + ctrlmat_step_deriv += ctrlmat_step_deriv_s + + return ctrlmat_step, ctrlmat_step_deriv + + +def _control_matrix_at_timestep_derivative_loop_2( + omega: Coefficients, + dt: Coefficients, + eigvals: ndarray, + eigvecs: ndarray, + basis_transformed, + c_opers_transformed: ndarray, + n_opers_transformed, + n_coeffs: Sequence[Coefficients], + n_coeffs_deriv: Sequence[Coefficients], + ctrlmat_step, + phase_factor, + deriv_integral, + exp_buf, + integral, + expr, + out +) -> Tuple[ndarray, ndarray]: + """".""" + d = len(eigvecs) + d2 = d**2 + n_ctrl = len(c_opers_transformed) + n_nops = len(n_opers_transformed) + n_omega = len(omega) + + deriv_integral = _derivative_integral(omega, eigvals, dt, out=deriv_integral) + + # Basically a tensor product + # K = np.einsum('taij,thkl->tahikjl', VBV, VHV) + # L = K.transpose(0, 1, 2, 4, 3, 6, 5) + K = util.tensor(n_opers_transformed[:, None], c_opers_transformed[None]) + L = util.tensor_transpose(K, (1, 0), [[d, d], [d, d]]) + k = np.diagonal(K.reshape(n_nops, n_ctrl, d, d, d, d), 0, -2, -3) + l = np.diagonal(L.reshape(n_nops, n_ctrl, d, d, d, d), 0, -2, -3) + i1 = np.diagonal(deriv_integral, 0, -2, -3) + i2 = np.diagonal(deriv_integral, 0, -1, -4) + # reshaping in F-major is significantly faster than C-major (~ factor 2-4) + M = np.einsum( + 'ahpm,opm->ahop', + l.reshape(n_nops, n_ctrl, d2, d, order='F'), + i1.reshape(n_omega, d2, d, order='F') + ).reshape(n_nops, n_ctrl, n_omega, d, d, order='F') + if d == 2: + # Life a bit simpler + mask = np.eye(d, dtype=bool) + M[..., mask] -= M[..., mask][..., ::-1] + M[..., ~mask] *= 2 + else: + M -= np.einsum( + 'ahpn,opn->ahop', + k.swapaxes(-2, -3).reshape(n_nops, n_ctrl, d2, d, order='F'), + i2.reshape(n_omega, d2, d, order='F') + ).reshape(n_nops, n_ctrl, n_omega, d, d, order='F').swapaxes(-1, -2) + + # Expand in basis transformed to eigenspace + ctrlmat_step_deriv = np.einsum('jnk,ahokn->ajho', 1j*basis_transformed, M) + ctrlmat_step_deriv *= phase_factor + + if n_coeffs_deriv is not None: + # equivalent contraction: 'ah,a,ako->akho', but this faster + ctrlmat_step_deriv += ((n_coeffs_deriv / n_coeffs[:, None])[:, None, :, None] + * ctrlmat_step[:, :, None]) + + return ctrlmat_step_deriv + + +def noise_operators_at_timestep_derivative( + omega: Coefficients, + dt: Coefficients, + eigvals: ndarray, + eigvecs: ndarray, + basis: Basis, + n_opers: Sequence[Operator], + n_coeffs: Sequence[Coefficients], + VHV: ndarray, + s_derivs: Optional[Sequence[Coefficients]] = None, + show_progressbar: Optional[bool] = None, + intermediates: Optional[Dict[str, ndarray]] = None +): + d = eigvecs.shape[-1] + n_dt = len(dt) + n_ctrl = len(VHV[0]) + n_nops = len(n_opers) + E = omega + + if intermediates is None: + noise_operators = np.empty((n_dt, len(E), n_nops, d, d), dtype=complex) + VBV = numeric._transform_hamiltonian(eigvecs, n_opers, n_coeffs).swapaxes(0, 1) + exp_buf, integral = np.empty((2, len(E), d, d), dtype=complex) + else: + noise_operators = intermediates['noise_operators_step'] + VBV = intermediates['n_opers_transformed'].swapaxes(0, 1) + + # R_g_d_s = np.empty((n_dt, len(E), n_nops, n_ctrl, d, d), dtype=complex) + # For calculating dR_g: d,d = p,q, d,d = m,n + integral_deriv = np.empty((n_dt, len(E), d, d, d, d), dtype=complex) + + expr_1 = oe.contract_expression('akl,okl->oakl', VBV[:, 0].shape, (len(E), d, d)) + expr_2 = oe.contract_expression('ij,...jk,lk', + eigvecs[0].shape, noise_operators[0].shape, eigvecs[0].shape, + optimize=[(0, 1), (0, 1)]) + for g in range(n_dt): + if intermediates is None: + integral = numeric._first_order_integral(E, eigvals[g], dt[g], exp_buf, integral) + noise_operators[g] = expr_1(VBV[g], integral, out=noise_operators[g]) + noise_operators[g] = expr_2(eigvecs[g], noise_operators[g], eigvecs[g].conj(), + out=noise_operators[g]) + + # Add the dependency of the susceptibilities + # s_derivs has shape (n_nops, n_ctrl, n_dt) + # VCV has shape (n_dt, d**2, d, d) + # VBV has shape (n_dt, n_nops, d, d) + # integral_Rg has shape (n_omega, d, d) + # R_g_d_s shall be of shape (n_dt, n_nops, d**2, n_ctrl, n_omega) + # if s_derivs is not None: + # Still need to think about n_coeffs here + # R_g_d_s[g] = s_derivs[:, None, :, g, None] * noise_operators[g] + + integral_deriv[g] = _derivative_integral(E, eigvals[g], dt[g], integral_deriv[g]) + + # (np.einsum('tcmn,tanp,topmn->tocamp', VHV, VBV, I1, optimize=['einsum_path', (0, 1), (0, 1)]) - + # np.einsum('tcnp,tamn,tompn->tocamp', VHV, VBV, I2, optimize=['einsum_path', (0, 1), (0, 1)])) + I1 = np.diagonal(integral_deriv, 0, -4, -1) + path = [(0, 1), (0, 1)] # optimization says [(0, 1), (0, 1)] + # M = oe.contract('tcmn,tanp,topmn->tocamp', VHV, VBV, I1, optimize=path) + M = oe.contract('tcmn,tanp,tompn->tocamp', VHV, VBV, I1, optimize=path) + if d == 2: + mask = np.eye(2, dtype=bool) + M[..., mask] -= M[..., mask][..., ::-1] + M[..., ~mask] *= 2 + else: + I2 = np.diagonal(integral_deriv, 0, -3, -2) + # M -= oe.contract('tcnp,tamn,tompn->tocamp', VHV, VBV, I2, optimize=path) + M -= oe.contract('tcnp,tamn,topmn->tocamp', VHV, VBV, I2, optimize=path) + + intermediates = intermediates or dict() + intermediates.setdefault('noise_operators_step', noise_operators) + intermediates.setdefault('n_opers_transformed', VBV.swapaxes(0, 1)) + + return 1j*M, intermediates + + +def calculate_derivative_of_control_matrix_from_scratch( + omega: Coefficients, + propagators: ndarray, + eigvals: ndarray, + eigvecs: ndarray, + basis: Basis, + t: Coefficients, + dt: Coefficients, + n_opers: Sequence[Operator], + n_coeffs: Sequence[Coefficients], + c_opers: Sequence[Operator], + all_identifiers: Sequence[str], + control_identifiers: Optional[Sequence[str]] = None, + s_derivs: Optional[Sequence[Coefficients]] = None) -> ndarray: + r""" + Calculate the derivative of the control matrix from scratch. + + Parameters + ---------- + omega : array_like, shape (n_omega,) + Frequencies, at which the pulse control matrix is to be evaluated. + propagators : array_like, shape (n_dt+1, d, d) + The propagators :math:`Q_g = P_g P_{g-1}\cdots P_0` as a (d, d) array + with *d* the dimension of the Hilbert space. + Q_Liou : ndarray, shape (n_dt+1, d**2, d**2) + The Liouville representation of the cumulative control propagators + U_c(t_g,0). + eigvals : array_like, shape (n_dt, d) + Eigenvalue vectors for each time pulse segment *g* with the first + axis counting the pulse segment, i.e. + ``HD == array([D_0, D_1, ...])``. + eigvecs : array_like, shape (n_dt, d, d) + Eigenvector matrices for each time pulse segment *g* with the first + axis counting the pulse segment, i.e. + ``HV == array([V_0, V_1, ...])``. + basis : Basis, shape (d**2, d, d) + The basis elements, in which the pulse control matrix will be expanded. + t : array_like, shape (n_dt+1), optional + The absolute times of the different segments. + dt : array_like, shape (n_dt) + Sequence duration, i.e. for the :math:`g`-th pulse + :math:`t_g - t_{g-1}`. + n_opers : array_like, shape (n_nops, d, d) + Noise operators :math:`B_\alpha`. + 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. + c_opers : array_like, shape (n_cops, d, d) + Control operators :math:`H_k`. + all_identifiers : array_like, shape (n_cops) + Identifiers of all control operators. + control_identifiers : Sequence[str], shape (n_ctrl), Optional + Sequence of strings with the control identifiers to distinguish between + accessible control and drift Hamiltonian. The default is None. + s_derivs : array_like, shape (n_nops, n_ctrl, n_dt) + The derivatives of the noise susceptibilities by the control amplitudes. + Defaults to None. + + Raises + ------ + ValueError + If the given identifiers *control_identifier* are not subset of the + total identifiers *all_identifiers* of all control operators. + + Returns + ------- + dR : array_like, shape (n_nops, d**2, n_dt, n_ctrl, n_omega) + Partial derivatives of the total control matrix with respect to each + control direction + :math:`\frac{\partial R_{\alpha k}(\omega)}{\partial u_h(t_{g'})}`. + The array's indexing has shape :math:`(\alpha,k,g^\prime,h,\omega)`. + + Notes + ----- + The derivative of the control matrix is calculated according to + + .. math :: + + \frac{\partial R_{\alpha k}(\omega)}{\partial u_h(t_{g'})} = + \sum_{g=1}^G \mathrm{e}^{\mathrm{i}\omega t_{g-1}}\cdot\left(\sum_j + \left[\frac{\partial R_{\alpha j}^{(g)}(\omega)} + {\partial u_h(t_{g'})} \cdot \mathcal{Q}_{jk}^{(g-1)} + + R_{\alpha j}^{(g)}(\omega) + \cdot\frac{\partial \mathcal{Q}_{jk}^{(g-1)}} + {\partial u_h(t_{g'})} \right] \right) + + See Also + -------- + :func:`liouville_derivative` + :func:`control_matrix_at_timestep_derivative` + """ + # Distinction between control and drift operators and only calculate the + # derivatives in control direction + + path = ['einsum_path', (0, 1), (0, 1)] + if (control_identifiers is None): + VHV = np.einsum('tji,hjk,tkl->thil', eigvecs.conj(), c_opers, eigvecs, + optimize=path) + elif (set(control_identifiers) <= set(all_identifiers)): + dict_id_oper = dict(zip(all_identifiers, c_opers)) + control = [dict_id_oper[element] for element in control_identifiers] + VHV = np.einsum('tji,hjk,tkl->thil', eigvecs.conj(), control, eigvecs, + optimize=path) + else: + raise ValueError('Given control identifiers have to be a \ + subset of (drift+control) Hamiltonian!') + + Q_Liou = superoperator.liouville_representation(propagators, basis) + # Get derivative of Liouville, control matrix at each time step, derivative + # derivative of control matrix at each time step + dL = liouville_derivative( + dt=dt, propagators=propagators, basis=basis, eigvecs=eigvecs, + eigvals=eigvals, transf_control_operators=VHV) + ctrlmat_data = control_matrix_at_timestep_derivative( + omega=omega, + dt=dt, + eigvals=eigvals, + eigvecs=eigvecs, + basis=basis, + n_opers=n_opers, + n_coeffs=n_coeffs, + transf_control_operators=VHV, + s_derivs=s_derivs + ) + ctrlmat_g = ctrlmat_data['R_g'] + ctrlmat_g_deriv = ctrlmat_data['dR_g'] + + # Calculate the derivative of the total control matrix + exp_factor = util.cexp(np.multiply.outer(t, omega)) + summand1 = np.einsum('to,tajho,tjk->aktho', exp_factor[:-1], + ctrlmat_g_deriv, Q_Liou[:-1], + optimize=['einsum_path', (1, 2), (0, 1)]) + summand2 = np.einsum('to,tajo,thsjk->aksho', exp_factor[1:-1], + ctrlmat_g[1:], dL[:-1], + optimize=['einsum_path', (0, 1), (0, 1)]) + + dR = summand1 + summand2 + return dR + + +def calculate_derivative_of_control_matrix_from_scratch_refactored( + omega: Coefficients, + propagators: ndarray, + eigvals: ndarray, + eigvecs: ndarray, + basis: Basis, + t: Coefficients, + dt: Coefficients, + n_opers: Sequence[Operator], + n_coeffs: Sequence[Coefficients], + c_opers: Sequence[Operator], + all_identifiers: Sequence[str], + control_identifiers: Optional[Sequence[str]] = None, + n_coeffs_deriv: Optional[Sequence[Coefficients]] = None) -> ndarray: + r""".""" + # Distinction between control and drift operators and only calculate the + # derivatives in control direction + + if control_identifiers is None: + control = c_opers + elif set(control_identifiers) <= set(all_identifiers): + dict_id_oper = dict(zip(all_identifiers, c_opers)) + control = [dict_id_oper[element] for element in control_identifiers] + else: + raise ValueError('Given control identifiers have to be a \ + subset of (drift+control) Hamiltonian!') + + c_opers_transformed = numeric._transform_hamiltonian(eigvecs, control).swapaxes(0, 1) + propagators_liouville = superoperator.liouville_representation(propagators[:-1], basis) + + propagators_liouville_deriv = liouville_derivative_vectorized(dt, propagators, basis, eigvecs, + eigvals, c_opers_transformed) + # ctrlmat_g, ctrlmat_g_deriv = control_matrix_at_timestep_derivative_refactored( + ctrlmat_g, ctrlmat_g_deriv = control_matrix_at_timestep_derivative_faster( + omega, dt, eigvals, eigvecs, basis, c_opers_transformed, n_opers, n_coeffs, n_coeffs_deriv + ) + + # Calculate the derivative of the total control matrix + exp_factor = util.cexp(np.multiply.outer(omega, t[:-1])) + # Equivalent (but slower) einsum contraction for first term: + # ctrlmat_deriv = np.einsum('ot,tajho,tjk->hotak', + # exp_factor, ctrlmat_g_deriv, propagators_liouville, + # optimize=['einsum_path', (0, 2), (0, 1)]) + ctrlmat_deriv = (ctrlmat_g_deriv.transpose(3, 4, 0, 1, 2) + @ (exp_factor[..., None, None] * propagators_liouville)) + ctrlmat_deriv += np.einsum('ot,tajo,thsjk->hosak', + exp_factor[:, 1:], ctrlmat_g[1:], propagators_liouville_deriv, + optimize=['einsum_path', (0, 1), (0, 1)]) + + return ctrlmat_deriv + + +def calculate_derivative_of_control_matrix_from_scratch_loop_1( + omega: Coefficients, + propagators: ndarray, + eigvals: ndarray, + eigvecs: ndarray, + basis: Basis, + t: Coefficients, + dt: Coefficients, + n_opers: Sequence[Operator], + n_coeffs: Sequence[Coefficients], + c_opers: Sequence[Operator], + intermediates, + all_identifiers: Sequence[str], + control_identifiers: Optional[Sequence[str]] = None, + n_coeffs_deriv: Optional[Sequence[Coefficients]] = None +) -> ndarray: + r""".""" + # Distinction between control and drift operators and only calculate the + # derivatives in control direction + if control_identifiers is None: + control = c_opers + elif set(control_identifiers) <= set(all_identifiers): + dict_id_oper = dict(zip(all_identifiers, c_opers)) + control = [dict_id_oper[element] for element in control_identifiers] + else: + raise ValueError('Given control identifiers have to be a \ + subset of (drift+control) Hamiltonian!') + + d = eigvecs.shape[-1] + d2 = d**2 + n_dt = len(dt) + n_ctrl = len(control) + n_nops = len(n_opers) + n_omega = len(omega) + + # Precompute some transformation into eigenspace of control Hamiltonian + c_opers_transformed = numeric._transform_hamiltonian(eigvecs, control).swapaxes(0, 1) + n_opers_transformed = numeric._transform_hamiltonian(eigvecs, n_opers, n_coeffs).swapaxes(0, 1) + basis_transformed = numeric._transform_by_unitary(eigvecs[:, None], basis[None], + out=np.empty((n_dt, d2, d, d), complex)) + propagators_liouville = superoperator.liouville_representation(propagators[:-1], basis) + + propagators_liouville_deriv = liouville_derivative_vectorized(dt, propagators, basis, eigvecs, + eigvals, c_opers_transformed) + exp_factor = util.cexp(np.multiply.outer(omega, t[:-1])) + + # Allocate memory + ctrlmat_deriv = np.empty((n_ctrl, n_omega, n_dt, n_nops, d2), dtype=complex) + ctrlmat_step = np.empty((n_dt, n_nops, d2, n_omega), dtype=complex) + deriv_integral = np.empty((n_omega, d, d, d, d), dtype=complex) + exp_buf, integral = np.empty((2, n_omega, d, d), dtype=complex) + ctrlmat_expr = oe.contract_expression('bcd,adc,hdc->abh', basis.shape, basis.shape, + integral.shape, optimize=[(0, 1), (0, 1)]) + for g in range(n_dt): + ctrlmat_step[g], ctrlmat_step_deriv = _control_matrix_at_timestep_derivative_loop_1( + omega, dt[g], eigvals[g], eigvecs[g], basis_transformed[g], c_opers_transformed[g], + n_opers_transformed[g], n_coeffs[:, g], n_coeffs_deriv[:, :, g], ctrlmat_step[g], + deriv_integral, exp_buf, integral, ctrlmat_expr + ) + + # Calculate the derivative of the total control matrix + ctrlmat_deriv[:, :, g] = ((exp_factor[:, g] * ctrlmat_step_deriv).transpose(2, 3, 0, 1) + @ propagators_liouville[g]) + + ctrlmat_deriv += oe.contract('ot,tajo,thsjk->hosak', + exp_factor[:, 1:], ctrlmat_step[1:], propagators_liouville_deriv, + optimize=[(0, 1), (0, 1)]) + + return ctrlmat_deriv + + +def calculate_derivative_of_control_matrix_from_scratch_loop_2( + omega: Coefficients, + propagators: ndarray, + eigvals: ndarray, + eigvecs: ndarray, + basis: Basis, + t: Coefficients, + dt: Coefficients, + n_opers: Sequence[Operator], + n_coeffs: Sequence[Coefficients], + c_opers: Sequence[Operator], + intermediates, + all_identifiers: Sequence[str], + control_identifiers: Optional[Sequence[str]] = None, + n_coeffs_deriv: Optional[Sequence[Coefficients]] = None +) -> ndarray: + r""".""" + # Distinction between control and drift operators and only calculate the + # derivatives in control direction + if control_identifiers is None: + control = c_opers + elif set(control_identifiers) <= set(all_identifiers): + dict_id_oper = dict(zip(all_identifiers, c_opers)) + control = [dict_id_oper[element] for element in control_identifiers] + else: + raise ValueError('Given control identifiers have to be a \ + subset of (drift+control) Hamiltonian!') + + d = eigvecs.shape[-1] + d2 = d**2 + n_dt = len(dt) + n_ctrl = len(control) + n_nops = len(n_opers) + n_omega = len(omega) + + # Precompute some transformation into eigenspace of control Hamiltonian + c_opers_transformed = numeric._transform_hamiltonian(eigvecs, control).swapaxes(0, 1) + n_opers_transformed = numeric._transform_hamiltonian(eigvecs, n_opers, n_coeffs).swapaxes(0, 1) + basis_transformed = numeric._transform_by_unitary(eigvecs[:, None], basis[None], + out=np.empty((n_dt, d2, d, d), complex)) + # U_deriv = liouville_derivative_vectorized_loop(dt, propagators, basis, eigvecs, eigvals, + # c_opers_transformed) + propagators_liouville = superoperator.liouville_representation(propagators[:-1], basis) + + # propagators_liouville_deriv = np.empty((n_ctrl, n_dt, d2, d2), dtype=complex) + # propagators_liouville_deriv = liouville_derivative_vectorized(dt, propagators, basis, eigvecs, + # eigvals, c_opers_transformed) + propagators_deriv = propagators_derivative_vectorized_completely(dt, propagators, basis, + eigvecs, eigvals, + c_opers_transformed) + liouville_deriv = np.empty((n_ctrl, n_dt, d2, d2), dtype=complex) + # Need to 'remove' the propagators from the control matrix at time step as computed by + # numeric.calculate_control_matrix_from_scratch + ctrlmat_step = (intermediates['control_matrix_step'].transpose(3, 0, 1, 2) + @ propagators_liouville.conj().swapaxes(-1, -2)).transpose(1, 2, 3, 0) + # exp_factor = util.cexp(np.multiply.outer(omega, t[:-1])) + + # Allocate memory + ctrlmat_deriv = np.empty((n_ctrl, n_omega, n_dt, n_nops, d2), dtype=complex) + deriv_integral = np.empty((n_omega, d, d, d, d), dtype=complex) + ctrlmat_step_deriv = np.empty((n_nops, d2, n_ctrl, n_omega), dtype=complex) + exp_buf, integral = np.empty((2, n_omega, d, d), dtype=complex) + ctrlmat_expr = oe.contract_expression('bcd,adc,hdc->abh', basis.shape, basis.shape, + integral.shape, optimize=[(0, 1), (0, 1)]) + liouville_deriv_expr = oe.contract_expression('hsba,jbc,cd,kda->hsjk', + propagators_deriv.conj()[:, 0].shape, + basis.shape, propagators[0].shape, basis.shape, + optimize=[(0, 2), (0, 2), (0, 1)]) + for g in range(n_dt): + phase_factor = util.cexp(omega * t[g]) + # Calculate the whole propagator derivative + # propagators_deriv[:, :g+1] = (propagators[g+1] + # @ propagators[1:g+2].conj().swapaxes(-1, -2) + # @ U_deriv[:, :g+1] + # @ propagators[:g+1]) + + # # Now calculate derivative of Liouville representation. Operation is faster using matmul, but + # # not nice to read. Full einsum contraction would be as follows (we insert new axes to induce + # # an outer product between basis elements steps using broadcasting): + # # liouville_deriv = np.einsum('hsba,jbc,cd,kda->hsjk', + # # deriv_U.conj(), basis, propagators[g+1], basis) + # propagators_liouville_deriv = np.einsum('hsba,jkba->hsjk', propagators_deriv.conj(), + # (basis @ propagators[g+1])[:, None] @ basis, + # out=propagators_liouville_deriv) + + ctrlmat_step_deriv = _control_matrix_at_timestep_derivative_loop_2( + omega, dt[g], eigvals[g], eigvecs[g], basis_transformed[g], c_opers_transformed[g], + n_opers_transformed[g], n_coeffs[:, g], n_coeffs_deriv[:, :, g], ctrlmat_step[g], + phase_factor, deriv_integral, exp_buf, integral, ctrlmat_expr, out=ctrlmat_step_deriv + ) + + # Calculate the derivative of the total control matrix (phase already included) + ctrlmat_deriv[:, :, g] = (ctrlmat_step_deriv.transpose(2, 3, 0, 1) + @ propagators_liouville[g]) + + if g < n_dt - 1: + liouville_deriv = liouville_derivative_vectorized_completely_loop( + propagators_deriv[:, g], basis, propagators[g+1], liouville_deriv_expr, + out=liouville_deriv + ) + ctrlmat_deriv += oe.contract('ajo,hsjk->hosak', + ctrlmat_step[g+1], liouville_deriv.real) + + # ctrlmat_deriv += oe.contract('ot,tajo,thsjk->hosak', + # exp_factor[:, 1:], ctrlmat_g[1:], propagators_liouville_deriv, + # optimize=[(0, 1), (0, 1)]) + # ctrlmat_deriv += oe.contract('otaj,thsjk->hosak', prefactor, propagators_liouville_deriv) + + return ctrlmat_deriv + + +def calculate_derivative_of_control_matrix_from_scratch_loop_3( + omega: Coefficients, + propagators: ndarray, + eigvals: ndarray, + eigvecs: ndarray, + basis: Basis, + t: Coefficients, + dt: Coefficients, + n_opers: Sequence[Operator], + n_coeffs: Sequence[Coefficients], + c_opers: Sequence[Operator], + intermediates, + all_identifiers: Sequence[str], + control_identifiers: Optional[Sequence[str]] = None, + n_coeffs_deriv: Optional[Sequence[Coefficients]] = None +) -> ndarray: + r""".""" + # Distinction between control and drift operators and only calculate the + # derivatives in control direction + if control_identifiers is None: + control = c_opers + elif set(control_identifiers) <= set(all_identifiers): + dict_id_oper = dict(zip(all_identifiers, c_opers)) + control = [dict_id_oper[element] for element in control_identifiers] + else: + raise ValueError('Given control identifiers have to be a \ + subset of (drift+control) Hamiltonian!') + + d = eigvecs.shape[-1] + d2 = d**2 + n_dt = len(dt) + n_ctrl = len(control) + n_nops = len(n_opers) + n_omega = len(omega) + + # Precompute some transformation into eigenspace of control Hamiltonian + c_opers_transformed = numeric._transform_hamiltonian(eigvecs, control).swapaxes(0, 1) + n_opers_transformed = numeric._transform_hamiltonian(eigvecs, n_opers, n_coeffs).swapaxes(0, 1) + basis_transformed = numeric._transform_by_unitary(eigvecs[:, None], basis[None], + out=np.empty((n_dt, d2, d, d), complex)) + # U_deriv = liouville_derivative_vectorized_loop(dt, propagators, basis, eigvecs, eigvals, + # c_opers_transformed) + propagators_liouville = superoperator.liouville_representation(propagators[:-1], basis) + + # propagators_liouville_deriv = np.empty((n_ctrl, n_dt, d2, d2), dtype=complex) + # propagators_liouville_deriv = liouville_derivative_vectorized(dt, propagators, basis, eigvecs, + # eigvals, c_opers_transformed) + propagators_deriv = propagators_derivative_vectorized_completely(dt, propagators, basis, + eigvecs, eigvals, + c_opers_transformed) + # Contraction 'jbc,tcd,kda->tjbka' + basis_propagators = np.tensordot(np.tensordot(propagators[1:-1], basis, axes=[1, 2]), + basis, axes=[1, 1]) + propagators_liouville_deriv = oe.contract('tjnkm,htsnm->tsjkh', + basis_propagators, propagators_deriv) + + # Need to 'remove' the propagators from the control matrix at time step as computed by + # numeric.calculate_control_matrix_from_scratch + ctrlmat_step = (intermediates['control_matrix_step'].transpose(3, 0, 1, 2) + @ propagators_liouville.conj().swapaxes(-1, -2)).transpose(1, 2, 3, 0) + # exp_factor = util.cexp(np.multiply.outer(omega, t[:-1])) + + # Allocate memory + ctrlmat_deriv = np.empty((n_ctrl, n_omega, n_dt, n_nops, d2), dtype=complex) + deriv_integral = np.empty((n_omega, d, d, d, d), dtype=complex) + ctrlmat_step_deriv = np.empty((n_nops, d2, n_ctrl, n_omega), dtype=complex) + exp_buf, integral = np.empty((2, n_omega, d, d), dtype=complex) + ctrlmat_expr = oe.contract_expression('bcd,adc,hdc->abh', basis.shape, basis.shape, + integral.shape, optimize=[(0, 1), (0, 1)]) + liouville_deriv_expr = oe.contract_expression('tajo,htnm,jnp,tpq,kqm->hoak', + ctrlmat_step[1:].shape, + propagators_deriv.conj()[:, :, 0].shape, + basis.shape, propagators[1:-1].shape, + basis.shape, + optimize=[(2, 3), (2, 3), (1, 2), (0, 1)]) + for g in range(n_dt): + phase_factor = util.cexp(omega * t[g]) + # Calculate the whole propagator derivative + # propagators_deriv[:, :g+1] = (propagators[g+1] + # @ propagators[1:g+2].conj().swapaxes(-1, -2) + # @ U_deriv[:, :g+1] + # @ propagators[:g+1]) + + # # Now calculate derivative of Liouville representation. Operation is faster using matmul, but + # # not nice to read. Full einsum contraction would be as follows (we insert new axes to induce + # # an outer product between basis elements steps using broadcasting): + # # liouville_deriv = np.einsum('hsba,jbc,cd,kda->hsjk', + # # deriv_U.conj(), basis, propagators[g+1], basis) + # propagators_liouville_deriv = np.einsum('hsba,jkba->hsjk', propagators_deriv.conj(), + # (basis @ propagators[g+1])[:, None] @ basis, + # out=propagators_liouville_deriv) + + ctrlmat_step_deriv = _control_matrix_at_timestep_derivative_loop_2( + omega, dt[g], eigvals[g], eigvecs[g], basis_transformed[g], c_opers_transformed[g], + n_opers_transformed[g], n_coeffs[:, g], n_coeffs_deriv[:, :, g], ctrlmat_step[g], + phase_factor, deriv_integral, exp_buf, integral, ctrlmat_expr, out=ctrlmat_step_deriv + ) + + # Calculate the derivative of the total control matrix (phase already included) + ctrlmat_deriv[:, :, g] = (ctrlmat_step_deriv.transpose(2, 3, 0, 1) + @ propagators_liouville[g]) + + propagators_liou_deriv = oe.contract('htba,jbc,tcd,kda->thjk', propagators_deriv[:, :, g], + basis, propagators[1:-1], basis) + ctrlmat_deriv[:, :, g] += oe.contract('tajo,thjk->hoak', ctrlmat_step[1:], + propagators_liou_deriv) + + # ctrlmat_deriv[:, :, g] += oe.contract('tajo,htnm,tjnkm->hoak', ctrlmat_step[1:], + # propagators_deriv.conj()[:, :, g], basis_propagators, + # optimize=[(1, 2), (0, 1)]) + + ctrlmat_deriv[:, :, g] += oe.contract('tajo,tjkh->hoak', + ctrlmat_step[1:], propagators_liouville_deriv[:, g]) + + # ctrlmat_deriv[:, :, g] += liouville_deriv_expr(ctrlmat_step[1:], + # propagators_deriv.conj()[:, :, 0], + # basis, propagators[1:-1], basis) + + # ctrlmat_deriv += oe.contract('ot,tajo,thsjk->hosak', + # exp_factor[:, 1:], ctrlmat_step[1:], propagators_liouville_deriv, + # optimize=[(0, 1), (0, 1)]) + # ctrlmat_deriv += oe.contract('otaj,thsjk->hosak', prefactor, propagators_liouville_deriv) + + return ctrlmat_deriv + + +def calculate_derivative_of_control_matrix_from_scratch_loop_4( + omega: Coefficients, + propagators: ndarray, + eigvals: ndarray, + eigvecs: ndarray, + basis: Basis, + t: Coefficients, + dt: Coefficients, + n_opers: Sequence[Operator], + n_coeffs: Sequence[Coefficients], + c_opers: Sequence[Operator], + intermediates, + all_identifiers: Sequence[str], + control_identifiers: Optional[Sequence[str]] = None, + n_coeffs_deriv: Optional[Sequence[Coefficients]] = None +) -> ndarray: + r""".""" + # Distinction between control and drift operators and only calculate the + # derivatives in control direction + if control_identifiers is None: + control = c_opers + elif set(control_identifiers) <= set(all_identifiers): + dict_id_oper = dict(zip(all_identifiers, c_opers)) + control = [dict_id_oper[element] for element in control_identifiers] + else: + raise ValueError('Given control identifiers have to be a \ + subset of (drift+control) Hamiltonian!') + + d = eigvecs.shape[-1] + d2 = d**2 + n_dt = len(dt) + n_ctrl = len(control) + n_nops = len(n_opers) + n_omega = len(omega) + + # Precompute some transformation into eigenspace of control Hamiltonian + c_opers_transformed = numeric._transform_hamiltonian(eigvecs, control).swapaxes(0, 1) + n_opers_transformed = numeric._transform_hamiltonian(eigvecs, n_opers, n_coeffs).swapaxes(0, 1) + basis_transformed = numeric._transform_by_unitary(eigvecs[:, None], basis[None], + out=np.empty((n_dt, d2, d, d), complex)) + propagators_liouville = superoperator.liouville_representation(propagators[:-1], basis) + + propagators_liouville_deriv = liouville_derivative_vectorized(dt, propagators, basis, eigvecs, + eigvals, c_opers_transformed) + ctrlmat_step = intermediates['control_matrix_step'] + prefactor = (ctrlmat_step[1:].transpose(-1, 0, 1, 2) + @ propagators_liouville[1:].conj().swapaxes(-1, -2)) + + # Allocate memory + ctrlmat_step_deriv = np.empty((n_nops, d2, n_ctrl, n_omega), dtype=complex) + ctrlmat_deriv = np.empty((n_ctrl, n_omega, n_dt, n_nops, d2), dtype=complex) + deriv_integral = np.empty((n_omega, d, d, d, d), dtype=complex) + exp_buf, integral = np.empty((2, n_omega, d, d), dtype=complex) + ctrlmat_expr = oe.contract_expression('bcd,adc,hdc->abh', basis.shape, basis.shape, + integral.shape, optimize=[(0, 1), (0, 1)]) + for g in range(n_dt): + ctrlmat_step_deriv = _control_matrix_at_timestep_derivative_loop_2( + omega, dt[g], eigvals[g], eigvecs[g], basis_transformed[g], c_opers_transformed[g], + n_opers_transformed[g], n_coeffs[:, g], n_coeffs_deriv[:, :, g], ctrlmat_step[g], + deriv_integral, exp_buf, integral, ctrlmat_expr, out=ctrlmat_step_deriv + ) + + # Calculate the derivative of the total control matrix + exp_factor = util.cexp(omega * t[g]) + ctrlmat_deriv[:, :, g] = ((exp_factor * ctrlmat_step_deriv).transpose(2, 3, 0, 1) + @ propagators_liouville[g]) + + ctrlmat_deriv += oe.contract('otaj,thsjk->hosak', prefactor, propagators_liouville_deriv) + + return ctrlmat_deriv + + +def calculate_derivative_of_noise_operators_from_scratch( + omega: Coefficients, + propagators: ndarray, + eigvals: ndarray, + eigvecs: ndarray, + basis: Basis, + t: Coefficients, + dt: Coefficients, + n_opers: Sequence[Operator], + n_coeffs: Sequence[Coefficients], + c_opers: Sequence[Operator], + all_identifiers: Sequence[str], + control_identifiers: Optional[Sequence[str]] = None, + s_derivs: Optional[Sequence[Coefficients]] = None, + intermediates: Optional[Dict[str, ndarray]] = None): + + if control_identifiers is None: + control = c_opers + elif set(control_identifiers) <= set(all_identifiers): + dict_id_oper = dict(zip(all_identifiers, c_opers)) + control = [dict_id_oper[element] for element in control_identifiers] + else: + raise ValueError('Given control identifiers have to be a \ + subset of (drift+control) Hamiltonian!') + + d = eigvecs.shape[-1] + VHV = numeric._transform_hamiltonian(eigvecs, control).swapaxes(0, 1) + VQ = numeric._propagate_eigenvectors(eigvecs, propagators[:-1]) + + exp_factor = util.cexp(np.multiply.outer(t[:-1], omega)) + M, intermediates = noise_operators_at_timestep_derivative(omega, dt, eigvals, eigvecs, basis, + n_opers, n_coeffs, VHV, + intermediates=intermediates) + + QVMVQ = numeric._transform_by_unitary(M, VQ[:, None, None, None], out=M) + QVMVQ *= exp_factor[..., None, None, None, None] + + noise_operators_step_pulse = intermediates['noise_operators_step_pulse'] + noise_operators_step_total = intermediates['noise_operators_step_total'] + + propagators_deriv = propagator_derivative(dt, propagators, eigvecs, eigvals, VHV) + propagators_deriv_factor = propagator_derivative_factor(dt, propagators, eigvecs, eigvals, VHV) + + noise_operators_deriv = np.zeros((len(omega), len(n_opers), d, len(c_opers), len(dt), d), + dtype=complex) + for g in range(len(dt)): + noise_operators_deriv[..., :g+1, :] += np.tensordot(noise_operators_step_total[g], + propagators_deriv_factor[:, g, :g+1], + axes=[-1, -2]) + noise_operators_deriv[..., :g+1, :] += np.tensordot(noise_operators_step_total[g], + propagators_deriv_factor[:, g, :g+1].conj(), + axes=[-2, -2]) + + noise_operators_deriv = np.zeros(QVMVQ.shape, dtype=complex) + noise_operators_deriv = np.einsum( + 'so,sji,soajk,hstkl->tohail', + exp_factor[:-1], propagators[1:-1].conj(), noise_operators_step_pulse[1:], propagators_deriv[:, :-1], + optimize=['einsum_path', (0, 1), (0, 2), (0, 1)] + ) + noise_operators_deriv += np.einsum( + 'so,hstji,soajk,skl->tohail', + exp_factor[:-1], propagators_deriv[:, :-1].conj(), noise_operators_step_pulse[1:], propagators[1:-1], + optimize=['einsum_path', (0, 1), (0, 2), (0, 1)] + ) + + noise_operators_deriv = np.zeros(QVMVQ.shape, dtype=complex) + noise_operators_deriv.real = 2 * np.einsum( + 'to,hstjk,skl,soali->tohaji', + exp_factor, propagators_deriv[:, :-1], propagators[1:-1], noise_operators_step_pulse[1:], + optimize=['einsum_path', (1, 2), (1, 2), (0, 1)] + ).real + noise_operators_deriv += QVMVQ + + return noise_operators_deriv + + +def calculate_filter_function_derivative(R: ndarray, deriv_R: ndarray) -> ndarray: + r""" + Compute the filter function derivative from the control matrix. + + Parameters + ---------- + R : array_like, shape (n_nops, d**2, n_omega) + The control matrix. + deriv_R: array_like, shape (n_nops, d**2, n_t, n_ctrl, n_omega) + The derivative of the control matrix. + + Returns + ------- + deriv_filter_function : ndarray, shape (n_nops, n_dt, n_ctrl, n_omega) + The regular filter functions' derivatives for variation in each control + contribution + :math:`\frac{\partial F_\alpha(\omega)}{\partial u_h(t_{g'})}`. + The array's indexing has shape :math:`(\alpha,g^\prime,h,\omega)`. + + Notes + ----- + The filter function derivative is calculated according to + + .. math :: + + \frac{\partial F_\alpha(\omega)}{\partial u_h(t_{g'})} + = 2 \mathrm{Re} \left( \sum_k R_{\alpha k}^\ast(\omega) + \frac{\partial R_{\alpha k}(\omega)} + {\partial u_h(t_{g'})} \right) + """ + summe = np.einsum('ako,aktho->atho', R.conj(), deriv_R) + return 2*summe.real + + +def calculate_filter_function_derivative_refactored(R: ndarray, R_deriv: ndarray) -> ndarray: + summe = np.einsum('ako,hotak->atho', R.conj(), R_deriv) + return 2*summe.real + + +def infidelity_derivative( + pulse: 'PulseSequence', + S: Union[Coefficients, Callable], + omega: Coefficients, + control_identifiers: Optional[Sequence[str]] = None, + s_derivs: Optional[Sequence[Coefficients]] = None +) -> ndarray: + r""" + Calculate the infidelity derivative. + + Calculate the entanglement infidelity derivative of the ``PulseSequence`` + *pulse*. + + Parameters + ---------- + pulse : PulseSequence + The ``PulseSequence`` instance, for which to calculate the infidelity + for. + S : array_like + The two-sided noise power spectral density in units of inverse + frequencies as an array of shape (n_omega,) or (n_nops, n_omega). In + the first case, the same spectrum is taken for all noise operators, in + the second, it is assumed that there are no correlations between + different noise sources and thus there is one spectrum for each noise + operator. + omega : array_like, shape (n_omega) + The frequencies at which the integration is to be carried out. + control_identifiers : Sequence[str], shape (n_ctrl) + Sequence of strings with the control identifiern to distinguish between + accessible control and drift Hamiltonian. + s_derivs : array_like, shape (n_nops, n_ctrl, n_dt) + The derivatives of the noise susceptibilities by the control amplitudes. + Defaults to None. + + Raises + ------ + ValueError + If the provided noise spectral density does not fit expected shape. + + Returns + ------- + deriv_infid : ndarray, shape (n_nops, n_dt, n_ctrl) + Array with the derivative of the infidelity for each noise source taken + for each control direction at each time step + :math:`\frac{\partial I_e}{\partial u_h(t_{g'})}`. The array's indexing + has shape :math:`(\alpha,g^\prime,h)`. + + Notes + ----- + The infidelity's derivative is given by + + .. math:: + + \frac{\partial I_e}{\partial u_h(t_{g'})} = \frac{1}{d} + \int_{-\infty}^\infty + \frac{d\omega}{2\pi} + S_\alpha(\omega) + \frac{\partial F_\alpha (\omega)} + {\partial u_h(t_{g'})} + + with :math:`S_{\alpha}(\omega)` the noise spectral density + and :math:`F_{\alpha}(\omega)` the canonical filter function for + noise source :math:`\alpha`. + + To convert to the average gate infidelity, use the + following relation given by Horodecki et al. [Hor99]_ and + Nielsen [Nie02]_: + + .. math:: + + \big\langle\mathcal{I}_\mathrm{avg}\big\rangle = \frac{d}{d+1} + \big\langle\mathcal{I}_\mathrm{e}\big\rangle. + + References + ---------- + .. [Hor99] + Horodecki, M., Horodecki, P., & Horodecki, R. (1999). General + teleportation channel, singlet fraction, and quasidistillation. + Physical Review A - Atomic, Molecular, and Optical Physics, 60(3), + 1888–1898. https://doi.org/10.1103/PhysRevA.60.1888 + + .. [Nie02] + Nielsen, M. A. (2002). A simple formula for the average gate + fidelity of a quantum dynamical operation. Physics Letters, + Section A: General, Atomic and Solid State Physics, 303(4), 249–252. + https://doi.org/10.1016/S0375-9601(02)01272-0 + """ + n_ops = len(pulse.n_opers) + S = np.asarray(S) + if(S.shape[0] == n_ops): + S_all = S + elif(S.shape[0] == 1): + S_all = np.array([S[0]]*n_ops) + elif(S.shape[0] == len(omega)): + S_all = np.array([S]*n_ops) + else: + raise ValueError('Not fitting shape of S.') + + filter_function_deriv = pulse.get_filter_function_derivative(omega, control_identifiers, + s_derivs) + + integrand = np.einsum('ao,atho->atho', S_all, filter_function_deriv) + infid_deriv = trapz(integrand, omega) / (2*np.pi*pulse.d) + + return infid_deriv diff --git a/filter_functions/pulse_sequence.py b/filter_functions/pulse_sequence.py index 8c4a832..2c4696f 100644 --- a/filter_functions/pulse_sequence.py +++ b/filter_functions/pulse_sequence.py @@ -793,52 +793,39 @@ def get_pulse_correlation_filter_function(self, which: str = 'fidelity') -> ndar def get_filter_function_derivative( self, omega: Coefficients, - contorl_identifier: Optional[Sequence[str]] = None, - s_derivs: Optional[Sequence[Coefficients]] = None + control_identifiers: Optional[Sequence[str]] = None, + n_coeffs_deriv: Optional[Sequence[Coefficients]] = None ) -> ndarray: - r""" - Calculate the filter function derivative of the pulse sequence. - - The filter function derivative is cached so it doesn't need to be - calculated twice for the same frequencies. + r"""Calculate the pulse sequence's filter function derivative. Parameters ---------- - omega : array_like, shape (n_omega,) - Frequencies at which the pulse control matrix is to be evaluated. - contorl_identifier : Sequence[str] - Sequence of strings with the control identifiern to distinguish - between control and drift Hamiltonian. The default is None. - s_derivs : array_like, shape (n_nops, n_ctrl, n_dt) + omega: array_like, shape (n_omega,) + Frequencies at which the pulse control matrix is to be + evaluated. + control_identifiers: Sequence[str] + Sequence of strings with the control identifiern to + distinguish between control and drift Hamiltonian. The + default is None. + n_coeffs_deriv: array_like, shape (n_nops, n_ctrl, n_dt) The derivatives of the noise susceptibilities by the control amplitudes. Defaults to None. Returns ------- - deriv_filter_function : ndarray, shape (n_nops, n_t, n_ctrl, n_omega) - The regular filter functions' derivatives for variation in each - control contribution. + filter_function_deriv: ndarray, shape (n_nops, n_t, n_ctrl, n_omega) + The regular filter functions' derivatives for variation in + each control contribution. """ + control_matrix = self.get_control_matrix(omega, cache_intermediates=True) - R = self.get_control_matrix(omega) - eigvals, eigvecs, propagators, all_id = ( - self.eigvals, self.eigvecs, self._propagators, - self.c_oper_identifiers) - basis, dt, t = self.basis, self.dt, self.t - n_coeffs, c_coeffs = self.n_coeffs, self.c_coeffs - c_opers, n_opers = self.c_opers, self.n_opers - - Q_Liou = liouville_representation(propagators, basis) - deriv_R = gradient.calculate_derivative_of_control_matrix_from_scratch( - omega=omega, propagators=propagators, Q_Liou=Q_Liou, eigvals=eigvals, eigvecs=eigvecs, - basis=basis, t=t, dt=dt, n_opers=n_opers, n_coeffs=n_coeffs, - c_opers=c_opers, all_identifiers=all_id, control_identifiers=contorl_identifier, - s_derivs=s_derivs) - - return gradient.calculate_canonical_filter_function_derivative( - R, deriv_R + control_matrix_deriv = gradient.calculate_derivative_of_control_matrix_from_scratch( + omega, self.propagators, self.eigvals, self.eigvecs, self.basis, self.t, self.dt, + self.n_opers, self.n_coeffs, self.c_opers, self.c_oper_identifiers, + control_identifiers, n_coeffs_deriv, self._intermediates ) + return gradient.calculate_filter_function_derivative(control_matrix, control_matrix_deriv) def get_total_phases(self, omega: Coefficients) -> ndarray: """Get the (cached) total phase factors for this pulse and omega.""" @@ -1024,8 +1011,10 @@ def cleanup(self, method: str = 'conservative') -> None: elif method == 'frequency dependent': attrs = filter_function_attrs.union({'_control_matrix', '_control_matrix_pc', '_total_phases'}) - # Remove frequency dependent control_matrix_step from intermediates + # Remove frequency dependent terms from intermediates self._intermediates.pop('control_matrix_step', None) + self._intermediates.pop('phase_factors', None) + self._intermediates.pop('first_order_integral', None) else: # method == all attrs = filter_function_attrs.union(default_attrs, concatenation_attrs) diff --git a/tests/gradient_testutil.py b/tests/gradient_testutil.py index 755d51c..e2c11bd 100644 --- a/tests/gradient_testutil.py +++ b/tests/gradient_testutil.py @@ -26,10 +26,12 @@ def deriv_2_exchange_interaction(eps): def one_over_f_noise(f): - return S_0 / f + spectrum = np.divide(S_0, f, where=(f != 0)) + spectrum[f == 0] = spectrum[np.abs(f).argmin()] + return spectrum -def create_sing_trip_pulse_seq(eps, dbz): +def create_sing_trip_pulse_seq(eps, dbz, *args): H_c = [ [sigma_z, exchange_interaction(eps[0]), 'control1'], @@ -48,7 +50,21 @@ def create_sing_trip_pulse_seq(eps, dbz): return pulse -def finite_diff_infid(u_ctrl_central, u_drift, pulse_sequence_builder, +def create_pulse_sequence(u_ctrl, u_drift, *args): + # hacky hacky! + if len(args): + d = args[0] + + basis = ff.Basis.ggm(d) + H_c = (list(zip(basis[1:len(u_ctrl)+1], u_ctrl, [f'c{i}' for i in range(len(u_ctrl)+1)])) + + list(zip(basis[len(u_ctrl)+1:], u_drift, [f'd{i}' for i in range(d**2-len(u_ctrl)+1)]))) + H_n = (list(zip(basis[1:], np.ones((d**2-1, u_drift.shape[1]))))) + dt = np.full(u_ctrl.shape[-1], fill_value=0.32) + + return ff.PulseSequence(H_c, H_n, dt, basis=basis) + + +def finite_diff_infid(u_ctrl_central, u_drift, d, pulse_sequence_builder, spectral_noise_density, n_freq_samples=200, c_id=None, delta_u=1e-6): """ @@ -67,13 +83,15 @@ def finite_diff_infid(u_ctrl_central, u_drift, pulse_sequence_builder, gradient: shape (n_nop, n_dt, n_ctrl) """ - pulse = pulse_sequence_builder(u_ctrl_central, u_drift) + pulse = pulse_sequence_builder(u_ctrl_central, u_drift, d) all_id = pulse.c_oper_identifiers if c_id is None: c_id = all_id[:len(u_ctrl_central)] - omega = ff.util.get_sample_frequencies( - pulse=pulse, n_samples=n_freq_samples, spacing='log') - S = spectral_noise_density(omega) + + # Make sure we test for zero frequency case (possible convergence issues) + omega = ff.util.get_sample_frequencies(pulse=pulse, n_samples=n_freq_samples, spacing='log', + include_quasistatic=True) + spectrum = spectral_noise_density(omega) gradient = np.empty((len(pulse.n_coeffs), len(pulse.dt), len(c_id))) @@ -81,22 +99,22 @@ def finite_diff_infid(u_ctrl_central, u_drift, pulse_sequence_builder, for k in range(len(c_id)): u_plus = u_ctrl_central.copy() u_plus[k, g] += delta_u - pulse_plus = pulse_sequence_builder(u_plus, u_drift) - infid_plus = ff.numeric.infidelity(pulse=pulse_plus, spectrum=S, + pulse_plus = pulse_sequence_builder(u_plus, u_drift, pulse.d) + infid_plus = ff.numeric.infidelity(pulse=pulse_plus, spectrum=spectrum, omega=omega) u_minus = u_ctrl_central.copy() u_minus[k, g] -= delta_u - pulse_minus = pulse_sequence_builder(u_minus, u_drift) - infid_minus = ff.numeric.infidelity(pulse=pulse_minus, spectrum=S, + pulse_minus = pulse_sequence_builder(u_minus, u_drift, pulse.d) + infid_minus = ff.numeric.infidelity(pulse=pulse_minus, spectrum=spectrum, omega=omega) gradient[:, g, k] = (infid_plus - infid_minus) / 2 / delta_u return gradient -def analytic_gradient(u_ctrl, u_drift, pulse_sequence_builder, - spectral_noise_density, s_derivs=None, n_freq_samples=200, +def analytic_gradient(u_ctrl, u_drift, d, pulse_sequence_builder, + spectral_noise_density, n_coeffs_deriv=None, n_freq_samples=200, c_id=None, ctrl_amp_deriv=None): """ Parameters @@ -107,7 +125,7 @@ def analytic_gradient(u_ctrl, u_drift, pulse_sequence_builder, spectral_noise_density: function handle n_freq_samples: int c_id: List of string - s_derivs: shape (n_nops, n_ctrl, n_dt) + n_coeffs_deriv: shape (n_nops, n_ctrl, n_dt) ctrl_amp_deriv: function handle Returns @@ -115,21 +133,15 @@ def analytic_gradient(u_ctrl, u_drift, pulse_sequence_builder, gradient: shape (n_nop, n_dt, n_ctrl) """ - - pulse = pulse_sequence_builder(u_ctrl, u_drift) + pulse = pulse_sequence_builder(u_ctrl, u_drift, d) omega = ff.util.get_sample_frequencies( - pulse=pulse, n_samples=n_freq_samples, spacing='log') - S = spectral_noise_density(omega) + pulse=pulse, n_samples=n_freq_samples, spacing='log', include_quasistatic=True) + spectrum = spectral_noise_density(omega) gradient = ff.gradient.infidelity_derivative( - pulse=pulse, S=S, omega=omega, control_identifiers=c_id, - s_derivs=s_derivs) + pulse=pulse, spectrum=spectrum, omega=omega, control_identifiers=c_id, + n_coeffs_deriv=n_coeffs_deriv) if ctrl_amp_deriv is not None: # gradient shape (n_nops, n_dt, n_ctrl) # deriv_exchange_interaction shape (n_ctrl, n_dt) - gradient = np.einsum("agk,kg -> agk", gradient, - ctrl_amp_deriv(u_ctrl)) + gradient = np.einsum("agk,kg -> agk", gradient, ctrl_amp_deriv(u_ctrl)) return gradient - - -def relative_norm_difference(m, n): - return np.linalg.norm(m - n) / (np.linalg.norm(m) + np.linalg.norm(n)) * 2 diff --git a/tests/test_gradient.py b/tests/test_gradient.py index 3b331a1..f5418ab 100644 --- a/tests/test_gradient.py +++ b/tests/test_gradient.py @@ -20,41 +20,60 @@ # # Contact email: j.teske@fz-juelich.de # ============================================================================= - import numpy as np -import tests.gradient_testutil as grad_util -from tests import testutil +import filter_functions as ff +from tests import gradient_testutil, testutil + -np.random.seed(0) -initial_pulse = np.random.rand(grad_util.n_time_steps) -initial_pulse = np.expand_dims(initial_pulse, 0) -u_drift = 1. * np.ones(grad_util.n_time_steps) +class GradientTest(testutil.TestCase): -s_derivs = grad_util.deriv_2_exchange_interaction(eps=initial_pulse) -s_derivs = np.expand_dims(s_derivs, 0) + def test_gradient_calculation_variable_noise_coefficients(self): + initial_pulse = np.random.rand(gradient_testutil.n_time_steps) + initial_pulse = np.expand_dims(initial_pulse, 0) + u_drift = 1. * np.ones(gradient_testutil.n_time_steps) -fin_diff_grad = grad_util.finite_diff_infid( - u_ctrl_central=initial_pulse, u_drift=u_drift, - pulse_sequence_builder=grad_util.create_sing_trip_pulse_seq, - spectral_noise_density=grad_util.one_over_f_noise, n_freq_samples=200, - c_id=['control1'], delta_u=1e-4 -) + n_coeffs_deriv = gradient_testutil.deriv_2_exchange_interaction(eps=initial_pulse) + n_coeffs_deriv = np.expand_dims(n_coeffs_deriv, 0) + fin_diff_grad = gradient_testutil.finite_diff_infid( + u_ctrl_central=initial_pulse, u_drift=u_drift, d=2, + pulse_sequence_builder=gradient_testutil.create_sing_trip_pulse_seq, + spectral_noise_density=gradient_testutil.one_over_f_noise, n_freq_samples=200, + c_id=['control1'], delta_u=1e-4 + ) + ana_grad = gradient_testutil.analytic_gradient( + u_ctrl=initial_pulse, u_drift=u_drift, d=2, + pulse_sequence_builder=gradient_testutil.create_sing_trip_pulse_seq, + spectral_noise_density=gradient_testutil.one_over_f_noise, + c_id=['control1'], n_coeffs_deriv=np.ones_like(n_coeffs_deriv), + ctrl_amp_deriv=gradient_testutil.deriv_exchange_interaction + ) + self.assertArrayAlmostEqual(ana_grad, fin_diff_grad, rtol=1e-6, atol=1e-10) -ana_grad = grad_util.analytic_gradient( - u_ctrl=initial_pulse, u_drift=u_drift, - pulse_sequence_builder=grad_util.create_sing_trip_pulse_seq, - spectral_noise_density=grad_util.one_over_f_noise, - c_id=['control1'], s_derivs=np.ones_like(s_derivs), - ctrl_amp_deriv=grad_util.deriv_exchange_interaction -) + def test_gradient_calculation_random_pulse(self): + for d, n_dt in zip(testutil.rng.integers(2, 5, 5), testutil.rng.integers(2, 8, 5)): + u_ctrl = testutil.rng.normal(0, 1, (testutil.rng.integers(1, 4), n_dt)) + u_drift = testutil.rng.normal(0, 1, (d**2-1-len(u_ctrl), n_dt)) -class GradientTest(testutil.TestCase): + fin_diff_grad = gradient_testutil.finite_diff_infid( + u_ctrl_central=u_ctrl, u_drift=u_drift, d=d, + pulse_sequence_builder=gradient_testutil.create_pulse_sequence, + spectral_noise_density=gradient_testutil.one_over_f_noise, n_freq_samples=200, + c_id=[f'c{i}' for i in range(len(u_ctrl))], delta_u=1e-4 + ) + ana_grad = gradient_testutil.analytic_gradient( + u_ctrl=u_ctrl, u_drift=u_drift, d=d, + pulse_sequence_builder=gradient_testutil.create_pulse_sequence, + spectral_noise_density=gradient_testutil.one_over_f_noise, + c_id=[f'c{i}' for i in range(len(u_ctrl))], n_coeffs_deriv=None + ) + self.assertArrayAlmostEqual(ana_grad, fin_diff_grad, rtol=1e-6, atol=1e-10) - def test_gradient_calculation_variable_noise_coefficients(self): - relativ_diff = grad_util.relative_norm_difference( - fin_diff_grad, ana_grad) - self.assertLess(relativ_diff, 5e-9) + def test_raises(self): + pulse = testutil.rand_pulse_sequence(2, 3) + omega = ff.util.get_sample_frequencies(pulse, n_samples=13) + with self.assertRaises(ValueError): + ff.infidelity_derivative(pulse, 1/omega, omega, control_identifiers=['long string']) diff --git a/tests/test_precision.py b/tests/test_precision.py index 000f317..56bf519 100644 --- a/tests/test_precision.py +++ b/tests/test_precision.py @@ -23,8 +23,8 @@ """ import numpy as np -from scipy import linalg as sla from scipy import integrate +from scipy import linalg as sla import filter_functions as ff from filter_functions import analytic, numeric, util From 6dd8d49256497998334756aad2e44f3136182947 Mon Sep 17 00:00:00 2001 From: Tobias Hangleiter Date: Tue, 2 Mar 2021 16:10:37 +0100 Subject: [PATCH 05/11] Drop vectorized calculation and add comments Looped calculation is faster in almost all cases. --- filter_functions/gradient.py | 221 ++++++++--------------------------- 1 file changed, 46 insertions(+), 175 deletions(-) diff --git a/filter_functions/gradient.py b/filter_functions/gradient.py index 52a0faf..58f3fc0 100644 --- a/filter_functions/gradient.py +++ b/filter_functions/gradient.py @@ -198,15 +198,18 @@ def control_matrix_at_timestep_derivative( dt: Coefficients, eigvals: ndarray, eigvecs: ndarray, - basis: Basis, + basis_transformed, c_opers_transformed: ndarray, - n_opers: Sequence[Operator], + n_opers_transformed: ndarray, n_coeffs: Sequence[Coefficients], - n_coeffs_deriv: Optional[Sequence[Coefficients]] = None, - intermediates: Optional[Dict[str, ndarray]] = None + n_coeffs_deriv: Sequence[Coefficients], + phase_factor: ndarray, + integral: ndarray, + deriv_integral: ndarray, + ctrlmat_step: ndarray, + ctrlmat_expr: ContractExpression, ) -> Tuple[ndarray, ndarray]: - r""" - Calculate the control matrices and corresponding derivatives. + r"""Calculate the control matrices and corresponding derivatives. Calculate control matrices at each time step and the corresponding partial derivatives of those with respect to control strength at @@ -223,29 +226,39 @@ def control_matrix_at_timestep_derivative( eigvals: array_like, shape (n_dt, d) Eigenvalue vectors for each time pulse segment *g* with the first axis counting the pulse segment, i.e. - ``HD == array([D_0, D_1, ...])``. + ``D == array([D_0, D_1, ...])``. eigvecs: array_like, shape (n_dt, d, d) Eigenvector matrices for each time pulse segment *g* with the first axis counting the pulse segment, i.e. - ``HV == array([V_0, V_1, ...])``. - basis: Basis, shape (d**2, d, d) - The basis elements, in which the pulse control matrix will be - expanded. - c_opers_transformed: array_like, shape (n_dt, n_ctrl, d, d) + ``V == array([V_0, V_1, ...])``. + basis_transformed: array_like, shape (d**2, d, d) + The basis elements in which the pulse control matrix will be + expanded transformed by the eigenvectors. + c_opers_transformed: array_like, shape (n_ctrl, d, d) The control operators transformed into the eigenspace of the - control Hamiltonian. The drift operators are ignored, if - identifiers for accessible control operators are provided. - n_opers: array_like, shape (n_nops, d, d) - Noise operators :math:`B_\alpha`. - n_coeffs: array_like, shape (n_nops, n_dt) + control Hamiltonian. + n_opers_transformed: array_like, shape (n_nops, d, d) + The noise operators transformed into the eigenspace of the + control Hamiltonian. + n_coeffs: array_like, shape (n_nops,) The sensitivities of the system to the noise operators given by - *n_opers* at the given time step. - n_coeffs_deriv: array_like, shape (n_nops, n_ctrl, n_dt) + *n_opers_transformed* at the given time step. + n_coeffs_deriv: array_like, shape (n_nops, n_ctrl,) The derivatives of the noise susceptibilities by the control amplitudes. Defaults to None. - intermediates: Dict[str, ndarray], optional - Optional dictionary containing intermediate results of the - calculation of the control matrix. + phase_factor: array_like, shape (n_omega,) + The phase factor :math:`e^{i\omega t_{g-1}}`. + integral: array_like, shape (n_omega, d, d) + The integral during the time step appearing in the regular + control matrix. + deriv_integral: array_like, shape (n_omega, d, d, d, d) + An array to write the integral during the time step appearing in + the derivative into. + ctrlmat_step: array_like, shape (n_nops, d**2, n_omega) + An array to write the control matrix during the time step into. + ctrlmat_expr: ContractExpression + An :class:`opt_einsum.contract.ContractExpression` to compute + the control matrix during the time step. Returns ------- @@ -314,99 +327,6 @@ def control_matrix_at_timestep_derivative( \frac{\exp(\mathrm{i}(\omega + \Delta\omega_{nm})\Delta t_g) - 1} {\mathrm{i}(\omega + \Delta\omega_{nm})} """ - d = eigvecs.shape[-1] - n_dt = len(dt) - n_ctrl = len(c_opers_transformed[0]) - n_nops = len(n_opers) - n_omega = len(omega) - - # Cannot grab this from cached intermediates because those include the propagators - basis_transformed = numeric._transform_by_unitary(eigvecs[:, None], basis[None], - out=np.empty((n_dt, d**2, d, d), complex)) - if not intermediates: - n_opers_transformed = numeric._transform_hamiltonian(eigvecs, n_opers, - n_coeffs).swapaxes(0, 1) - exp_buf, integral = np.empty((2, n_omega, d, d), dtype=complex) - else: - n_opers_transformed = intermediates['n_opers_transformed'].swapaxes(0, 1) - first_order_integral = intermediates['first_order_integral'] - - deriv_integral = np.empty((n_dt, n_omega, d, d, d, d), dtype=complex) - ctrlmat_g = np.empty((n_dt, n_nops, d**2, n_omega), dtype=complex) - ctrlmat_g_deriv_s = np.empty((n_dt, n_nops, d**2, n_ctrl, n_omega), dtype=complex) - # Expression for calculating the control matrix during single time step - expr = oe.contract_expression('bcd,adc,hdc->abh', basis.shape, basis.shape, - (n_omega, d, d), optimize=[(0, 1), (0, 1)]) - for g in range(n_dt): - deriv_integral[g] = _derivative_integral(omega, eigvals[g], dt[g], deriv_integral[g]) - if not intermediates: - integral = numeric._first_order_integral(omega, eigvals[g], dt[g], exp_buf, integral) - else: - integral = first_order_integral[g] - - ctrlmat_g[g] = expr(basis_transformed[g], n_opers_transformed[g], integral, - out=ctrlmat_g[g]) - - if n_coeffs_deriv is not None: - # equivalent contraction: 'ah,a,ako->akho', but this faster - ctrlmat_g_deriv_s[g] = ( - (n_coeffs_deriv[..., g] / n_coeffs[:, g, None])[:, None, :, None] - * ctrlmat_g[g, :, :, None] - ) - - # Basically a tensor product - # K = np.einsum('taij,thkl->tahikjl', VBV, VHV) - # L = K.transpose(0, 1, 2, 4, 3, 6, 5) - K = util.tensor(n_opers_transformed[:, :, None], c_opers_transformed[:, None]) - L = util.tensor_transpose(K, (1, 0), [[d, d], [d, d]]) - k = np.diagonal(K.reshape(n_dt, n_nops, n_ctrl, d, d, d, d), 0, -2, -3) - l = np.diagonal(L.reshape(n_dt, n_nops, n_ctrl, d, d, d, d), 0, -2, -3) - i1 = np.diagonal(deriv_integral, 0, -2, -3) - i2 = np.diagonal(deriv_integral, 0, -1, -4) - # Reshaping in F-major is faster than not (~ factor 2-4) - M = np.einsum( - 'tahpm,topm->tahop', - l.reshape(n_dt, n_nops, n_ctrl, d**2, d, order='F'), - i1.reshape(n_dt, n_omega, d**2, d, order='F') - ).reshape(n_dt, n_nops, n_ctrl, n_omega, d, d, order='F') - if d == 2: - # Life a bit simpler - mask = np.eye(d, dtype=bool) - M[..., mask] -= M[..., mask][..., ::-1] - M[..., ~mask] *= 2 - else: - M -= np.einsum( - 'tahpn,topn->tahop', - k.swapaxes(-2, -3).reshape(n_dt, n_nops, n_ctrl, d**2, d, order='F'), - i2.reshape(n_dt, n_omega, d**2, d, order='F') - ).reshape(n_dt, n_nops, n_ctrl, n_omega, d, d, order='F').swapaxes(-1, -2) - - # Expand in basis transformed to eigenspace - ctrlmat_g_deriv = oe.contract('tjnk,tahokn->tajho', 1j*basis_transformed, M) - - if n_coeffs_deriv is not None: - ctrlmat_g_deriv += ctrlmat_g_deriv_s - - return ctrlmat_g, ctrlmat_g_deriv - - -def control_matrix_at_timestep_derivative_loop( - omega: Coefficients, - dt: Coefficients, - eigvals: ndarray, - eigvecs: ndarray, - basis_transformed, - c_opers_transformed: ndarray, - n_opers_transformed, - n_coeffs: Sequence[Coefficients], - n_coeffs_deriv: Sequence[Coefficients], - ctrlmat_step, - phase_factor, - deriv_integral, - integral, - ctrlmat_expr -) -> Tuple[ndarray, ndarray]: - """".""" d = len(eigvecs) d2 = d**2 n_ctrl = len(c_opers_transformed) @@ -444,7 +364,8 @@ def control_matrix_at_timestep_derivative_loop( i2.reshape(n_omega, d2, d, order='F') ).reshape(n_nops, n_ctrl, n_omega, d, d, order='F').swapaxes(-1, -2) - # Expand in basis transformed to eigenspace + # Expand in basis transformed to eigenspace. Include phase factor and + # factor 1j here to make use of optimized contraction order ctrlmat_step_deriv = oe.contract('o,jnk,ahokn->ajho', phase_factor, 1j*basis_transformed, M, optimize=[(1, 2), (0, 1)]) @@ -550,62 +471,9 @@ def calculate_derivative_of_control_matrix_from_scratch( See Also -------- - :func:`liouville_derivative` - :func:`control_matrix_at_timestep_derivative` + :func:`_liouville_derivative` + :func:`_control_matrix_at_timestep_derivative` """ - # Distinction between control and drift operators and only calculate the - # derivatives in control direction - try: - idx = util.get_indices_from_identifiers(all_identifiers, control_identifiers) - except ValueError as err: - raise ValueError('Given control identifiers have to be a subset of (drift+control) ' + - 'Hamiltonian!') from err - - c_opers_transformed = numeric._transform_hamiltonian(eigvecs, c_opers[idx]).swapaxes(0, 1) - propagators_liouville = superoperator.liouville_representation(propagators[:-1], basis) - propagators_liouville_deriv = liouville_derivative(dt, propagators, basis, eigvecs, eigvals, - c_opers_transformed) - - ctrlmat_g, ctrlmat_g_deriv = control_matrix_at_timestep_derivative( - omega, dt, eigvals, eigvecs, basis, c_opers_transformed, n_opers, n_coeffs, n_coeffs_deriv, - intermediates - ) - - if intermediates: - phase_factors = intermediates['phase_factors'].T - else: - phase_factors = util.cexp(np.multiply.outer(omega, t[:-1])) - - # Equivalent (but slower) einsum contraction for first term: - # ctrlmat_deriv = np.einsum('ot,tajho,tjk->hotak', - # phase_factors, ctrlmat_g_deriv, propagators_liouville, - # optimize=['einsum_path', (0, 2), (0, 1)]) - ctrlmat_deriv = (ctrlmat_g_deriv.transpose(3, 4, 0, 1, 2) - @ (phase_factors[..., None, None] * propagators_liouville)) - ctrlmat_deriv += np.einsum('ot,tajo,thsjk->hosak', - phase_factors[:, 1:], ctrlmat_g[1:], propagators_liouville_deriv, - optimize=['einsum_path', (0, 1), (0, 1)]) - - return ctrlmat_deriv - - -def calculate_derivative_of_control_matrix_from_scratch_loop( - omega: Coefficients, - propagators: ndarray, - eigvals: ndarray, - eigvecs: ndarray, - basis: Basis, - t: Coefficients, - dt: Coefficients, - n_opers: Sequence[Operator], - n_coeffs: Sequence[Coefficients], - c_opers: Sequence[Operator], - all_identifiers: Sequence[str], - control_identifiers: Optional[Sequence[str]] = None, - n_coeffs_deriv: Optional[Sequence[Coefficients]] = None, - intermediates: Dict[str, ndarray] = None -) -> ndarray: - r""".""" # Distinction between control and drift operators and only # calculate the derivatives in control direction try: @@ -620,6 +488,7 @@ def calculate_derivative_of_control_matrix_from_scratch_loop( n_nops = len(n_opers) n_omega = len(omega) + # Precompute some transformations or grab from cache if possible basis_transformed = numeric._transform_by_unitary(eigvecs[:, None], basis[None], out=np.empty((n_dt, d**2, d, d), complex)) c_opers_transformed = numeric._transform_hamiltonian(eigvecs, c_opers[idx]).swapaxes(0, 1) @@ -637,6 +506,8 @@ def calculate_derivative_of_control_matrix_from_scratch_loop( deriv_integral = np.empty((n_omega, d, d, d, d), dtype=complex) ctrlmat_deriv = np.empty((n_ctrl, n_omega, n_dt, n_nops, d**2), dtype=complex) ctrlmat_step = np.empty((n_dt, n_nops, d**2, n_omega), dtype=complex) + # Optimized expression that is passed to control_matrix_at_timestep_derivative + # in each iteration ctrlmat_expr = oe.contract_expression('o,icd,adc,odc->aio', (len(omega),), basis.shape, n_opers.shape, (len(omega), d, d), optimize=[(0, 3), (0, 1), (0, 1)]) @@ -644,21 +515,21 @@ def calculate_derivative_of_control_matrix_from_scratch_loop( if intermediates is None: integral = numeric._first_order_integral(omega, eigvals[g], dt[g], exp_buf, integral) else: - # Grab from cache integral = intermediates['first_order_integral'][g] - phase_factor = util.cexp(omega * t[g]) n_coeff_deriv = n_coeffs_deriv if n_coeffs_deriv is None else n_coeffs_deriv[:, :, g] ctrlmat_step[g], ctrlmat_step_deriv = control_matrix_at_timestep_derivative_loop( omega, dt[g], eigvals[g], eigvecs[g], basis_transformed[g], c_opers_transformed[g], - n_opers_transformed[g], n_coeffs[:, g], n_coeff_deriv, ctrlmat_step[g], - phase_factor, deriv_integral, integral, ctrlmat_expr + n_opers_transformed[g], n_coeffs[:, g], n_coeff_deriv, util.cexp(omega*t[g]), integral, + deriv_integral, ctrlmat_step[g], ctrlmat_expr ) + # Phase factor already part of ctrlmat_step_deriv ctrlmat_deriv[:, :, g] = (ctrlmat_step_deriv.transpose(2, 3, 0, 1) @ propagators_liouville[g]) # opt_einsum a lot faster here + # Phase factor again already part of ctrlmat_step ctrlmat_deriv += oe.contract('tajo,thsjk->hosak', ctrlmat_step[1:], propagators_liouville_deriv) From e0f02ade3c2a9caa8cb950a71228abf4667f0a8a Mon Sep 17 00:00:00 2001 From: Tobias Hangleiter Date: Tue, 2 Mar 2021 16:12:27 +0100 Subject: [PATCH 06/11] Make lower-level functions private liouville_derivative and control_matrix_at_timestep_derivative do not really make sense on their own. --- filter_functions/gradient.py | 25 ++++++++++++------------- 1 file changed, 12 insertions(+), 13 deletions(-) diff --git a/filter_functions/gradient.py b/filter_functions/gradient.py index 58f3fc0..7216044 100644 --- a/filter_functions/gradient.py +++ b/filter_functions/gradient.py @@ -44,11 +44,6 @@ Functions --------- -:func:`liouville_derivative` - Calculate the derivatives of the control propagators in Liouville - representation. -:func:`control_matrix_at_timestep_derivative` - Calculate the control matrices and corresponding derivatives. :func:`calculate_derivative_of_control_matrix_from_scratch` Calculate the derivative of the control matrix from scratch. :func:`calculate_canonical_filter_function_derivative` @@ -60,14 +55,14 @@ import numpy as np import opt_einsum as oe +from opt_einsum.contract import ContractExpression from numpy import ndarray from . import numeric, superoperator, util from .basis import Basis from .types import Coefficients, Operator -__all__ = ['liouville_derivative', 'control_matrix_at_timestep_derivative', - 'calculate_derivative_of_control_matrix_from_scratch', +__all__ = ['calculate_derivative_of_control_matrix_from_scratch', 'calculate_filter_function_derivative', 'infidelity_derivative'] @@ -105,8 +100,8 @@ def _derivative_integral(E: Coefficients, eigvals: Coefficients, dt: float, return out -def liouville_derivative(dt: Coefficients, propagators: ndarray, basis: Basis, eigvecs: ndarray, - eigvals: ndarray, c_opers_transformed: ndarray) -> ndarray: +def _liouville_derivative(dt: Coefficients, propagators: ndarray, basis: Basis, eigvecs: ndarray, + eigvals: ndarray, c_opers_transformed: ndarray) -> ndarray: r""" Calculate the derivatives of the control propagators in Liouville representation. @@ -193,7 +188,7 @@ def liouville_derivative(dt: Coefficients, propagators: ndarray, basis: Basis, e return liouville_deriv -def control_matrix_at_timestep_derivative( +def _control_matrix_at_timestep_derivative( omega: Coefficients, dt: Coefficients, eigvals: ndarray, @@ -500,8 +495,8 @@ def calculate_derivative_of_control_matrix_from_scratch( n_opers_transformed = intermediates['n_opers_transformed'].swapaxes(0, 1) propagators_liouville = superoperator.liouville_representation(propagators[:-1], basis) - propagators_liouville_deriv = liouville_derivative(dt, propagators, basis, eigvecs, eigvals, - c_opers_transformed) + propagators_liouville_deriv = _liouville_derivative(dt, propagators, basis, eigvecs, eigvals, + c_opers_transformed) deriv_integral = np.empty((n_omega, d, d, d, d), dtype=complex) ctrlmat_deriv = np.empty((n_ctrl, n_omega, n_dt, n_nops, d**2), dtype=complex) @@ -519,7 +514,11 @@ def calculate_derivative_of_control_matrix_from_scratch( n_coeff_deriv = n_coeffs_deriv if n_coeffs_deriv is None else n_coeffs_deriv[:, :, g] - ctrlmat_step[g], ctrlmat_step_deriv = control_matrix_at_timestep_derivative_loop( + # ctrlmat_step is computed from scratch because the quantity in + # the cache (from numeric.calculate_control_matrix_from_scratch) + # contains the Liouville propagators already and it is expensive + # to remove them by multiplying with the transpose. + ctrlmat_step[g], ctrlmat_step_deriv = _control_matrix_at_timestep_derivative( omega, dt[g], eigvals[g], eigvecs[g], basis_transformed[g], c_opers_transformed[g], n_opers_transformed[g], n_coeffs[:, g], n_coeff_deriv, util.cexp(omega*t[g]), integral, deriv_integral, ctrlmat_step[g], ctrlmat_expr From 767d685560efc97932b79c0ba99da154c286ab44 Mon Sep 17 00:00:00 2001 From: Tobias Hangleiter Date: Tue, 2 Mar 2021 16:47:04 +0100 Subject: [PATCH 07/11] Remove erroneously added old file --- filter_functions/gradient_.py | 1855 --------------------------------- pytest.ini | 6 +- 2 files changed, 5 insertions(+), 1856 deletions(-) delete mode 100644 filter_functions/gradient_.py diff --git a/filter_functions/gradient_.py b/filter_functions/gradient_.py deleted file mode 100644 index 4ef2dfa..0000000 --- a/filter_functions/gradient_.py +++ /dev/null @@ -1,1855 +0,0 @@ -# -*- coding: utf-8 -*- -# ============================================================================= -# filter_functions -# Copyright (C) 2019 Quantum Technology Group, RWTH Aachen University -# -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program. If not, see . -# -# This module is an extension of the filter_functions package written by -# Tobias Hangleiter. Its implementation was center of the Bachelor thesis -# 'Filter Function Derivative for Quantum Optimal Control' by Isabel Le. -# -# Contact email: isabel.le@rwth-aachen.de -# -# The code has been extended by Julian Teske such that the noise -# susceptibilities or noise coefficients can depend on the control -# amplitudes as well. -# -# Contact email: j.teske@fz-juelich.de -# ============================================================================= -""" -This module implements functions to calculate filter function and -infidelity derivatives. - -Throughout this documentation the following notation will be used: - - n_dt denotes the number of time steps, - - n_cops the number of all control operators, - - n_ctrl the number of accessible control operators (if identifiers - are provided, otherwise n_ctrl=n_cops), - - n_nops the number of noise operators, - - n_omega the number of frequency samples, and - - d the dimension of the system. - -Functions ---------- -:func:`liouville_derivative` - Calculate the derivatives of the control propagators in Liouville - representation. -:func:`control_matrix_at_timestep_derivative` - Calculate the control matrices and corresponding derivatives. -:func:`calculate_derivative_of_control_matrix_from_scratch` - Calculate the derivative of the control matrix from scratch. -:func:`calculate_canonical_filter_function_derivative` - Compute the filter function derivative from the control matrix. -:func:`infidelity_derivative` - Calculate the infidelity derivative. -""" -from typing import Callable, Dict, Optional, Sequence, Tuple, Union - -import numpy as np -import opt_einsum as oe -from numpy import ndarray -from scipy.integrate import trapz - -from . import numeric, superoperator, util -from .basis import Basis -from .types import Coefficients, Operator - -__all__ = ['liouville_derivative', 'control_matrix_at_timestep_derivative', - 'calculate_derivative_of_control_matrix_from_scratch', - 'calculate_filter_function_derivative', 'infidelity_derivative'] - - -def _derivative_integral_(E, eigvals, dt, deriv_integral): - """Calculate I_nmpq, but return I_pqnm""" - # Any combination of \omega_p-\omega_q (dR_g) - dE = np.subtract.outer(eigvals, eigvals) - mask_dE = np.abs(dE) < 1e-7 - # Any combination of \omega+\omega_n-\omega_m (dR_g) - EdE = np.add.outer(E, dE) - mask_EdE = np.abs(EdE) < 1e-7 - # Any combination of \omega+\omega_n-\omega_m+\omega_p-\omega_q (dR_g) - EdEdE = np.add.outer(EdE, dE[~mask_dE]) - - # calculation if omega_diff = 0 - tmp1 = np.divide(util.cexp(EdE*dt), EdE, where=~mask_EdE) - tmp2 = tmp1 - np.divide(1, EdE, where=~mask_EdE) - - tmp1 *= -1j * dt - tmp1 += np.divide(tmp2, EdE, where=~mask_EdE) - tmp1[mask_EdE] = dt**2 / 2 - - # tmp1 = util.cexp(EdE*dt) / EdE - # tmp2 = tmp1 - 1 / EdE - # tmp1 = -1j * dt * tmp1 + tmp2 / EdE - - deriv_integral[:, mask_dE] = tmp1[:, None] - - # # calculation if omega_diff != 0 - # tmp1 = np.divide(1 - util.cexp(EdEdE*dt), EdEdE, where=~mask_ndE_nEdE) - # tmp1[mask_ndE_EdE] = -1j * dt - # tmp1 += tmp2[:, None, None] - # tmp1 = np.divide(tmp1, dE, where=~mask_dE, out=tmp1) - # deriv_integral[:, ~mask_dE] = tmp1[:, ~mask_dE] - # calculation if omega_diff != 0 - - # pq on last axis in EdEdE, pq on 1,2nd axes in deriv_integral - deriv_integral[:, ~mask_dE] = ((1 - util.cexp(EdEdE*dt)) / EdEdE / dE[~mask_dE]).transpose(0, 3, 1, 2) - deriv_integral[:, ~mask_dE] += np.divide.outer(tmp2, dE[~mask_dE]).transpose(0, 3, 1, 2) - return deriv_integral - - -def _derivative_integral(E, eigvals, dt, out): - dE = np.subtract.outer(eigvals, eigvals) - mask_dE = np.abs(dE) < 1e-7 - EdE = np.add.outer(E, dE) - mask_EdE = np.abs(EdE) < 1e-7 - EdEdE = np.add.outer(EdE, dE[~mask_dE]) - mask_EdEdE = np.abs(EdEdE) < 1e-7 - - # Omega_pq == 0 - tmp1 = np.divide(util.cexp(EdE*dt), EdE, where=~mask_EdE) - tmp2 = tmp1 - np.divide(1, EdE, where=~mask_EdE) - tmp2[mask_EdE] = 1j * dt - - tmp1 *= -1j * dt - tmp1 += np.divide(tmp2, EdE, where=~mask_EdE) - tmp1[mask_EdE] = dt**2 / 2 - - out[:, mask_dE] = tmp1[:, None] - - # Omega_pq != 0 - tmp1 = np.divide(1 - util.cexp(EdEdE*dt), EdEdE, where=~mask_EdEdE) - tmp1[mask_EdEdE] = -1j * dt - tmp1 += tmp2[..., None] - - out[:, ~mask_dE] = (tmp1 / dE[~mask_dE]).transpose(0, 3, 1, 2) - - return out - - -def liouville_derivative( - dt: Coefficients, - propagators: ndarray, - basis: Basis, - eigvecs: ndarray, - eigvals: ndarray, - transf_control_operators: ndarray) -> ndarray: - r""" - Calculate the derivatives of the control propagators in Liouville - representation. - - Parameters - ---------- - dt : array_like, shape (n_dt) - Sequence duration, i.e. for the :math:`g`-th pulse - :math:`t_g - t_{g-1}`. - propagators : array_like, shape (n_dt+1, d, d) - The propagators :math:`Q_g = P_g P_{g-1}\cdots P_0` as a (d, d) array - with *d* the dimension of the Hilbert space. - basis : Basis, shape (d**2, d, d) - The basis elements, in which the pulse control matrix will be expanded. - eigvecs : array_like, shape (n_dt, d, d) - Eigenvector matrices for each time pulse segment *g* with the first - axis counting the pulse segment, i.e. - ``HV == array([V_0, V_1, ...])``. - eigvals : array_like, shape (n_dt, d) - Eigenvalue vectors for each time pulse segment *g* with the first - axis counting the pulse segment, i.e. - ``HD == array([D_0, D_1, ...])``. - transf_control_operators : array_like, shape (n_dt, c_ctrl, d, d) - The control operators transformed into the eigenspace of the control - Hamiltonian. The drift operators are ignored, if identifiers for - accessible control operators are provided. - - Returns - ------- - dL: array_like, shape (n_dt, n_ctrl, n_dt, d**2, d**2) - The derivative of the control propagators in Liouville representation - :math:`\frac{\partial \mathcal{Q}_{jk}^{(g)}} - {\partial u_h(t_{g^\prime})}`. - The array's indexing has shape :math:`(g,h,g^\prime,j,k)`. - - Notes - ----- - The derivatives of the control propagators in Liouville representation are - calculated according to - - .. math:: - - \frac{\partial\mathcal{Q}_{jk}^{(g-1)}}{\partial u_h(t_{g^\prime})} &= - \Theta_{g-1}(g^\prime) \mathrm{tr}\Big( - \frac{\partial U_c^\dagger(t_{g-1},0)}{\partial u_h(t_{g^\prime})} - C_j U_c(t_{g-1},0) C_k\Big)\\ - &+\Theta_{g-1}(g^\prime)\mathrm{tr}\Big(U_c^\dagger(t_{g-1},0)C_j - \frac{\partial U_c(t_{g-1},0)} - {\partial u_h(t_{g^\prime})}C_k - \Big), - - where :math:`\Theta_{g-1}(g^\prime)` being 1 if :math:`g^\prime < g-1` and - zero otherwise. - - """ - n = len(dt) - d = basis.shape[-1] - - # Allocate memory - A_mat = np.empty((d, d), dtype=complex) - partial_U = np.empty((n, transf_control_operators.shape[1], d, d), - dtype=complex) - deriv_U = np.zeros((n, n, transf_control_operators.shape[1], d, d), - dtype=complex) - - for g in (range(n)): - omega_diff = np.subtract.outer(eigvals[g], eigvals[g]) - mask = (abs(omega_diff) < 1e-10) - A_mat[mask] = dt[g] # if the integral diverges - A_mat[~mask] = (util.cexp(omega_diff[~mask]*dt[g]) - 1) \ - / (1j * omega_diff[~mask]) - # Calculate dU(t_g,t_{g-1})/du_h within one time step - partial_U[g] = -1j * np.einsum( - 'ia,ja,jk,hkl,kl,ml->him', propagators[g + 1], - propagators[g].conj(), eigvecs[g], transf_control_operators[g], - A_mat, eigvecs[g].conj(), optimize=['einsum_path', (3, 4), (0, 1), - (0, 3), (0, 1), (0, 1)]) - # Calculate the whole propagator derivative - for g_prime in range(g+1): # condition g' <= g-1 - deriv_U[g, g_prime] = np.einsum( - 'ij,kj,hkl,lm->him', propagators[g + 1], - propagators[g_prime + 1].conj(), partial_U[g_prime], - propagators[g_prime], optimize=['einsum_path', (0, 1), - (0, 1), (0, 1)]) - - # Now calculate derivative of Liouville representation - # Calculate traces first to save memory - sum1 = np.einsum( - 'tshba,jbc,tcd,kda->thsjk', deriv_U.conj(), basis, propagators[1:], - basis, optimize=['einsum_path', (1, 2), (0, 2), (0, 1)]) - sum2 = np.einsum( - 'tba,jbc,tshcd,kda->thsjk', propagators[1:].conj(), basis, deriv_U, - basis, optimize=['einsum_path', (0, 1), (0, 2), (0, 1)]) - liouville_deriv = sum1 + sum2 - - return liouville_deriv - - -def liouville_derivative_faster( - dt: Coefficients, - propagators: ndarray, - basis: Basis, - eigvecs: ndarray, - eigvals: ndarray, - transf_control_operators: ndarray) -> ndarray: - """.""" - n = len(dt) - d = basis.shape[-1] - - # Allocate memory - A_mat = np.empty((d, d), dtype=complex) - omega_diff = np.empty((d, d), dtype=float) - partial_U = np.empty((n, transf_control_operators.shape[1], d, d), - dtype=complex) - deriv_U = np.zeros((n, n, transf_control_operators.shape[1], d, d), - dtype=complex) - - for g in (range(n)): - omega_diff = np.subtract.outer(eigvals[g], eigvals[g], out=omega_diff) - mask = omega_diff == 0 - A_mat[mask] = dt[g] # if the integral diverges - A_mat[~mask] = (util.cexp(omega_diff[~mask]*dt[g]) - 1) / (1j * omega_diff[~mask]) - # Calculate dU(t_g,t_{g-1})/du_h within one time step - partial_U[g] = -1j * (propagators[g+1] @ propagators[g].conj().T @ eigvecs[g] @ - (transf_control_operators[g] * A_mat) @ eigvecs[g].conj().T) - # # Calculate the whole propagator derivative - deriv_U[g, :g+1] = (propagators[g+1] @ propagators[1:g+2].conj().swapaxes(-1, -2) @ - partial_U[:g+1].swapaxes(0, 1) @ propagators[:g+1]).swapaxes(0, 1) - - # Now calculate derivative of Liouville representation - # Calculate traces first to save memory - liouville_deriv = np.einsum( - 'tshba,jbc,tcd,kda->thsjk', deriv_U.conj(), basis, propagators[1:], - basis, optimize=['einsum_path', (1, 2), (0, 2), (0, 1)]) - - return 2*liouville_deriv.real - - -def liouville_derivative_vectorized( - dt: Coefficients, - propagators: ndarray, - basis: Basis, - eigvecs: ndarray, - eigvals: ndarray, - VHV: ndarray -) -> ndarray: - """".""" - n, d = eigvecs.shape[:2] - # omega_i - omega_j - omega_diff = eigvals[:, :, None] - eigvals[:, None] - dt_broadcast = np.broadcast_to(dt[:, None, None], omega_diff.shape) - # mask = omega_diff == 0 - mask = np.broadcast_to(np.eye(d, dtype=bool), omega_diff.shape) - A_mat = np.empty(omega_diff.shape, dtype=complex) - A_mat[mask] = dt_broadcast[mask] - # A_mat[~mask] = (util.cexp((omega_diff[~mask]*dt_broadcast[~mask])) - 1) / (1j * omega_diff[~mask]) - A_mat[~mask] = np.expm1(1j*(omega_diff[~mask]*dt_broadcast[~mask])) / (1j * omega_diff[~mask]) - U_deriv = -1j * (propagators[1:] - @ propagators[:-1].conj().swapaxes(-1, -2) - @ eigvecs - @ (A_mat * VHV.swapaxes(0, 1)) - @ eigvecs.conj().swapaxes(-1, -2)) - - # Calculate the whole propagator derivative - propagators_deriv = np.zeros((VHV.shape[1], n-1, n, d, d), dtype=complex) - for g in range(n-1): - propagators_deriv[:, g, :g+1] = (propagators[g+1] - @ propagators[1:g+2].conj().swapaxes(-1, -2) - @ U_deriv[:, :g+1] - @ propagators[:g+1]) - - # liouville_deriv = np.einsum('htsba,tjkba->thsjk', propagators_deriv.conj(), - # (basis @ propagators[1:-1, None])[:, :, None] @ basis).real - liouville_deriv = oe.contract('htsba,jbc,tcd,kda->thsjk', - propagators_deriv.conj(), basis, propagators[1:-1], basis, - optimize=[(0, 2), (0, 2), (0, 1)]).real - liouville_deriv *= 2 - - return liouville_deriv - - -def liouville_derivative_vectorized_completely( - dt: Coefficients, - propagators: ndarray, - basis: Basis, - eigvecs: ndarray, - eigvals: ndarray, - VHV: ndarray -) -> ndarray: - n, d = eigvecs.shape[:2] - # omega_i - omega_j - omega_diff = eigvals[:, :, None] - eigvals[:, None] - dt_broadcast = np.broadcast_to(dt[:, None, None], omega_diff.shape) - # mask = omega_diff == 0 - mask = np.broadcast_to(np.eye(d, dtype=bool), omega_diff.shape) - A_mat = np.empty(omega_diff.shape, dtype=complex) - A_mat[mask] = dt_broadcast[mask] - A_mat[~mask] = np.expm1(1j*(omega_diff[~mask]*dt_broadcast[~mask])) / (1j * omega_diff[~mask]) - U_deriv = -1j * (propagators[1:] - @ propagators[:-1].conj().swapaxes(-1, -2) - @ eigvecs - @ (A_mat * VHV.swapaxes(0, 1)) - @ eigvecs.conj().swapaxes(-1, -2)) - - # opt_einsum performs some magic with intermediate terms here, faster than explicit loop - propagators_deriv = oe.contract('tab,scb,hscd,sde->htsae', propagators[1:-1], - propagators[1:].conj(), U_deriv, propagators[:-1], - optimize=[(1, 2), (1, 2), (0, 1)]) - # Derivative is zero for times in the future of the control step - propagators_deriv[:, ~np.tri(n - 1, n, dtype=bool)] = 0 - - liouville_deriv = oe.contract('htsba,jbc,tcd,kda->thsjk', - propagators_deriv.conj(), basis, propagators[1:-1], basis, - optimize=[(0, 2), (0, 2), (0, 1)]).real - # oe.contract('xtba,jbc,tcd,kda->xtjk', propagators_deriv.transpose(0,2,1,3,4).reshape(5*100,99,4,4, order='F').conj(), basis, propagators[1:-1], basis, optimize=p).real.reshape(5,100,99,16,16,order='F') - liouville_deriv *= 2 - - return liouville_deriv - - -def liouville_derivative_vectorized_loop( - dt: Coefficients, - propagators: ndarray, - basis: Basis, - eigvecs: ndarray, - eigvals: ndarray, - c_opers_transformed: ndarray -) -> ndarray: - """".""" - n, d = eigvecs.shape[:2] - # omega_i - omega_j - omega_diff = eigvals[:, :, None] - eigvals[:, None] - dt_broadcast = np.broadcast_to(dt[:, None, None], omega_diff.shape) - # mask = omega_diff == 0 - mask = np.broadcast_to(np.eye(d, dtype=bool), omega_diff.shape) - integral = np.empty(omega_diff.shape, dtype=complex) - integral[mask] = dt_broadcast[mask] - # integral[~mask] = (util.cexp((omega_diff[~mask]*dt_broadcast[~mask])) - 1) / (1j * omega_diff[~mask]) - integral[~mask] = (np.expm1(1j*(omega_diff[~mask]*dt_broadcast[~mask])) - / (1j * omega_diff[~mask])) - - return -1j * (propagators[1:] - @ propagators[:-1].conj().swapaxes(-1, -2) - @ eigvecs - @ (integral * c_opers_transformed.swapaxes(0, 1)) - @ eigvecs.conj().swapaxes(-1, -2)) - - -def liouville_derivative_vectorized_completely_loop(propagator_deriv, basis, propagator, expr, - out): - liouville_deriv = expr(propagator_deriv.conj(), basis, propagator, basis, out=out) - liouville_deriv.real *= 2 - return liouville_deriv - - -def propagators_derivative_vectorized_completely( - dt: Coefficients, - propagators: ndarray, - basis: Basis, - eigvecs: ndarray, - eigvals: ndarray, - c_opers_transformed: ndarray -) -> ndarray: - - n, d = eigvecs.shape[:2] - # omega_i - omega_j - omega_diff = eigvals[:, :, None] - eigvals[:, None] - dt_broadcast = np.broadcast_to(dt[:, None, None], omega_diff.shape) - # mask = omega_diff == 0 - mask = np.broadcast_to(np.eye(d, dtype=bool), omega_diff.shape) - A_mat = np.empty(omega_diff.shape, dtype=complex) - A_mat[mask] = dt_broadcast[mask] - A_mat[~mask] = np.expm1(1j*(omega_diff[~mask]*dt_broadcast[~mask])) / (1j * omega_diff[~mask]) - U_deriv = -1j * (propagators[1:] - @ propagators[:-1].conj().swapaxes(-1, -2) - @ eigvecs - @ (A_mat * c_opers_transformed.swapaxes(0, 1)) - @ eigvecs.conj().swapaxes(-1, -2)) - - # opt_einsum performs some magic with intermediate terms here, faster than explicit loop - propagators_deriv = oe.contract('tab,scb,hscd,sde->htsae', propagators[1:-1], - propagators[1:].conj(), U_deriv, propagators[:-1], - optimize=[(1, 2), (1, 2), (0, 1)]) - # Derivative is zero for times in the future of the control step - propagators_deriv[:, ~np.tri(n - 1, n, dtype=bool)] = 0 - return propagators_deriv - - -def propagator_derivative( - dt: Coefficients, - propagators: ndarray, - eigvecs: ndarray, - eigvals: ndarray, - transf_control_operators: ndarray) -> ndarray: - """".""" - n, d = eigvecs.shape[:2] - # omega_i - omega_j - omega_diff = eigvals[:, :, None] - eigvals[:, None] - dt_broadcast = np.broadcast_to(dt[:, None, None], omega_diff.shape) - # mask = omega_diff == 0 - mask = np.broadcast_to(np.eye(d, dtype=bool), omega_diff.shape) - A_mat = np.empty(omega_diff.shape, dtype=complex) - A_mat[mask] = dt_broadcast[mask] - # A_mat[~mask] = (util.cexp((omega_diff[~mask]*dt_broadcast[~mask])) - 1) / (1j * omega_diff[~mask]) - A_mat[~mask] = np.expm1(1j*(omega_diff[~mask]*dt_broadcast[~mask])) / (1j * omega_diff[~mask]) - partial_U = -1j * (propagators[1:] @ propagators[:-1].conj().swapaxes(-1, -2) - @ eigvecs - @ (transf_control_operators.swapaxes(0, 1) * A_mat) - @ eigvecs.conj().swapaxes(-1, -2)) - - # Calculate the whole propagator derivative - derivative = np.zeros((transf_control_operators.shape[1], n, n, d, d), dtype=complex) - for g in range(n): - derivative[:, g, :g+1] = (propagators[g+1] @ propagators[1:g+2].conj().swapaxes(-1, -2) @ - partial_U[:, :g+1] @ propagators[:g+1]) - - return derivative - - -def propagator_derivative_factor( - dt: Coefficients, - propagators: ndarray, - eigvecs: ndarray, - eigvals: ndarray, - VHV: ndarray) -> ndarray: - """".""" - n, d = eigvecs.shape[:2] - # omega_i - omega_j - omega_diff = eigvals[:, :, None] - eigvals[:, None] - dt_broadcast = np.broadcast_to(dt[:, None, None], omega_diff.shape) - # mask = omega_diff == 0 - mask = np.broadcast_to(np.eye(d, dtype=bool), omega_diff.shape) - A_mat = np.empty(omega_diff.shape, dtype=complex) - A_mat[mask] = dt_broadcast[mask] - # A_mat[~mask] = (util.cexp((omega_diff[~mask]*dt_broadcast[~mask])) - 1) / (1j * omega_diff[~mask]) - A_mat[~mask] = np.expm1(1j*(omega_diff[~mask]*dt_broadcast[~mask])) / (1j * omega_diff[~mask]) - partial_U = -1j * (propagators[1:] @ propagators[:-1].conj().swapaxes(-1, -2) - @ eigvecs - @ (VHV.swapaxes(0, 1) * A_mat) - @ eigvecs.conj().swapaxes(-1, -2)) - - # Calculate the whole propagator derivative - derivative = np.zeros((VHV.shape[1], n, n, d, d), dtype=complex) - for g in range(n): - derivative[:, g, :g+1] = (propagators[1:g+2].conj().swapaxes(-1, -2) - @ partial_U[:, :g+1] - @ propagators[:g+1]) - - return derivative - - -def control_matrix_at_timestep_derivative( - omega: Coefficients, - dt: Coefficients, - eigvals: ndarray, - eigvecs: ndarray, - basis: Basis, - n_opers: Sequence[Operator], - n_coeffs: Sequence[Coefficients], - transf_control_operators: ndarray, - s_derivs: Optional[Sequence[Coefficients]] = None) -> dict: - r""" - Calculate the control matrices and corresponding derivatives. - - Calculate control matrices at each time step and the corresponding partial - derivatives of those with respect to control strength at each time step. - - Parameters - ---------- - omega : array_like, shape (n_omega) - Frequencies, at which the pulse control matrix is to be evaluated. - dt : array_like, shape (n_dt) - Sequence duration, i.e. for the :math:`g`-th pulse - :math:`t_g - t_{g-1}`. - eigvals : array_like, shape (n_dt, d) - Eigenvalue vectors for each time pulse segment *g* with the first - axis counting the pulse segment, i.e. - ``HD == array([D_0, D_1, ...])``. - eigvecs : array_like, shape (n_dt, d, d) - Eigenvector matrices for each time pulse segment *g* with the first - axis counting the pulse segment, i.e. - ``HV == array([V_0, V_1, ...])``. - basis : Basis, shape (d**2, d, d) - The basis elements, in which the pulse control matrix will be expanded. - n_opers : array_like, shape (n_nops, d, d) - Noise operators :math:`B_\alpha`. - 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. - transf_control_operators : array_like, shape (n_dt, n_ctrl, d, d) - The control operators transformed into the eigenspace of the control - Hamiltonian. The drift operators are ignored, if identifiers for - accessible control operators are provided. - s_derivs : array_like, shape (n_nops, n_ctrl, n_dt) - The derivatives of the noise susceptibilities by the control amplitudes. - Defaults to None. - - Returns - ------- - ctrlmat_data : dict {'R_g': R_g, 'dR_g': gradient} - * **R_g** *(array_like, shape (n_dt, n_nops, d**2, n_omega))* - The control matrix at each time step - :math:`\mathcal{R}_{\alpha j}^{(g)}(\omega)` is identified with R_g. - The array's indexing has shape :math:`(g,\alpha,j,\omega)`. - - * **dR_g** *(array_like, shape (n_dt, n_nops, d**2, n_ctrl, n_omega))* - The corresponding derivative with respect to the control strength - :math:`\frac{\partial\mathcal{R}_{\alpha j}^{(g)}(\omega)} - {\partial u_h(t_{g^\prime})}` is identified with dR_g. The array's - indexing has shape :math:`(g^\prime,\alpha,j,h,\omega)`. Here, only one - time axis is needed, since the derivative is zero for - :math:`g\neq g^\prime`. - - - Notes - ----- - The control matrix at each time step is evaluated according to - - .. math:: - - \mathcal{R}_{\alpha j}^{(g)}(\omega) = s_\alpha^{(g)}\mathrm{tr} - \left([\bar{B}_\alpha \circ I_1^{(g)}(\omega)] \bar{C}_j \right), - - where - - .. math:: - - I_{1,nm}^{(g)}(\omega) = \frac{\exp(\mathrm{i}(\omega + \omega_n^{(g)} - - \omega_m^{(g)}) \Delta t_g) - 1} - {\mathrm{i}(\omega + \omega_n^{(g)} - \omega_m^{(g)})} - - The derivative of the control matrix with respect to the control strength - at different time steps is calculated according to - - .. math:: - - \frac{\partial \mathcal{R}_{\alpha j}^{(g)}(\omega)} - {\partial u_h(t_{g^\prime})} = - \mathrm{i}\delta_{gg^\prime} s_\alpha^{(g)} \mathrm{tr} - \left( \bar{B}_{\alpha} \cdot \mathbb{M} \right) - + \frac{\partial s_\alpha^{(g)}}{u_h (t_{g^\prime})} \text{tr} - \left( (\overline{B}_{\alpha} \circ I_1^{(g)}{}(\omega)) \cdot - \overline{C}_{j}) \right). - - We assume that the noise susceptibility :math:`s` only depends locally - on the time i.e. :math:`\partial_{u(t_g)} s(t_{g^\prime}) - = \delta_{gg^\prime} \partial_{u(t_g)} s(t_g)` - If denoting :math:`\Delta\omega_{ij} = \omega_i^{(g)} - \omega_j^{(g)}` - the integral part is encapsulated in - - .. math:: - - \mathbb{M}_{mn} = \left[ \bar{C}_j, \mathbb{I}^{(mn)} - \circ \bar{H}_h \right]_{mn}, - - with - - .. math:: - - \mathbb{I}_{pq}^{(mn)} &= \delta_{pq} \left( - \frac{\Delta t_g \cdot - \exp(\mathrm{i}(\omega + \Delta\omega_{nm})\Delta t_g)} - {\mathrm{i}(\omega + \Delta\omega_{nm})} - + \frac{\exp(\mathrm{i}(\omega + \Delta\omega_{nm})\Delta t_g) - 1} - {(\omega + \Delta\omega_{nm})^2}\right)\\ - &+ \frac{1-\delta_{pq}}{\mathrm{i}\Delta\omega_{pq}} \cdot - \frac{\exp(\mathrm{i}(\omega + \Delta\omega_{nm} - + \Delta\omega_{pq})\Delta t_g) - 1} - {\mathrm{i}(\omega + \Delta\omega_{nm} + \Delta\omega_{pq})}\\ - &- \frac{1-\delta_{pq}}{\mathrm{i}\Delta\omega_{pq}} \cdot - \frac{\exp(\mathrm{i}(\omega + \Delta\omega_{nm})\Delta t_g) - 1} - {\mathrm{i}(\omega + \Delta\omega_{nm})} - """ - d = eigvecs.shape[-1] - n_dt = len(dt) - E = omega - - # Precompute some transformation into eigenspace of control Hamiltonian - path = ['einsum_path', (0, 1), (0, 1)] - VBV = np.einsum('gji,ajk,gkl->gail', eigvecs.conj(), n_opers, eigvecs, - optimize=path) - VCV = np.einsum('tnm,jnk,tkl->tjml', eigvecs.conj(), basis, eigvecs, - optimize=path) - - # Allocate memory - R_g = np.empty((n_dt, len(n_opers), len(basis), len(E)), dtype=complex) - R_g_d_s = np.empty( - (n_dt, len(n_opers), len(basis), len(transf_control_operators[0]), - len(E)), dtype=complex) - # For calculating dR_g: d,d = p,q, d,d = m,n - integral_deriv = np.empty((n_dt, len(E), d, d, d, d), dtype=complex) - - for g in range(n_dt): - # Any combination of \omega_m-\omega_n (R_g), \omega_p-\omega_q (dR_g) - dE = np.subtract.outer(eigvals[g], eigvals[g]) - # Any combination of \omega+\omega_m-\omega_n (R_g) or - # \omega-\omega_m+\omega_n (dR_g) - EdE = np.add.outer(E, dE) - - # 1) Calculation of the control matrix R_g at each time step - integral_Rg = np.empty((len(E), d, d), dtype=complex) - # Mask the integral to avoid convergence problems - mask_Rg = np.abs(EdE*dt[g]) <= 1e-7 - integral_Rg[mask_Rg] = dt[g] - integral_Rg[~mask_Rg] = (util.cexp(EdE[~mask_Rg]*dt[g]) - 1) \ - / (1j*(EdE[~mask_Rg])) - - R_g[g] = np.einsum('a,bcd,adc,hdc->abh', n_coeffs[:, g], VCV[g], - VBV[g], integral_Rg, - optimize=['einsum_path', (0, 2), (0, 2), (0, 1)]) - - # Add the dependency of the susceptibilities - # s_derivs has shape (n_nops, n_ctrl, n_dt) - # VCV has shape (n_dt, d**2, d, d) - # VBV has shape (n_dt, n_nops, d, d) - # integral_Rg has shape (n_omega, d, d) - # R_g_d_s shall be of shape (n_dt, n_nops, d**2, n_ctrl, n_omega) - if s_derivs is not None: - R_g_d_s[g] = np.einsum( - 'ae,bcd,adc,hdc->abeh', s_derivs[:, :, g], VCV[g], VBV[g], - integral_Rg, optimize=['einsum_path', (0, 2), (0, 2), (0, 1)]) - - # 2) Calculation of derivatives of control matrices at each time step - mask_deriv = (abs(dE) < 1e-15) # case: \omega_p-\omega_q = 0 - # calculation if omega_diff = 0 - n_case = sum(sum(mask_deriv)) - a = dt[g]*util.cexp(EdE*dt[g]) / (1j * EdE) \ - + (util.cexp(EdE*dt[g]) - 1) / (EdE)**2 - integral_deriv[g, :, mask_deriv] = np.concatenate(([[a]*n_case]), - axis=0) - # calculation if omega_diff != 0 - b1 = - (util.cexp(np.add.outer(EdE, dE[~mask_deriv])*dt[g]) - 1) \ - / (np.add.outer(EdE, dE[~mask_deriv])) / dE[~mask_deriv] - b2 = + np.divide.outer(((util.cexp(EdE*dt[g]) - 1) / EdE), dE[~mask_deriv]) - integral_deriv[g, :, ~mask_deriv] = (b1 + b2).transpose(3, 0, 1, 2) - - # Computation of the derivative/ gradient - I_circ_H = np.einsum('toijnm,thij->tohijnm', integral_deriv, - transf_control_operators) - M_mat = (np.einsum('tjmk,tohknnm->tojhmn', VCV, I_circ_H) - - np.einsum('tohmknm,tjkn->tojhmn', I_circ_H, VCV)) - gradient = 1j * np.einsum('at,tamn,tojhnm->tajho', n_coeffs, VBV, M_mat, - optimize=['einsum_path', (1, 2), (0, 1)]) - if s_derivs is not None: - gradient += R_g_d_s - ctrlmat_data = {'R_g': R_g, 'dR_g': gradient} - return ctrlmat_data - - -def control_matrix_at_timestep_derivative_refactored( - omega: Coefficients, - dt: Coefficients, - eigvals: ndarray, - eigvecs: ndarray, - basis: Basis, - c_opers_transformed: ndarray, - n_opers: Sequence[Operator], - n_coeffs: Sequence[Coefficients], - n_coeffs_deriv: Optional[Sequence[Coefficients]] = None, - intermediates: Optional[Dict[str, ndarray]] = None -) -> Tuple[ndarray, ndarray]: - """".""" - d = eigvecs.shape[-1] - n_dt = len(dt) - n_ctrl = len(c_opers_transformed[0]) - n_nops = len(n_opers) - E = omega - - # Precompute some transformation into eigenspace of control Hamiltonian - path = ['einsum_path', (0, 1), (0, 1)] - VBV = np.einsum('gji,ajk,gkl->gail', eigvecs.conj(), n_opers, eigvecs, optimize=path) - VCV = np.einsum('tnm,jnk,tkl->tjml', eigvecs.conj(), basis, eigvecs, optimize=path) - - # Allocate memory - R_g = np.empty((n_dt, n_nops, d**2, len(E)), dtype=complex) - buffer = np.empty((n_nops, d**2, len(E)), dtype=complex) - R_g_d_s = np.empty((n_dt, n_nops, d**2, n_ctrl, len(E)), dtype=complex) - # For calculating dR_g: d,d = p,q, d,d = m,n - integral_deriv = np.empty((n_dt, len(E), d, d, d, d), dtype=complex) - exp_buf, integral_Rg = np.empty((2, len(E), d, d), dtype=complex) - - expr = oe.contract_expression('bcd,adc,hdc->abh', VCV[0].shape, VBV[0].shape, - integral_Rg.shape, optimize=[(0, 1), (0, 1)]) - for g in range(n_dt): - integral_Rg = numeric._first_order_integral(E, eigvals[g], dt[g], exp_buf, integral_Rg) - - buffer = expr(VCV[g], VBV[g], integral_Rg, out=buffer) - R_g[g] = n_coeffs[:, g, None, None] * buffer - - # Add the dependency of the susceptibilities - # n_coeffs_deriv has shape (n_nops, n_ctrl, n_dt) - # VCV has shape (n_dt, d**2, d, d) - # VBV has shape (n_dt, n_nops, d, d) - # integral_Rg has shape (n_omega, d, d) - # R_g_d_s shall be of shape (n_dt, n_nops, d**2, n_ctrl, n_omega) - if n_coeffs_deriv is not None: - R_g_d_s[g] = n_coeffs_deriv[:, None, :, g, None] * buffer[:, :, None] - - integral_deriv[g] = _derivative_integral(E, eigvals[g], dt[g], integral_deriv[g]) - - # Computation of the derivative/ gradient - # I_circ_H = np.einsum('toijnm,thij->tohijnm', integral_deriv, c_opers_transformed) - I_circ_H = integral_deriv[:, :, None] * c_opers_transformed[:, None, ..., None, None] - if d == 2: - # Life a bit simpler - mask = np.eye(d, dtype=bool) - M_mat = np.zeros((n_dt, len(E), d**2, n_ctrl, d, d), dtype=complex) - M_mat[:, :, 1:] = np.einsum('tjmk,tohknnm->tojhmn', VCV[:, 1:], I_circ_H) - M_mat[:, :, 1:, :, mask] -= M_mat[:, :, 1:, :, mask][..., ::-1] - M_mat[:, :, 1:, :, ~mask] *= 2 - # M_mat -= M_mat.conj().swapaxes(-1, -2) - else: - M_mat = np.einsum('tjmk,tohknnm->tojhmn', VCV, I_circ_H) - M_mat -= np.einsum('tjkn,tohmknm->tojhmn', VCV, I_circ_H) - - ctrlmat_g_deriv = 1j*np.einsum('tamn,tojhnm->tajho', n_coeffs.T[..., None, None] * VBV, M_mat) - - if n_coeffs_deriv is not None: - ctrlmat_g_deriv += R_g_d_s - - return R_g, ctrlmat_g_deriv - - -def control_matrix_at_timestep_derivative_faster( - omega: Coefficients, - dt: Coefficients, - eigvals: ndarray, - eigvecs: ndarray, - basis: Basis, - c_opers_transformed: ndarray, - n_opers: Sequence[Operator], - n_coeffs: Sequence[Coefficients], - n_coeffs_deriv: Optional[Sequence[Coefficients]] = None -) -> Tuple[ndarray, ndarray]: - """".""" - d = eigvecs.shape[-1] - d2 = d**2 - n_dt = len(dt) - n_ctrl = len(c_opers_transformed[0]) - n_nops = len(n_opers) - n_omega = len(omega) - E = omega - - # Precompute some transformation into eigenspace of control Hamiltonian - n_opers_transformed = numeric._transform_hamiltonian(eigvecs, n_opers, n_coeffs).swapaxes(0, 1) - basis_transformed = numeric._transform_by_unitary(eigvecs[:, None], basis[None], - out=np.empty((n_dt, d2, d, d), complex)) - - # Allocate memory - ctrlmat_g = np.empty((n_dt, n_nops, d2, len(E)), dtype=complex) - ctrlmat_g_deriv_s = np.empty((n_dt, n_nops, d2, n_ctrl, len(E)), dtype=complex) - # For calculating dR_g: d,d = p,q, d,d = m,n - deriv_integral = np.empty((n_dt, len(E), d, d, d, d), dtype=complex) - exp_buf, integral = np.empty((2, len(E), d, d), dtype=complex) - - expr = oe.contract_expression('bcd,adc,hdc->abh', basis.shape, basis.shape, - integral.shape, optimize=[(0, 1), (0, 1)]) - for g in range(n_dt): - integral = numeric._first_order_integral(E, eigvals[g], dt[g], exp_buf, integral) - deriv_integral[g] = _derivative_integral(E, eigvals[g], dt[g], deriv_integral[g]) - - ctrlmat_g[g] = expr(basis_transformed[g], n_opers_transformed[g], integral, - out=ctrlmat_g[g]) - if n_coeffs_deriv is not None: - # equivalent contraction: 'ah,a,ako->akho', but this faster - ctrlmat_g_deriv_s[g] = ( - (n_coeffs_deriv[..., g] / n_coeffs[:, g, None])[:, None, :, None] - * ctrlmat_g[g, :, :, None] - ) - - # Basically a tensor product - # K = np.einsum('taij,thkl->tahikjl', VBV, VHV) - # L = K.transpose(0, 1, 2, 4, 3, 6, 5) - K = util.tensor(n_opers_transformed[:, :, None], c_opers_transformed[:, None]) - L = util.tensor_transpose(K, (1, 0), [[d, d], [d, d]]) - k = np.diagonal(K.reshape(n_dt, n_nops, n_ctrl, d, d, d, d), 0, -2, -3) - l = np.diagonal(L.reshape(n_dt, n_nops, n_ctrl, d, d, d, d), 0, -2, -3) - i1 = np.diagonal(deriv_integral, 0, -2, -3) - i2 = np.diagonal(deriv_integral, 0, -1, -4) - # Let magic ensue. Reshaping in F-major is faster than not (~ factor 2-4) - M = np.einsum( - 'tahpm,topm->tahop', - l.reshape(n_dt, n_nops, n_ctrl, d2, d, order='F'), - i1.reshape(n_dt, n_omega, d2, d, order='F') - ).reshape(n_dt, n_nops, n_ctrl, n_omega, d, d, order='F') - if d == 2: - # Life a bit simpler - mask = np.eye(d, dtype=bool) - M[..., mask] -= M[..., mask][..., ::-1] - M[..., ~mask] *= 2 - else: - M -= np.einsum( - 'tahpn,topn->tahop', - k.swapaxes(-2, -3).reshape(n_dt, n_nops, n_ctrl, d2, d, order='F'), - i2.reshape(n_dt, n_omega, d2, d, order='F') - ).reshape(n_dt, n_nops, n_ctrl, n_omega, d, d, order='F').swapaxes(-1, -2) - - # Expand in basis transformed to eigenspace - ctrlmat_g_deriv = np.einsum('tjnk,tahokn->tajho', 1j*basis_transformed, M) - - if n_coeffs_deriv is not None: - ctrlmat_g_deriv += ctrlmat_g_deriv_s - - return ctrlmat_g, ctrlmat_g_deriv - - -def _control_matrix_at_timestep_derivative_loop_1( - omega: Coefficients, - dt: Coefficients, - eigvals: ndarray, - eigvecs: ndarray, - basis_transformed, - c_opers_transformed: ndarray, - n_opers_transformed, - n_coeffs: Sequence[Coefficients], - n_coeffs_deriv: Sequence[Coefficients], - ctrlmat_step, - deriv_integral, - exp_buf, - integral, - expr -) -> Tuple[ndarray, ndarray]: - """".""" - d = len(eigvecs) - d2 = d**2 - n_ctrl = len(c_opers_transformed) - n_nops = len(n_opers_transformed) - n_omega = len(omega) - - integral = numeric._first_order_integral(omega, eigvals, dt, exp_buf, integral) - deriv_integral = _derivative_integral(omega, eigvals, dt, deriv_integral) - ctrlmat_step = expr(basis_transformed, n_opers_transformed, integral, out=ctrlmat_step) - if n_coeffs_deriv is not None: - # equivalent contraction: 'ah,a,ako->akho', but this faster - ctrlmat_step_deriv_s = ((n_coeffs_deriv / n_coeffs[:, None])[:, None, :, None] - * ctrlmat_step[:, :, None]) - - # Basically a tensor product - # K = np.einsum('taij,thkl->tahikjl', VBV, VHV) - # L = K.transpose(0, 1, 2, 4, 3, 6, 5) - K = util.tensor(n_opers_transformed[:, None], c_opers_transformed[None]) - L = util.tensor_transpose(K, (1, 0), [[d, d], [d, d]]) - k = np.diagonal(K.reshape(n_nops, n_ctrl, d, d, d, d), 0, -2, -3) - l = np.diagonal(L.reshape(n_nops, n_ctrl, d, d, d, d), 0, -2, -3) - i1 = np.diagonal(deriv_integral, 0, -2, -3) - i2 = np.diagonal(deriv_integral, 0, -1, -4) - # Let magic ensue. Reshaping in F-major is faster than not (~ factor 2-4) - M = np.einsum( - 'ahpm,opm->ahop', - l.reshape(n_nops, n_ctrl, d2, d, order='F'), - i1.reshape(n_omega, d2, d, order='F') - ).reshape(n_nops, n_ctrl, n_omega, d, d, order='F') - if d == 2: - # Life a bit simpler - mask = np.eye(d, dtype=bool) - M[..., mask] -= M[..., mask][..., ::-1] - M[..., ~mask] *= 2 - else: - M -= np.einsum( - 'ahpn,opn->ahop', - k.swapaxes(-2, -3).reshape(n_nops, n_ctrl, d2, d, order='F'), - i2.reshape(n_omega, d2, d, order='F') - ).reshape(n_nops, n_ctrl, n_omega, d, d, order='F').swapaxes(-1, -2) - - # Expand in basis transformed to eigenspace - ctrlmat_step_deriv = np.einsum('jnk,ahokn->ajho', 1j*basis_transformed, M) - - if n_coeffs_deriv is not None: - ctrlmat_step_deriv += ctrlmat_step_deriv_s - - return ctrlmat_step, ctrlmat_step_deriv - - -def _control_matrix_at_timestep_derivative_loop_2( - omega: Coefficients, - dt: Coefficients, - eigvals: ndarray, - eigvecs: ndarray, - basis_transformed, - c_opers_transformed: ndarray, - n_opers_transformed, - n_coeffs: Sequence[Coefficients], - n_coeffs_deriv: Sequence[Coefficients], - ctrlmat_step, - phase_factor, - deriv_integral, - exp_buf, - integral, - expr, - out -) -> Tuple[ndarray, ndarray]: - """".""" - d = len(eigvecs) - d2 = d**2 - n_ctrl = len(c_opers_transformed) - n_nops = len(n_opers_transformed) - n_omega = len(omega) - - deriv_integral = _derivative_integral(omega, eigvals, dt, out=deriv_integral) - - # Basically a tensor product - # K = np.einsum('taij,thkl->tahikjl', VBV, VHV) - # L = K.transpose(0, 1, 2, 4, 3, 6, 5) - K = util.tensor(n_opers_transformed[:, None], c_opers_transformed[None]) - L = util.tensor_transpose(K, (1, 0), [[d, d], [d, d]]) - k = np.diagonal(K.reshape(n_nops, n_ctrl, d, d, d, d), 0, -2, -3) - l = np.diagonal(L.reshape(n_nops, n_ctrl, d, d, d, d), 0, -2, -3) - i1 = np.diagonal(deriv_integral, 0, -2, -3) - i2 = np.diagonal(deriv_integral, 0, -1, -4) - # reshaping in F-major is significantly faster than C-major (~ factor 2-4) - M = np.einsum( - 'ahpm,opm->ahop', - l.reshape(n_nops, n_ctrl, d2, d, order='F'), - i1.reshape(n_omega, d2, d, order='F') - ).reshape(n_nops, n_ctrl, n_omega, d, d, order='F') - if d == 2: - # Life a bit simpler - mask = np.eye(d, dtype=bool) - M[..., mask] -= M[..., mask][..., ::-1] - M[..., ~mask] *= 2 - else: - M -= np.einsum( - 'ahpn,opn->ahop', - k.swapaxes(-2, -3).reshape(n_nops, n_ctrl, d2, d, order='F'), - i2.reshape(n_omega, d2, d, order='F') - ).reshape(n_nops, n_ctrl, n_omega, d, d, order='F').swapaxes(-1, -2) - - # Expand in basis transformed to eigenspace - ctrlmat_step_deriv = np.einsum('jnk,ahokn->ajho', 1j*basis_transformed, M) - ctrlmat_step_deriv *= phase_factor - - if n_coeffs_deriv is not None: - # equivalent contraction: 'ah,a,ako->akho', but this faster - ctrlmat_step_deriv += ((n_coeffs_deriv / n_coeffs[:, None])[:, None, :, None] - * ctrlmat_step[:, :, None]) - - return ctrlmat_step_deriv - - -def noise_operators_at_timestep_derivative( - omega: Coefficients, - dt: Coefficients, - eigvals: ndarray, - eigvecs: ndarray, - basis: Basis, - n_opers: Sequence[Operator], - n_coeffs: Sequence[Coefficients], - VHV: ndarray, - s_derivs: Optional[Sequence[Coefficients]] = None, - show_progressbar: Optional[bool] = None, - intermediates: Optional[Dict[str, ndarray]] = None -): - d = eigvecs.shape[-1] - n_dt = len(dt) - n_ctrl = len(VHV[0]) - n_nops = len(n_opers) - E = omega - - if intermediates is None: - noise_operators = np.empty((n_dt, len(E), n_nops, d, d), dtype=complex) - VBV = numeric._transform_hamiltonian(eigvecs, n_opers, n_coeffs).swapaxes(0, 1) - exp_buf, integral = np.empty((2, len(E), d, d), dtype=complex) - else: - noise_operators = intermediates['noise_operators_step'] - VBV = intermediates['n_opers_transformed'].swapaxes(0, 1) - - # R_g_d_s = np.empty((n_dt, len(E), n_nops, n_ctrl, d, d), dtype=complex) - # For calculating dR_g: d,d = p,q, d,d = m,n - integral_deriv = np.empty((n_dt, len(E), d, d, d, d), dtype=complex) - - expr_1 = oe.contract_expression('akl,okl->oakl', VBV[:, 0].shape, (len(E), d, d)) - expr_2 = oe.contract_expression('ij,...jk,lk', - eigvecs[0].shape, noise_operators[0].shape, eigvecs[0].shape, - optimize=[(0, 1), (0, 1)]) - for g in range(n_dt): - if intermediates is None: - integral = numeric._first_order_integral(E, eigvals[g], dt[g], exp_buf, integral) - noise_operators[g] = expr_1(VBV[g], integral, out=noise_operators[g]) - noise_operators[g] = expr_2(eigvecs[g], noise_operators[g], eigvecs[g].conj(), - out=noise_operators[g]) - - # Add the dependency of the susceptibilities - # s_derivs has shape (n_nops, n_ctrl, n_dt) - # VCV has shape (n_dt, d**2, d, d) - # VBV has shape (n_dt, n_nops, d, d) - # integral_Rg has shape (n_omega, d, d) - # R_g_d_s shall be of shape (n_dt, n_nops, d**2, n_ctrl, n_omega) - # if s_derivs is not None: - # Still need to think about n_coeffs here - # R_g_d_s[g] = s_derivs[:, None, :, g, None] * noise_operators[g] - - integral_deriv[g] = _derivative_integral(E, eigvals[g], dt[g], integral_deriv[g]) - - # (np.einsum('tcmn,tanp,topmn->tocamp', VHV, VBV, I1, optimize=['einsum_path', (0, 1), (0, 1)]) - - # np.einsum('tcnp,tamn,tompn->tocamp', VHV, VBV, I2, optimize=['einsum_path', (0, 1), (0, 1)])) - I1 = np.diagonal(integral_deriv, 0, -4, -1) - path = [(0, 1), (0, 1)] # optimization says [(0, 1), (0, 1)] - # M = oe.contract('tcmn,tanp,topmn->tocamp', VHV, VBV, I1, optimize=path) - M = oe.contract('tcmn,tanp,tompn->tocamp', VHV, VBV, I1, optimize=path) - if d == 2: - mask = np.eye(2, dtype=bool) - M[..., mask] -= M[..., mask][..., ::-1] - M[..., ~mask] *= 2 - else: - I2 = np.diagonal(integral_deriv, 0, -3, -2) - # M -= oe.contract('tcnp,tamn,tompn->tocamp', VHV, VBV, I2, optimize=path) - M -= oe.contract('tcnp,tamn,topmn->tocamp', VHV, VBV, I2, optimize=path) - - intermediates = intermediates or dict() - intermediates.setdefault('noise_operators_step', noise_operators) - intermediates.setdefault('n_opers_transformed', VBV.swapaxes(0, 1)) - - return 1j*M, intermediates - - -def calculate_derivative_of_control_matrix_from_scratch( - omega: Coefficients, - propagators: ndarray, - eigvals: ndarray, - eigvecs: ndarray, - basis: Basis, - t: Coefficients, - dt: Coefficients, - n_opers: Sequence[Operator], - n_coeffs: Sequence[Coefficients], - c_opers: Sequence[Operator], - all_identifiers: Sequence[str], - control_identifiers: Optional[Sequence[str]] = None, - s_derivs: Optional[Sequence[Coefficients]] = None) -> ndarray: - r""" - Calculate the derivative of the control matrix from scratch. - - Parameters - ---------- - omega : array_like, shape (n_omega,) - Frequencies, at which the pulse control matrix is to be evaluated. - propagators : array_like, shape (n_dt+1, d, d) - The propagators :math:`Q_g = P_g P_{g-1}\cdots P_0` as a (d, d) array - with *d* the dimension of the Hilbert space. - Q_Liou : ndarray, shape (n_dt+1, d**2, d**2) - The Liouville representation of the cumulative control propagators - U_c(t_g,0). - eigvals : array_like, shape (n_dt, d) - Eigenvalue vectors for each time pulse segment *g* with the first - axis counting the pulse segment, i.e. - ``HD == array([D_0, D_1, ...])``. - eigvecs : array_like, shape (n_dt, d, d) - Eigenvector matrices for each time pulse segment *g* with the first - axis counting the pulse segment, i.e. - ``HV == array([V_0, V_1, ...])``. - basis : Basis, shape (d**2, d, d) - The basis elements, in which the pulse control matrix will be expanded. - t : array_like, shape (n_dt+1), optional - The absolute times of the different segments. - dt : array_like, shape (n_dt) - Sequence duration, i.e. for the :math:`g`-th pulse - :math:`t_g - t_{g-1}`. - n_opers : array_like, shape (n_nops, d, d) - Noise operators :math:`B_\alpha`. - 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. - c_opers : array_like, shape (n_cops, d, d) - Control operators :math:`H_k`. - all_identifiers : array_like, shape (n_cops) - Identifiers of all control operators. - control_identifiers : Sequence[str], shape (n_ctrl), Optional - Sequence of strings with the control identifiers to distinguish between - accessible control and drift Hamiltonian. The default is None. - s_derivs : array_like, shape (n_nops, n_ctrl, n_dt) - The derivatives of the noise susceptibilities by the control amplitudes. - Defaults to None. - - Raises - ------ - ValueError - If the given identifiers *control_identifier* are not subset of the - total identifiers *all_identifiers* of all control operators. - - Returns - ------- - dR : array_like, shape (n_nops, d**2, n_dt, n_ctrl, n_omega) - Partial derivatives of the total control matrix with respect to each - control direction - :math:`\frac{\partial R_{\alpha k}(\omega)}{\partial u_h(t_{g'})}`. - The array's indexing has shape :math:`(\alpha,k,g^\prime,h,\omega)`. - - Notes - ----- - The derivative of the control matrix is calculated according to - - .. math :: - - \frac{\partial R_{\alpha k}(\omega)}{\partial u_h(t_{g'})} = - \sum_{g=1}^G \mathrm{e}^{\mathrm{i}\omega t_{g-1}}\cdot\left(\sum_j - \left[\frac{\partial R_{\alpha j}^{(g)}(\omega)} - {\partial u_h(t_{g'})} \cdot \mathcal{Q}_{jk}^{(g-1)} - + R_{\alpha j}^{(g)}(\omega) - \cdot\frac{\partial \mathcal{Q}_{jk}^{(g-1)}} - {\partial u_h(t_{g'})} \right] \right) - - See Also - -------- - :func:`liouville_derivative` - :func:`control_matrix_at_timestep_derivative` - """ - # Distinction between control and drift operators and only calculate the - # derivatives in control direction - - path = ['einsum_path', (0, 1), (0, 1)] - if (control_identifiers is None): - VHV = np.einsum('tji,hjk,tkl->thil', eigvecs.conj(), c_opers, eigvecs, - optimize=path) - elif (set(control_identifiers) <= set(all_identifiers)): - dict_id_oper = dict(zip(all_identifiers, c_opers)) - control = [dict_id_oper[element] for element in control_identifiers] - VHV = np.einsum('tji,hjk,tkl->thil', eigvecs.conj(), control, eigvecs, - optimize=path) - else: - raise ValueError('Given control identifiers have to be a \ - subset of (drift+control) Hamiltonian!') - - Q_Liou = superoperator.liouville_representation(propagators, basis) - # Get derivative of Liouville, control matrix at each time step, derivative - # derivative of control matrix at each time step - dL = liouville_derivative( - dt=dt, propagators=propagators, basis=basis, eigvecs=eigvecs, - eigvals=eigvals, transf_control_operators=VHV) - ctrlmat_data = control_matrix_at_timestep_derivative( - omega=omega, - dt=dt, - eigvals=eigvals, - eigvecs=eigvecs, - basis=basis, - n_opers=n_opers, - n_coeffs=n_coeffs, - transf_control_operators=VHV, - s_derivs=s_derivs - ) - ctrlmat_g = ctrlmat_data['R_g'] - ctrlmat_g_deriv = ctrlmat_data['dR_g'] - - # Calculate the derivative of the total control matrix - exp_factor = util.cexp(np.multiply.outer(t, omega)) - summand1 = np.einsum('to,tajho,tjk->aktho', exp_factor[:-1], - ctrlmat_g_deriv, Q_Liou[:-1], - optimize=['einsum_path', (1, 2), (0, 1)]) - summand2 = np.einsum('to,tajo,thsjk->aksho', exp_factor[1:-1], - ctrlmat_g[1:], dL[:-1], - optimize=['einsum_path', (0, 1), (0, 1)]) - - dR = summand1 + summand2 - return dR - - -def calculate_derivative_of_control_matrix_from_scratch_refactored( - omega: Coefficients, - propagators: ndarray, - eigvals: ndarray, - eigvecs: ndarray, - basis: Basis, - t: Coefficients, - dt: Coefficients, - n_opers: Sequence[Operator], - n_coeffs: Sequence[Coefficients], - c_opers: Sequence[Operator], - all_identifiers: Sequence[str], - control_identifiers: Optional[Sequence[str]] = None, - n_coeffs_deriv: Optional[Sequence[Coefficients]] = None) -> ndarray: - r""".""" - # Distinction between control and drift operators and only calculate the - # derivatives in control direction - - if control_identifiers is None: - control = c_opers - elif set(control_identifiers) <= set(all_identifiers): - dict_id_oper = dict(zip(all_identifiers, c_opers)) - control = [dict_id_oper[element] for element in control_identifiers] - else: - raise ValueError('Given control identifiers have to be a \ - subset of (drift+control) Hamiltonian!') - - c_opers_transformed = numeric._transform_hamiltonian(eigvecs, control).swapaxes(0, 1) - propagators_liouville = superoperator.liouville_representation(propagators[:-1], basis) - - propagators_liouville_deriv = liouville_derivative_vectorized(dt, propagators, basis, eigvecs, - eigvals, c_opers_transformed) - # ctrlmat_g, ctrlmat_g_deriv = control_matrix_at_timestep_derivative_refactored( - ctrlmat_g, ctrlmat_g_deriv = control_matrix_at_timestep_derivative_faster( - omega, dt, eigvals, eigvecs, basis, c_opers_transformed, n_opers, n_coeffs, n_coeffs_deriv - ) - - # Calculate the derivative of the total control matrix - exp_factor = util.cexp(np.multiply.outer(omega, t[:-1])) - # Equivalent (but slower) einsum contraction for first term: - # ctrlmat_deriv = np.einsum('ot,tajho,tjk->hotak', - # exp_factor, ctrlmat_g_deriv, propagators_liouville, - # optimize=['einsum_path', (0, 2), (0, 1)]) - ctrlmat_deriv = (ctrlmat_g_deriv.transpose(3, 4, 0, 1, 2) - @ (exp_factor[..., None, None] * propagators_liouville)) - ctrlmat_deriv += np.einsum('ot,tajo,thsjk->hosak', - exp_factor[:, 1:], ctrlmat_g[1:], propagators_liouville_deriv, - optimize=['einsum_path', (0, 1), (0, 1)]) - - return ctrlmat_deriv - - -def calculate_derivative_of_control_matrix_from_scratch_loop_1( - omega: Coefficients, - propagators: ndarray, - eigvals: ndarray, - eigvecs: ndarray, - basis: Basis, - t: Coefficients, - dt: Coefficients, - n_opers: Sequence[Operator], - n_coeffs: Sequence[Coefficients], - c_opers: Sequence[Operator], - intermediates, - all_identifiers: Sequence[str], - control_identifiers: Optional[Sequence[str]] = None, - n_coeffs_deriv: Optional[Sequence[Coefficients]] = None -) -> ndarray: - r""".""" - # Distinction between control and drift operators and only calculate the - # derivatives in control direction - if control_identifiers is None: - control = c_opers - elif set(control_identifiers) <= set(all_identifiers): - dict_id_oper = dict(zip(all_identifiers, c_opers)) - control = [dict_id_oper[element] for element in control_identifiers] - else: - raise ValueError('Given control identifiers have to be a \ - subset of (drift+control) Hamiltonian!') - - d = eigvecs.shape[-1] - d2 = d**2 - n_dt = len(dt) - n_ctrl = len(control) - n_nops = len(n_opers) - n_omega = len(omega) - - # Precompute some transformation into eigenspace of control Hamiltonian - c_opers_transformed = numeric._transform_hamiltonian(eigvecs, control).swapaxes(0, 1) - n_opers_transformed = numeric._transform_hamiltonian(eigvecs, n_opers, n_coeffs).swapaxes(0, 1) - basis_transformed = numeric._transform_by_unitary(eigvecs[:, None], basis[None], - out=np.empty((n_dt, d2, d, d), complex)) - propagators_liouville = superoperator.liouville_representation(propagators[:-1], basis) - - propagators_liouville_deriv = liouville_derivative_vectorized(dt, propagators, basis, eigvecs, - eigvals, c_opers_transformed) - exp_factor = util.cexp(np.multiply.outer(omega, t[:-1])) - - # Allocate memory - ctrlmat_deriv = np.empty((n_ctrl, n_omega, n_dt, n_nops, d2), dtype=complex) - ctrlmat_step = np.empty((n_dt, n_nops, d2, n_omega), dtype=complex) - deriv_integral = np.empty((n_omega, d, d, d, d), dtype=complex) - exp_buf, integral = np.empty((2, n_omega, d, d), dtype=complex) - ctrlmat_expr = oe.contract_expression('bcd,adc,hdc->abh', basis.shape, basis.shape, - integral.shape, optimize=[(0, 1), (0, 1)]) - for g in range(n_dt): - ctrlmat_step[g], ctrlmat_step_deriv = _control_matrix_at_timestep_derivative_loop_1( - omega, dt[g], eigvals[g], eigvecs[g], basis_transformed[g], c_opers_transformed[g], - n_opers_transformed[g], n_coeffs[:, g], n_coeffs_deriv[:, :, g], ctrlmat_step[g], - deriv_integral, exp_buf, integral, ctrlmat_expr - ) - - # Calculate the derivative of the total control matrix - ctrlmat_deriv[:, :, g] = ((exp_factor[:, g] * ctrlmat_step_deriv).transpose(2, 3, 0, 1) - @ propagators_liouville[g]) - - ctrlmat_deriv += oe.contract('ot,tajo,thsjk->hosak', - exp_factor[:, 1:], ctrlmat_step[1:], propagators_liouville_deriv, - optimize=[(0, 1), (0, 1)]) - - return ctrlmat_deriv - - -def calculate_derivative_of_control_matrix_from_scratch_loop_2( - omega: Coefficients, - propagators: ndarray, - eigvals: ndarray, - eigvecs: ndarray, - basis: Basis, - t: Coefficients, - dt: Coefficients, - n_opers: Sequence[Operator], - n_coeffs: Sequence[Coefficients], - c_opers: Sequence[Operator], - intermediates, - all_identifiers: Sequence[str], - control_identifiers: Optional[Sequence[str]] = None, - n_coeffs_deriv: Optional[Sequence[Coefficients]] = None -) -> ndarray: - r""".""" - # Distinction between control and drift operators and only calculate the - # derivatives in control direction - if control_identifiers is None: - control = c_opers - elif set(control_identifiers) <= set(all_identifiers): - dict_id_oper = dict(zip(all_identifiers, c_opers)) - control = [dict_id_oper[element] for element in control_identifiers] - else: - raise ValueError('Given control identifiers have to be a \ - subset of (drift+control) Hamiltonian!') - - d = eigvecs.shape[-1] - d2 = d**2 - n_dt = len(dt) - n_ctrl = len(control) - n_nops = len(n_opers) - n_omega = len(omega) - - # Precompute some transformation into eigenspace of control Hamiltonian - c_opers_transformed = numeric._transform_hamiltonian(eigvecs, control).swapaxes(0, 1) - n_opers_transformed = numeric._transform_hamiltonian(eigvecs, n_opers, n_coeffs).swapaxes(0, 1) - basis_transformed = numeric._transform_by_unitary(eigvecs[:, None], basis[None], - out=np.empty((n_dt, d2, d, d), complex)) - # U_deriv = liouville_derivative_vectorized_loop(dt, propagators, basis, eigvecs, eigvals, - # c_opers_transformed) - propagators_liouville = superoperator.liouville_representation(propagators[:-1], basis) - - # propagators_liouville_deriv = np.empty((n_ctrl, n_dt, d2, d2), dtype=complex) - # propagators_liouville_deriv = liouville_derivative_vectorized(dt, propagators, basis, eigvecs, - # eigvals, c_opers_transformed) - propagators_deriv = propagators_derivative_vectorized_completely(dt, propagators, basis, - eigvecs, eigvals, - c_opers_transformed) - liouville_deriv = np.empty((n_ctrl, n_dt, d2, d2), dtype=complex) - # Need to 'remove' the propagators from the control matrix at time step as computed by - # numeric.calculate_control_matrix_from_scratch - ctrlmat_step = (intermediates['control_matrix_step'].transpose(3, 0, 1, 2) - @ propagators_liouville.conj().swapaxes(-1, -2)).transpose(1, 2, 3, 0) - # exp_factor = util.cexp(np.multiply.outer(omega, t[:-1])) - - # Allocate memory - ctrlmat_deriv = np.empty((n_ctrl, n_omega, n_dt, n_nops, d2), dtype=complex) - deriv_integral = np.empty((n_omega, d, d, d, d), dtype=complex) - ctrlmat_step_deriv = np.empty((n_nops, d2, n_ctrl, n_omega), dtype=complex) - exp_buf, integral = np.empty((2, n_omega, d, d), dtype=complex) - ctrlmat_expr = oe.contract_expression('bcd,adc,hdc->abh', basis.shape, basis.shape, - integral.shape, optimize=[(0, 1), (0, 1)]) - liouville_deriv_expr = oe.contract_expression('hsba,jbc,cd,kda->hsjk', - propagators_deriv.conj()[:, 0].shape, - basis.shape, propagators[0].shape, basis.shape, - optimize=[(0, 2), (0, 2), (0, 1)]) - for g in range(n_dt): - phase_factor = util.cexp(omega * t[g]) - # Calculate the whole propagator derivative - # propagators_deriv[:, :g+1] = (propagators[g+1] - # @ propagators[1:g+2].conj().swapaxes(-1, -2) - # @ U_deriv[:, :g+1] - # @ propagators[:g+1]) - - # # Now calculate derivative of Liouville representation. Operation is faster using matmul, but - # # not nice to read. Full einsum contraction would be as follows (we insert new axes to induce - # # an outer product between basis elements steps using broadcasting): - # # liouville_deriv = np.einsum('hsba,jbc,cd,kda->hsjk', - # # deriv_U.conj(), basis, propagators[g+1], basis) - # propagators_liouville_deriv = np.einsum('hsba,jkba->hsjk', propagators_deriv.conj(), - # (basis @ propagators[g+1])[:, None] @ basis, - # out=propagators_liouville_deriv) - - ctrlmat_step_deriv = _control_matrix_at_timestep_derivative_loop_2( - omega, dt[g], eigvals[g], eigvecs[g], basis_transformed[g], c_opers_transformed[g], - n_opers_transformed[g], n_coeffs[:, g], n_coeffs_deriv[:, :, g], ctrlmat_step[g], - phase_factor, deriv_integral, exp_buf, integral, ctrlmat_expr, out=ctrlmat_step_deriv - ) - - # Calculate the derivative of the total control matrix (phase already included) - ctrlmat_deriv[:, :, g] = (ctrlmat_step_deriv.transpose(2, 3, 0, 1) - @ propagators_liouville[g]) - - if g < n_dt - 1: - liouville_deriv = liouville_derivative_vectorized_completely_loop( - propagators_deriv[:, g], basis, propagators[g+1], liouville_deriv_expr, - out=liouville_deriv - ) - ctrlmat_deriv += oe.contract('ajo,hsjk->hosak', - ctrlmat_step[g+1], liouville_deriv.real) - - # ctrlmat_deriv += oe.contract('ot,tajo,thsjk->hosak', - # exp_factor[:, 1:], ctrlmat_g[1:], propagators_liouville_deriv, - # optimize=[(0, 1), (0, 1)]) - # ctrlmat_deriv += oe.contract('otaj,thsjk->hosak', prefactor, propagators_liouville_deriv) - - return ctrlmat_deriv - - -def calculate_derivative_of_control_matrix_from_scratch_loop_3( - omega: Coefficients, - propagators: ndarray, - eigvals: ndarray, - eigvecs: ndarray, - basis: Basis, - t: Coefficients, - dt: Coefficients, - n_opers: Sequence[Operator], - n_coeffs: Sequence[Coefficients], - c_opers: Sequence[Operator], - intermediates, - all_identifiers: Sequence[str], - control_identifiers: Optional[Sequence[str]] = None, - n_coeffs_deriv: Optional[Sequence[Coefficients]] = None -) -> ndarray: - r""".""" - # Distinction between control and drift operators and only calculate the - # derivatives in control direction - if control_identifiers is None: - control = c_opers - elif set(control_identifiers) <= set(all_identifiers): - dict_id_oper = dict(zip(all_identifiers, c_opers)) - control = [dict_id_oper[element] for element in control_identifiers] - else: - raise ValueError('Given control identifiers have to be a \ - subset of (drift+control) Hamiltonian!') - - d = eigvecs.shape[-1] - d2 = d**2 - n_dt = len(dt) - n_ctrl = len(control) - n_nops = len(n_opers) - n_omega = len(omega) - - # Precompute some transformation into eigenspace of control Hamiltonian - c_opers_transformed = numeric._transform_hamiltonian(eigvecs, control).swapaxes(0, 1) - n_opers_transformed = numeric._transform_hamiltonian(eigvecs, n_opers, n_coeffs).swapaxes(0, 1) - basis_transformed = numeric._transform_by_unitary(eigvecs[:, None], basis[None], - out=np.empty((n_dt, d2, d, d), complex)) - # U_deriv = liouville_derivative_vectorized_loop(dt, propagators, basis, eigvecs, eigvals, - # c_opers_transformed) - propagators_liouville = superoperator.liouville_representation(propagators[:-1], basis) - - # propagators_liouville_deriv = np.empty((n_ctrl, n_dt, d2, d2), dtype=complex) - # propagators_liouville_deriv = liouville_derivative_vectorized(dt, propagators, basis, eigvecs, - # eigvals, c_opers_transformed) - propagators_deriv = propagators_derivative_vectorized_completely(dt, propagators, basis, - eigvecs, eigvals, - c_opers_transformed) - # Contraction 'jbc,tcd,kda->tjbka' - basis_propagators = np.tensordot(np.tensordot(propagators[1:-1], basis, axes=[1, 2]), - basis, axes=[1, 1]) - propagators_liouville_deriv = oe.contract('tjnkm,htsnm->tsjkh', - basis_propagators, propagators_deriv) - - # Need to 'remove' the propagators from the control matrix at time step as computed by - # numeric.calculate_control_matrix_from_scratch - ctrlmat_step = (intermediates['control_matrix_step'].transpose(3, 0, 1, 2) - @ propagators_liouville.conj().swapaxes(-1, -2)).transpose(1, 2, 3, 0) - # exp_factor = util.cexp(np.multiply.outer(omega, t[:-1])) - - # Allocate memory - ctrlmat_deriv = np.empty((n_ctrl, n_omega, n_dt, n_nops, d2), dtype=complex) - deriv_integral = np.empty((n_omega, d, d, d, d), dtype=complex) - ctrlmat_step_deriv = np.empty((n_nops, d2, n_ctrl, n_omega), dtype=complex) - exp_buf, integral = np.empty((2, n_omega, d, d), dtype=complex) - ctrlmat_expr = oe.contract_expression('bcd,adc,hdc->abh', basis.shape, basis.shape, - integral.shape, optimize=[(0, 1), (0, 1)]) - liouville_deriv_expr = oe.contract_expression('tajo,htnm,jnp,tpq,kqm->hoak', - ctrlmat_step[1:].shape, - propagators_deriv.conj()[:, :, 0].shape, - basis.shape, propagators[1:-1].shape, - basis.shape, - optimize=[(2, 3), (2, 3), (1, 2), (0, 1)]) - for g in range(n_dt): - phase_factor = util.cexp(omega * t[g]) - # Calculate the whole propagator derivative - # propagators_deriv[:, :g+1] = (propagators[g+1] - # @ propagators[1:g+2].conj().swapaxes(-1, -2) - # @ U_deriv[:, :g+1] - # @ propagators[:g+1]) - - # # Now calculate derivative of Liouville representation. Operation is faster using matmul, but - # # not nice to read. Full einsum contraction would be as follows (we insert new axes to induce - # # an outer product between basis elements steps using broadcasting): - # # liouville_deriv = np.einsum('hsba,jbc,cd,kda->hsjk', - # # deriv_U.conj(), basis, propagators[g+1], basis) - # propagators_liouville_deriv = np.einsum('hsba,jkba->hsjk', propagators_deriv.conj(), - # (basis @ propagators[g+1])[:, None] @ basis, - # out=propagators_liouville_deriv) - - ctrlmat_step_deriv = _control_matrix_at_timestep_derivative_loop_2( - omega, dt[g], eigvals[g], eigvecs[g], basis_transformed[g], c_opers_transformed[g], - n_opers_transformed[g], n_coeffs[:, g], n_coeffs_deriv[:, :, g], ctrlmat_step[g], - phase_factor, deriv_integral, exp_buf, integral, ctrlmat_expr, out=ctrlmat_step_deriv - ) - - # Calculate the derivative of the total control matrix (phase already included) - ctrlmat_deriv[:, :, g] = (ctrlmat_step_deriv.transpose(2, 3, 0, 1) - @ propagators_liouville[g]) - - propagators_liou_deriv = oe.contract('htba,jbc,tcd,kda->thjk', propagators_deriv[:, :, g], - basis, propagators[1:-1], basis) - ctrlmat_deriv[:, :, g] += oe.contract('tajo,thjk->hoak', ctrlmat_step[1:], - propagators_liou_deriv) - - # ctrlmat_deriv[:, :, g] += oe.contract('tajo,htnm,tjnkm->hoak', ctrlmat_step[1:], - # propagators_deriv.conj()[:, :, g], basis_propagators, - # optimize=[(1, 2), (0, 1)]) - - ctrlmat_deriv[:, :, g] += oe.contract('tajo,tjkh->hoak', - ctrlmat_step[1:], propagators_liouville_deriv[:, g]) - - # ctrlmat_deriv[:, :, g] += liouville_deriv_expr(ctrlmat_step[1:], - # propagators_deriv.conj()[:, :, 0], - # basis, propagators[1:-1], basis) - - # ctrlmat_deriv += oe.contract('ot,tajo,thsjk->hosak', - # exp_factor[:, 1:], ctrlmat_step[1:], propagators_liouville_deriv, - # optimize=[(0, 1), (0, 1)]) - # ctrlmat_deriv += oe.contract('otaj,thsjk->hosak', prefactor, propagators_liouville_deriv) - - return ctrlmat_deriv - - -def calculate_derivative_of_control_matrix_from_scratch_loop_4( - omega: Coefficients, - propagators: ndarray, - eigvals: ndarray, - eigvecs: ndarray, - basis: Basis, - t: Coefficients, - dt: Coefficients, - n_opers: Sequence[Operator], - n_coeffs: Sequence[Coefficients], - c_opers: Sequence[Operator], - intermediates, - all_identifiers: Sequence[str], - control_identifiers: Optional[Sequence[str]] = None, - n_coeffs_deriv: Optional[Sequence[Coefficients]] = None -) -> ndarray: - r""".""" - # Distinction between control and drift operators and only calculate the - # derivatives in control direction - if control_identifiers is None: - control = c_opers - elif set(control_identifiers) <= set(all_identifiers): - dict_id_oper = dict(zip(all_identifiers, c_opers)) - control = [dict_id_oper[element] for element in control_identifiers] - else: - raise ValueError('Given control identifiers have to be a \ - subset of (drift+control) Hamiltonian!') - - d = eigvecs.shape[-1] - d2 = d**2 - n_dt = len(dt) - n_ctrl = len(control) - n_nops = len(n_opers) - n_omega = len(omega) - - # Precompute some transformation into eigenspace of control Hamiltonian - c_opers_transformed = numeric._transform_hamiltonian(eigvecs, control).swapaxes(0, 1) - n_opers_transformed = numeric._transform_hamiltonian(eigvecs, n_opers, n_coeffs).swapaxes(0, 1) - basis_transformed = numeric._transform_by_unitary(eigvecs[:, None], basis[None], - out=np.empty((n_dt, d2, d, d), complex)) - propagators_liouville = superoperator.liouville_representation(propagators[:-1], basis) - - propagators_liouville_deriv = liouville_derivative_vectorized(dt, propagators, basis, eigvecs, - eigvals, c_opers_transformed) - ctrlmat_step = intermediates['control_matrix_step'] - prefactor = (ctrlmat_step[1:].transpose(-1, 0, 1, 2) - @ propagators_liouville[1:].conj().swapaxes(-1, -2)) - - # Allocate memory - ctrlmat_step_deriv = np.empty((n_nops, d2, n_ctrl, n_omega), dtype=complex) - ctrlmat_deriv = np.empty((n_ctrl, n_omega, n_dt, n_nops, d2), dtype=complex) - deriv_integral = np.empty((n_omega, d, d, d, d), dtype=complex) - exp_buf, integral = np.empty((2, n_omega, d, d), dtype=complex) - ctrlmat_expr = oe.contract_expression('bcd,adc,hdc->abh', basis.shape, basis.shape, - integral.shape, optimize=[(0, 1), (0, 1)]) - for g in range(n_dt): - ctrlmat_step_deriv = _control_matrix_at_timestep_derivative_loop_2( - omega, dt[g], eigvals[g], eigvecs[g], basis_transformed[g], c_opers_transformed[g], - n_opers_transformed[g], n_coeffs[:, g], n_coeffs_deriv[:, :, g], ctrlmat_step[g], - deriv_integral, exp_buf, integral, ctrlmat_expr, out=ctrlmat_step_deriv - ) - - # Calculate the derivative of the total control matrix - exp_factor = util.cexp(omega * t[g]) - ctrlmat_deriv[:, :, g] = ((exp_factor * ctrlmat_step_deriv).transpose(2, 3, 0, 1) - @ propagators_liouville[g]) - - ctrlmat_deriv += oe.contract('otaj,thsjk->hosak', prefactor, propagators_liouville_deriv) - - return ctrlmat_deriv - - -def calculate_derivative_of_noise_operators_from_scratch( - omega: Coefficients, - propagators: ndarray, - eigvals: ndarray, - eigvecs: ndarray, - basis: Basis, - t: Coefficients, - dt: Coefficients, - n_opers: Sequence[Operator], - n_coeffs: Sequence[Coefficients], - c_opers: Sequence[Operator], - all_identifiers: Sequence[str], - control_identifiers: Optional[Sequence[str]] = None, - s_derivs: Optional[Sequence[Coefficients]] = None, - intermediates: Optional[Dict[str, ndarray]] = None): - - if control_identifiers is None: - control = c_opers - elif set(control_identifiers) <= set(all_identifiers): - dict_id_oper = dict(zip(all_identifiers, c_opers)) - control = [dict_id_oper[element] for element in control_identifiers] - else: - raise ValueError('Given control identifiers have to be a \ - subset of (drift+control) Hamiltonian!') - - d = eigvecs.shape[-1] - VHV = numeric._transform_hamiltonian(eigvecs, control).swapaxes(0, 1) - VQ = numeric._propagate_eigenvectors(eigvecs, propagators[:-1]) - - exp_factor = util.cexp(np.multiply.outer(t[:-1], omega)) - M, intermediates = noise_operators_at_timestep_derivative(omega, dt, eigvals, eigvecs, basis, - n_opers, n_coeffs, VHV, - intermediates=intermediates) - - QVMVQ = numeric._transform_by_unitary(M, VQ[:, None, None, None], out=M) - QVMVQ *= exp_factor[..., None, None, None, None] - - noise_operators_step_pulse = intermediates['noise_operators_step_pulse'] - noise_operators_step_total = intermediates['noise_operators_step_total'] - - propagators_deriv = propagator_derivative(dt, propagators, eigvecs, eigvals, VHV) - propagators_deriv_factor = propagator_derivative_factor(dt, propagators, eigvecs, eigvals, VHV) - - noise_operators_deriv = np.zeros((len(omega), len(n_opers), d, len(c_opers), len(dt), d), - dtype=complex) - for g in range(len(dt)): - noise_operators_deriv[..., :g+1, :] += np.tensordot(noise_operators_step_total[g], - propagators_deriv_factor[:, g, :g+1], - axes=[-1, -2]) - noise_operators_deriv[..., :g+1, :] += np.tensordot(noise_operators_step_total[g], - propagators_deriv_factor[:, g, :g+1].conj(), - axes=[-2, -2]) - - noise_operators_deriv = np.zeros(QVMVQ.shape, dtype=complex) - noise_operators_deriv = np.einsum( - 'so,sji,soajk,hstkl->tohail', - exp_factor[:-1], propagators[1:-1].conj(), noise_operators_step_pulse[1:], propagators_deriv[:, :-1], - optimize=['einsum_path', (0, 1), (0, 2), (0, 1)] - ) - noise_operators_deriv += np.einsum( - 'so,hstji,soajk,skl->tohail', - exp_factor[:-1], propagators_deriv[:, :-1].conj(), noise_operators_step_pulse[1:], propagators[1:-1], - optimize=['einsum_path', (0, 1), (0, 2), (0, 1)] - ) - - noise_operators_deriv = np.zeros(QVMVQ.shape, dtype=complex) - noise_operators_deriv.real = 2 * np.einsum( - 'to,hstjk,skl,soali->tohaji', - exp_factor, propagators_deriv[:, :-1], propagators[1:-1], noise_operators_step_pulse[1:], - optimize=['einsum_path', (1, 2), (1, 2), (0, 1)] - ).real - noise_operators_deriv += QVMVQ - - return noise_operators_deriv - - -def calculate_filter_function_derivative(R: ndarray, deriv_R: ndarray) -> ndarray: - r""" - Compute the filter function derivative from the control matrix. - - Parameters - ---------- - R : array_like, shape (n_nops, d**2, n_omega) - The control matrix. - deriv_R: array_like, shape (n_nops, d**2, n_t, n_ctrl, n_omega) - The derivative of the control matrix. - - Returns - ------- - deriv_filter_function : ndarray, shape (n_nops, n_dt, n_ctrl, n_omega) - The regular filter functions' derivatives for variation in each control - contribution - :math:`\frac{\partial F_\alpha(\omega)}{\partial u_h(t_{g'})}`. - The array's indexing has shape :math:`(\alpha,g^\prime,h,\omega)`. - - Notes - ----- - The filter function derivative is calculated according to - - .. math :: - - \frac{\partial F_\alpha(\omega)}{\partial u_h(t_{g'})} - = 2 \mathrm{Re} \left( \sum_k R_{\alpha k}^\ast(\omega) - \frac{\partial R_{\alpha k}(\omega)} - {\partial u_h(t_{g'})} \right) - """ - summe = np.einsum('ako,aktho->atho', R.conj(), deriv_R) - return 2*summe.real - - -def calculate_filter_function_derivative_refactored(R: ndarray, R_deriv: ndarray) -> ndarray: - summe = np.einsum('ako,hotak->atho', R.conj(), R_deriv) - return 2*summe.real - - -def infidelity_derivative( - pulse: 'PulseSequence', - S: Union[Coefficients, Callable], - omega: Coefficients, - control_identifiers: Optional[Sequence[str]] = None, - s_derivs: Optional[Sequence[Coefficients]] = None -) -> ndarray: - r""" - Calculate the infidelity derivative. - - Calculate the entanglement infidelity derivative of the ``PulseSequence`` - *pulse*. - - Parameters - ---------- - pulse : PulseSequence - The ``PulseSequence`` instance, for which to calculate the infidelity - for. - S : array_like - The two-sided noise power spectral density in units of inverse - frequencies as an array of shape (n_omega,) or (n_nops, n_omega). In - the first case, the same spectrum is taken for all noise operators, in - the second, it is assumed that there are no correlations between - different noise sources and thus there is one spectrum for each noise - operator. - omega : array_like, shape (n_omega) - The frequencies at which the integration is to be carried out. - control_identifiers : Sequence[str], shape (n_ctrl) - Sequence of strings with the control identifiern to distinguish between - accessible control and drift Hamiltonian. - s_derivs : array_like, shape (n_nops, n_ctrl, n_dt) - The derivatives of the noise susceptibilities by the control amplitudes. - Defaults to None. - - Raises - ------ - ValueError - If the provided noise spectral density does not fit expected shape. - - Returns - ------- - deriv_infid : ndarray, shape (n_nops, n_dt, n_ctrl) - Array with the derivative of the infidelity for each noise source taken - for each control direction at each time step - :math:`\frac{\partial I_e}{\partial u_h(t_{g'})}`. The array's indexing - has shape :math:`(\alpha,g^\prime,h)`. - - Notes - ----- - The infidelity's derivative is given by - - .. math:: - - \frac{\partial I_e}{\partial u_h(t_{g'})} = \frac{1}{d} - \int_{-\infty}^\infty - \frac{d\omega}{2\pi} - S_\alpha(\omega) - \frac{\partial F_\alpha (\omega)} - {\partial u_h(t_{g'})} - - with :math:`S_{\alpha}(\omega)` the noise spectral density - and :math:`F_{\alpha}(\omega)` the canonical filter function for - noise source :math:`\alpha`. - - To convert to the average gate infidelity, use the - following relation given by Horodecki et al. [Hor99]_ and - Nielsen [Nie02]_: - - .. math:: - - \big\langle\mathcal{I}_\mathrm{avg}\big\rangle = \frac{d}{d+1} - \big\langle\mathcal{I}_\mathrm{e}\big\rangle. - - References - ---------- - .. [Hor99] - Horodecki, M., Horodecki, P., & Horodecki, R. (1999). General - teleportation channel, singlet fraction, and quasidistillation. - Physical Review A - Atomic, Molecular, and Optical Physics, 60(3), - 1888–1898. https://doi.org/10.1103/PhysRevA.60.1888 - - .. [Nie02] - Nielsen, M. A. (2002). A simple formula for the average gate - fidelity of a quantum dynamical operation. Physics Letters, - Section A: General, Atomic and Solid State Physics, 303(4), 249–252. - https://doi.org/10.1016/S0375-9601(02)01272-0 - """ - n_ops = len(pulse.n_opers) - S = np.asarray(S) - if(S.shape[0] == n_ops): - S_all = S - elif(S.shape[0] == 1): - S_all = np.array([S[0]]*n_ops) - elif(S.shape[0] == len(omega)): - S_all = np.array([S]*n_ops) - else: - raise ValueError('Not fitting shape of S.') - - filter_function_deriv = pulse.get_filter_function_derivative(omega, control_identifiers, - s_derivs) - - integrand = np.einsum('ao,atho->atho', S_all, filter_function_deriv) - infid_deriv = trapz(integrand, omega) / (2*np.pi*pulse.d) - - return infid_deriv diff --git a/pytest.ini b/pytest.ini index f9973a4..e859e2d 100644 --- a/pytest.ini +++ b/pytest.ini @@ -6,4 +6,8 @@ addopts = --verbose --cov=filter_functions --cov-config=.coveragerc - --cov-report=xml \ No newline at end of file + --cov-report=xml + --benchmark-autosave + --benchmark-columns='mean,stddev,min,max,rounds' + --benchmark-sort='mean' + From 35efbf588481d5dcc4f8ae25bb7fb085ddcf20e8 Mon Sep 17 00:00:00 2001 From: Tobias Hangleiter Date: Tue, 2 Mar 2021 16:51:22 +0100 Subject: [PATCH 08/11] Remove pytest-benchmark config from pytest.ini --- pytest.ini | 4 ---- 1 file changed, 4 deletions(-) diff --git a/pytest.ini b/pytest.ini index e859e2d..de2e972 100644 --- a/pytest.ini +++ b/pytest.ini @@ -7,7 +7,3 @@ addopts = --cov=filter_functions --cov-config=.coveragerc --cov-report=xml - --benchmark-autosave - --benchmark-columns='mean,stddev,min,max,rounds' - --benchmark-sort='mean' - From a4d39eb7bd9c921acc1f51d88fdcbe51e6b61a11 Mon Sep 17 00:00:00 2001 From: Tobias Hangleiter Date: Tue, 2 Mar 2021 17:26:45 +0100 Subject: [PATCH 09/11] Factor out oper parsing, add qopt.DenseOperator detection Used in both Basis and PulseSequence constructors --- filter_functions/basis.py | 12 +------- filter_functions/pulse_sequence.py | 24 ++------------- filter_functions/util.py | 49 ++++++++++++++++++++++++++++++ 3 files changed, 52 insertions(+), 33 deletions(-) diff --git a/filter_functions/basis.py b/filter_functions/basis.py index 7f37877..08f2da9 100644 --- a/filter_functions/basis.py +++ b/filter_functions/basis.py @@ -173,21 +173,11 @@ def __new__(cls, basis_array: Sequence, traceless: Optional[bool] = None, except AttributeError: pass - basis = np.empty((len(basis_array), *basis_array[0].shape), dtype=complex) + basis = util.parse_operators(basis_array, 'basis_array') if basis.shape[0] > np.product(basis.shape[1:]): raise ValueError('Given overcomplete set of basis matrices. ' 'Not linearly independent.') - for i, elem in enumerate(basis_array): - if isinstance(elem, ndarray): # numpy array - basis[i] = elem - elif hasattr(elem, 'full'): # qutip.Qobj - basis[i] = elem.full() - elif hasattr(elem, 'todense'): # sparse array - basis[i] = elem.todense() - else: - raise TypeError('At least one element invalid type!') - basis = basis.view(cls) basis.btype = btype or 'Custom' basis.d = basis.shape[-1] diff --git a/filter_functions/pulse_sequence.py b/filter_functions/pulse_sequence.py index 2c4696f..e48d065 100644 --- a/filter_functions/pulse_sequence.py +++ b/filter_functions/pulse_sequence.py @@ -1134,28 +1134,8 @@ def _parse_Hamiltonian(H: Hamiltonian, n_dt: int, H_str: str) -> Tuple[Sequence[ coeffs = args[0] identifiers = list(args[1]) - # Parse opers and convert to squeezed ndarray if possible - parsed_opers = [] - for oper in opers: - if isinstance(oper, ndarray): - parsed_opers.append(oper.squeeze()) - elif hasattr(oper, 'full'): - # qutip.Qobj - parsed_opers.append(oper.full()) - elif hasattr(oper, 'todense'): - # sparse object - parsed_opers.append(oper.todense()) - else: - raise TypeError(f'Expected operators in {H_str} to be NumPy arrays or QuTiP Qobjs!') - - # Check correct dimensions for the operators - if set(oper.ndim for oper in parsed_opers) != {2}: - raise ValueError(f'Expected all operators in {H_str} to be two-dimensional!') - - if len(set(*set(oper.shape for oper in parsed_opers))) != 1: - raise ValueError(f'Expected operators in {H_str} to be square!') - - parsed_opers = np.asarray(parsed_opers) + # Parse opers and convert to ndarray + parsed_opers = util.parse_operators(opers, H_str) if not all(hasattr(coeff, '__len__') for coeff in coeffs): raise TypeError(f'Expected coefficients in {H_str} to be a sequence') diff --git a/filter_functions/util.py b/filter_functions/util.py index 3f0bb26..3dab060 100644 --- a/filter_functions/util.py +++ b/filter_functions/util.py @@ -183,6 +183,55 @@ def wrapper(*args, **kwargs): parse_which_FF_parameter = parse_optional_parameters({'which': ('fidelity', 'generalized')}) +def parse_operators(opers: Sequence[Operator], err_loc: str) -> List[ndarray]: + """Parse a sequence of operators and convert to ndarray. + + Parameters + ---------- + opers: Sequence[Operator] + Sequence of operators. + err_loc: str + Some cosmetics for the exceptions to be raised. + + Raises + ------ + TypeError + If any operator is not a valid type. + ValueError + If not all operators are 2d and square. + + Returns + ------- + parse_opers: ndarray, shape (len(opers), *opers[0].shape) + The parsed ndarray. + + """ + parsed_opers = [] + for oper in opers: + if isinstance(oper, ndarray): + parsed_opers.append(oper.squeeze()) + elif hasattr(oper, 'full'): + # qutip.Qobj + parsed_opers.append(oper.full()) + elif hasattr(oper, 'todense'): + # sparse object + parsed_opers.append(oper.todense()) + elif hasattr(oper, 'data') and hasattr(oper, 'dexp'): + # qopt DenseMatrix + parsed_opers.append(oper.data) + else: + raise TypeError(f'Expected operators in {err_loc} to be NumPy arrays or QuTiP Qobjs!') + + # Check correct dimensions of the operators + if set(oper.ndim for oper in parsed_opers) != {2}: + raise ValueError(f'Expected all operators in {err_loc} to be two-dimensional!') + + if len(set(*set(oper.shape for oper in parsed_opers))) != 1: + raise ValueError(f'Expected operators in {err_loc} to be square!') + + return np.asarray(parsed_opers) + + def _tensor_product_shape(shape_A: Sequence[int], shape_B: Sequence[int], rank: int): """Get shape of the tensor product between A and B of rank rank""" broadcast_shape = () From 2718040b9e5ec6c5886f7cddb28b82ad61cb3924 Mon Sep 17 00:00:00 2001 From: Tobias Hangleiter Date: Tue, 2 Mar 2021 17:28:38 +0100 Subject: [PATCH 10/11] Be more lenient in random gradient test --- tests/test_gradient.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_gradient.py b/tests/test_gradient.py index f5418ab..ea3c182 100644 --- a/tests/test_gradient.py +++ b/tests/test_gradient.py @@ -70,7 +70,7 @@ def test_gradient_calculation_random_pulse(self): spectral_noise_density=gradient_testutil.one_over_f_noise, c_id=[f'c{i}' for i in range(len(u_ctrl))], n_coeffs_deriv=None ) - self.assertArrayAlmostEqual(ana_grad, fin_diff_grad, rtol=1e-6, atol=1e-10) + self.assertArrayAlmostEqual(ana_grad, fin_diff_grad, rtol=1e-6, atol=1e-8) def test_raises(self): pulse = testutil.rand_pulse_sequence(2, 3) From 29f637d493ba038cfb4101bd5880474283bb08b9 Mon Sep 17 00:00:00 2001 From: Tobias Hangleiter Date: Tue, 2 Mar 2021 19:47:11 +0100 Subject: [PATCH 11/11] Also compare modulo 2pi in Qutip test (cf 42dbb98) --- tests/test_util.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/test_util.py b/tests/test_util.py index c23cb1b..f5003d6 100644 --- a/tests/test_util.py +++ b/tests/test_util.py @@ -615,16 +615,16 @@ def test_oper_equiv(self): result = util.oper_equiv(psi, psi*np.exp(1j*phase)) self.assertTrue(result[0]) - self.assertAlmostEqual(result[1], phase, places=5) + self.assertAlmostEqual(np.mod(result[1], 2*np.pi), np.mod(phase, 2*np.pi), places=5) result = util.oper_equiv(psi*np.exp(1j*phase), psi) self.assertTrue(result[0]) - self.assertAlmostEqual(result[1], -phase, places=5) + self.assertAlmostEqual(np.mod(result[1], 2*np.pi), np.mod(-phase, 2*np.pi), places=5) result = util.oper_equiv(U, U*np.exp(1j*phase)) self.assertTrue(result[0]) - self.assertAlmostEqual(result[1], phase, places=5) + self.assertAlmostEqual(np.mod(result[1], 2*np.pi), np.mod(phase, 2*np.pi), places=5) result = util.oper_equiv(U*np.exp(1j*phase), U) self.assertTrue(result[0]) - self.assertAlmostEqual(result[1], -phase, places=5) + self.assertAlmostEqual(np.mod(result[1], 2*np.pi), np.mod(-phase, 2*np.pi), places=5)