Remove memory_parsimonious option for frequency shifts

Doesn't make much sense since we cannot factor the second order filter
function into control matrices as in the first order case.

In princinple, one can still loop the integration over one basis axis,
but this is a bit tedious and not pretty to implement, so we leave it
for now.

Author: thangleiter

filter_functions/numeric.py

@@ -1037,9 +1037,11 @@ def calculate_cumulant_function(
 
     if second_order:
         if frequency_shifts is None:
+            if memory_parsimonious:
+                warn('Memory parsimonious calculation not implemented for frequency shifts.')
+
             frequency_shifts = calculate_frequency_shifts(pulse, spectrum, omega,
-                                                          n_oper_identifiers, show_progressbar,
-                                                          memory_parsimonious)
+                                                          n_oper_identifiers, show_progressbar)
 
         if frequency_shifts.shape != decay_amplitudes.shape:
             raise ValueError('Frequency shifts not same shape as decay amplitudes')
@@ -1236,8 +1238,7 @@ def calculate_frequency_shifts(
         spectrum: ndarray,
         omega: Coefficients,
         n_oper_identifiers: Optional[Sequence[str]] = None,
-        show_progressbar: bool = False,
-        memory_parsimonious: bool = False
+        show_progressbar: bool = False
 ) -> ndarray:
     r"""
     Get the frequency shifts :math:`\Delta_{\alpha\beta, kl}` for noise
@@ -1262,18 +1263,6 @@ def calculate_frequency_shifts(
         the frequency shifts. The default is all.
     show_progressbar: bool, optional
         Show a progress bar for the calculation.
-    memory_parsimonious: bool, optional
-        For large dimensions, the integrand
-
-        .. math::
-
-            F_{\alpha\beta, kl}^{(2)}(\omega)S_{\alpha\beta}(\omega)
-
-        can consume quite a large amount of memory if set up for all
-        :math:`\alpha,\beta,k,l` at once. If ``True``, it is only set up
-        and integrated for a single :math:`k` at a time and looped over.
-        This is slower but requires much less memory. The default is
-        ``False``.
 
     Raises
     ------
@@ -1309,15 +1298,12 @@ def calculate_frequency_shifts(
     calculate_pulse_correlation_filter_function
     """
     idx = util.get_indices_from_identifiers(pulse, n_oper_identifiers, 'noise')
-    if not memory_parsimonious:
-        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', filter_function=filter_function_2)
-        frequency_shifts = util.integrate(integrand, omega)/(2*np.pi)
-        return frequency_shifts
-
-    raise NotImplementedError
+    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',
+                               filter_function=filter_function_2)
+    frequency_shifts = util.integrate(integrand, omega)/(2*np.pi)
+    return frequency_shifts
 
 
 @util.parse_which_FF_parameter
@@ -1387,7 +1373,6 @@ def calculate_second_order_filter_function(
         n_coeffs: Sequence[Coefficients],
         dt: Coefficients,
         intermediates: Optional[Dict[str, ndarray]] = None,
-        memory_parsimonious: bool = False,
         show_progressbar: bool = False
 ) -> ndarray:
     r"""Calculate the second order filter function for frequency shifts.
@@ -1423,21 +1408,6 @@ def calculate_second_order_filter_function(
         Intermediate terms of the calculation of the control matrix that
         can be reused here. If None (default), they are computed from
         scratch.
-    memory_parsimonious: bool, optional
-
-        .. warning:: Not implemented.
-
-        For large dimensions, the integrand
-
-        .. math::
-
-            F_{\alpha\beta, kl}^{(2)}(\omega)S_{\alpha\beta}(\omega)
-
-        can consume quite a large amount of memory if set up for all
-        :math:`\alpha,\beta,k,l` at once. If ``True``, it is only set up
-        and integrated for a single :math:`k` at a time and looped over.
-        This is slower but requires much less memory. The default is
-        ``False``.
     show_progressbar: bool, optional
         Show a progress bar for the calculation.
 

filter_functions/pulse_sequence.py

@@ -705,8 +705,7 @@ def cache_filter_function(
                 # order == 2
                 filter_function = numeric.calculate_second_order_filter_function(
                     self.eigvals, self.eigvecs, self.propagators, omega, self.basis,
-                    self.n_opers, self.n_coeffs, self.dt, memory_parsimonious=False,
-                    show_progressbar=show_progressbar, intermediates=self._intermediates
+                    self.n_opers, self.n_coeffs, self.dt, self._intermediates, show_progressbar
                 )
 
         self.omega = omega

tests/test_core.py

@@ -675,9 +675,8 @@ def test_second_order_filter_function(self):
             # Test caching
             F_2 = pulse.get_filter_function(omega, order=2)
             F_3 = numeric.calculate_second_order_filter_function(
-                pulse.eigvals, pulse.eigvecs, pulse.propagators, omega, pulse.basis,
-                pulse.n_opers, pulse.n_coeffs, pulse.dt, memory_parsimonious=False,
-                show_progressbar=False, intermediates=None
+                pulse.eigvals, pulse.eigvecs, pulse.propagators, omega, pulse.basis, pulse.n_opers,
+                pulse.n_coeffs, pulse.dt, show_progressbar=False, intermediates=None
             )
             # Make sure first and second order are of same order of magnitude
             rel = np.linalg.norm(F) / np.linalg.norm(F_1)
@@ -853,15 +852,6 @@ def test_calculate_decay_amplitudes(self):
             with self.assertRaises(ValueError):
                 numeric.calculate_decay_amplitudes(pulse, np.tile(spectrum, [1]*i), omega)
 
-    def test_calculate_frequency_shifts(self):
-        """Test raises of numeric.calculate_frequency_shifts"""
-        pulse = testutil.rand_pulse_sequence(2, 1, 1, 1)
-        omega = rng.standard_normal(43)
-        spectrum = rng.standard_normal(43)
-        with self.assertRaises(NotImplementedError):
-            numeric.calculate_frequency_shifts(pulse, spectrum, omega,
-                                               memory_parsimonious=True)
-
     def test_cumulant_function(self):
         pulse = testutil.rand_pulse_sequence(2, 1, 1, 1)
         omega = rng.standard_normal(43)
@@ -895,6 +885,11 @@ def test_cumulant_function(self):
             numeric.calculate_cumulant_function(pulse, spectrum, omega, second_order=True,
                                                 decay_amplitudes=Gamma[1:])
 
+        with self.assertWarns(UserWarning):
+            # Memory parsimonious only works for decay amplitudes
+            numeric.calculate_cumulant_function(pulse, spectrum, omega, second_order=True,
+                                                memory_parsimonious=True)
+
         for d in [2, *rng.integers(2, 7, 5)]:
             pulse = testutil.rand_pulse_sequence(d, 3, 2, 2)
             omega = util.get_sample_frequencies(pulse, n_samples=42)