n_coeffs_deriv sorting (#83)

Sort n_coeffs_deriv by identifiers

Author: thangleiter

filter_functions/gradient.py

@@ -262,7 +262,7 @@ def _control_matrix_at_timestep_derivative(
         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)}
+        strength :math:`\frac{\partial\mathcal{B}_{\alpha j}^{(g)}(\omega)}`
 
     Notes
     -----
@@ -384,8 +384,6 @@ def calculate_derivative_of_control_matrix_from_scratch(
         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: Optional[Dict[str, ndarray]] = None
 ) -> ndarray:
@@ -421,28 +419,16 @@ def calculate_derivative_of_control_matrix_from_scratch(
     n_coeffs: array_like, shape (n_nops, n_dt)
         The sensitivities of the system to the noise operators given by
         *n_opers* at the given time step.
-    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.
+    c_opers: array_like, shape (n_ctrl, d, d)
+        Control operators :math:`H_k` with respect to which the
+        derivative is computed.
     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.
-
     Returns
     -------
     ctrlmat_deriv: ndarray, shape (n_ctrl, n_omega, n_dt, n_nops, d**2)
@@ -470,24 +456,16 @@ def calculate_derivative_of_control_matrix_from_scratch(
     _liouville_derivative
     _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
-
     d = eigvecs.shape[-1]
     n_dt = len(dt)
-    n_ctrl = len(idx)
+    n_ctrl = len(c_opers)
     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)
+    c_opers_transformed = numeric._transform_hamiltonian(eigvecs, c_opers).swapaxes(0, 1)
     if not intermediates:
         # None or empty
         n_opers_transformed = numeric._transform_hamiltonian(eigvecs, n_opers,
@@ -575,6 +553,7 @@ def infidelity_derivative(
         spectrum: Coefficients,
         omega: Coefficients,
         control_identifiers: Optional[Sequence[str]] = None,
+        n_oper_identifiers: Optional[Sequence[str]] = None,
         n_coeffs_deriv: Optional[Sequence[Coefficients]] = None
 ) -> ndarray:
     r"""Calculate the entanglement infidelity derivative of the
@@ -599,11 +578,29 @@ def infidelity_derivative(
     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.
+        Sequence of strings with the control identifiers to
+        distinguish between control and drift Hamiltonian. The
+        default is None, in which case the derivative is computed
+        for all known non-noise operators.
+    n_oper_identifiers: Sequence[str], shape (n_nops,)
+        Sequence of strings with the noise identifiers for which to
+        compute the derivative contribution. The default is None, in
+        which case it is computed for all known noise operators.
     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.
+        amplitudes. The rows and columns should be in the same order
+        as the corresponding identifiers above. Defaults to None, in
+        which case the coefficients are assumed to be constant and
+        hence their derivative vanishing.
+
+        .. warning::
+
+            Internally, control and noise terms of the Hamiltonian
+            are stored alphanumerically sorted by their identifiers.
+            If the noise and/or control identifiers above are not
+            explicitly given, the rows and/or columns of this
+            parameter need to be sorted in the same fashion.
+
 
     Raises
     ------
@@ -616,7 +613,9 @@ def infidelity_derivative(
     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'})}`.
+        :math:`\frac{\partial I_e}{\partial u_h(t_{g'})}`. Sorted in
+        the same fashion as `n_coeffs_deriv` or, if not given,
+        alphanumerically by the identifiers.
 
     Notes
     -----
@@ -658,7 +657,9 @@ def infidelity_derivative(
         https://doi.org/10.1016/S0375-9601(02)01272-0
     """
     spectrum = numeric._parse_spectrum(spectrum, omega, range(len(pulse.n_opers)))
-    filter_function_deriv = pulse.get_filter_function_derivative(omega, control_identifiers,
+    filter_function_deriv = pulse.get_filter_function_derivative(omega,
+                                                                 control_identifiers,
+                                                                 n_oper_identifiers,
                                                                  n_coeffs_deriv)
 
     integrand = np.einsum('ao,atho->atho', spectrum, filter_function_deriv)

filter_functions/pulse_sequence.py

@@ -160,17 +160,23 @@ class PulseSequence:
     Attributes
     ----------
     c_opers: ndarray, shape (n_cops, d, d)
-        Control operators
+        Control operators. Note that they are stored sorted by their
+        corresponding identifiers.
     n_opers: ndarray, shape (n_nops, d, d)
-        Noise operators
+        Noise operators. Note that they are stored sorted by their
+        corresponding identifiers.
     c_oper_identifers: sequence of str
-        Identifiers for the control operators of the system
+        Identifiers for the control operators of the system. Stored
+        sorted.
     n_oper_identifers: sequence of str
-        Identifiers for the noise operators of the system
+        Identifiers for the noise operators of the system. Stored
+        sorted.
     c_coeffs: ndarray, shape (n_cops, n_dt)
-        Control parameters in units of :math:`\hbar`
+        Control parameters in units of :math:`\hbar`. Note that they
+        are stored sorted by their corresponding identifiers.
     n_coeffs: ndarray, shape (n_nops, n_dt)
-        Noise sensitivities in units of :math:`\hbar`
+        Noise sensitivities in units of :math:`\hbar`. Note that they
+        are stored sorted by their corresponding identifiers.
     dt: ndarray, shape (n_dt,)
         Time steps
     t: ndarray, shape (n_dt + 1,)
@@ -247,8 +253,10 @@ def __init__(self, *args, **kwargs) -> None:
         self.basis = None
 
         # Parse the input arguments and set attributes
-        attributes = ('c_opers', 'c_oper_identifiers', 'c_coeffs', 'n_opers',
-                      'n_oper_identifiers', 'n_coeffs', 'dt', 'd', 'basis')
+        # TODO: Jesus, this is in need of refactoring.
+        attributes = ('c_opers', 'c_oper_identifiers', 'c_coeffs',
+                      'n_opers', 'n_oper_identifiers', 'n_coeffs',
+                      'dt', 'd', 'basis')
         if not args:
             # Bypass args parsing and directly set necessary attributes
             values = (kwargs[attr] for attr in attributes)
@@ -855,6 +863,7 @@ def get_filter_function_derivative(
             self,
             omega: Coefficients,
             control_identifiers: Optional[Sequence[str]] = None,
+            n_oper_identifiers: Optional[Sequence[str]] = None,
             n_coeffs_deriv: Optional[Sequence[Coefficients]] = None
     ) -> ndarray:
         r"""Calculate the pulse sequence's filter function derivative.
@@ -864,27 +873,70 @@ def get_filter_function_derivative(
         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
+        control_identifiers: Sequence[str], shape (n_ctrl,)
+            Sequence of strings with the control identifiers to
             distinguish between control and drift Hamiltonian. The
-            default is None.
+            default is None, in which case the derivative is computed
+            for all known non-noise operators.
+        n_oper_identifiers: Sequence[str], shape (n_nops,)
+            Sequence of strings with the noise identifiers for which to
+            compute the derivative contribution. The default is None, in
+            which case it is computed for all known noise operators.
         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.
+            amplitudes. The rows and columns should be in the same order
+            as the corresponding identifiers above. Defaults to None, in
+            which case the coefficients are assumed to be constant and
+            hence their derivative vanishing.
+
+            .. warning::
+
+                Internally, control and noise terms of the Hamiltonian
+                are stored alphanumerically sorted by their identifiers.
+                If the noise and/or control identifiers above are not
+                explicitly given, the rows and/or columns of this
+                parameter need to be sorted in the same fashion.
 
         Returns
         -------
         filter_function_deriv: ndarray, shape (n_nops, n_t, n_ctrl, n_omega)
             The regular filter functions' derivatives for variation in
-            each control contribution.
+            each control contribution. Sorted in the same fashion as
+            `n_coeffs_deriv` or, if not given, alphanumerically by the
+            identifiers.
 
         """
-        control_matrix = self.get_control_matrix(omega, cache_intermediates=True)
+        c_idx = util.get_indices_from_identifiers(self.c_oper_identifiers, control_identifiers)
+        n_idx = util.get_indices_from_identifiers(self.n_oper_identifiers, n_oper_identifiers)
+
+        if n_coeffs_deriv is not None:
+            # TODO 05/22: walrus once support for 3.7 is dropped.
+            actual_shape = np.shape(n_coeffs_deriv)
+            required_shape = (len(n_idx), len(c_idx), len(self))
+            if actual_shape != required_shape:
+                raise ValueError(f'Expected n_coeffs_deriv to be of shape {required_shape}, '
+                                 f'not {actual_shape}. Did you forget to specify identifiers?')
+            else:
+                # Do nothing; n_coeffs_deriv specifies the sorting order. If identifiers
+                # are given, we sort everything else by them. If not, n_coeffs_deriv is
+                # expected to be sorted in accordance with the opers and coeffs.
+                pass
 
+        # Check if we can pass on intermediates.
+        intermediates = dict()
+        # TODO 05/22: walrus once support for 3.7 is dropped.
+        n_opers_transformed = self._intermediates.get('n_opers_transformed')
+        first_order_integral = self._intermediates.get('first_order_integral')
+        if n_opers_transformed is not None:
+            intermediates['n_opers_transformed'] = n_opers_transformed[n_idx]
+        if first_order_integral is not None:
+            intermediates['first_order_integral'] = first_order_integral
+
+        control_matrix = self.get_control_matrix(omega, cache_intermediates=True)[n_idx]
         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
+            self.n_opers[n_idx], self.n_coeffs[n_idx], self.c_opers[c_idx], n_coeffs_deriv,
+            intermediates
         )
         return gradient.calculate_filter_function_derivative(control_matrix, control_matrix_deriv)
 

tests/gradient_testutil.py

@@ -21,10 +21,6 @@ def deriv_exchange_interaction(eps):
     return j_0 / eps_0 * np.exp(eps / eps_0)
 
 
-def deriv_2_exchange_interaction(eps):
-    return j_0 / eps_0 ** 2 * np.exp(eps / eps_0)
-
-
 def one_over_f_noise(f):
     spectrum = np.divide(S_0, f, where=(f != 0))
     spectrum[f == 0] = spectrum[np.abs(f).argmin()]
@@ -40,8 +36,6 @@ def create_sing_trip_pulse_seq(eps, dbz, *args):
 
     H_n = [
         [sigma_z, deriv_exchange_interaction(eps[0])]
-        # [sigma_z, np.ones_like(eps[0])]
-
     ]
 
     dt = time_step * np.ones(n_time_steps)
@@ -84,18 +78,17 @@ def finite_diff_infid(u_ctrl_central, u_drift, d, pulse_sequence_builder,
 
     """
     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)]
+        c_id = pulse.c_oper_identifiers[:len(u_ctrl_central)]
 
     # 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)))
+    gradient = np.empty(pulse.n_coeffs.shape + (len(c_id),))
 
-    for g in range(len(pulse.dt)):
+    for g in range(len(pulse)):
         for k in range(len(c_id)):
             u_plus = u_ctrl_central.copy()
             u_plus[k, g] += delta_u