diff --git a/doc/source/examples/advanced_concatenation.ipynb b/doc/source/examples/advanced_concatenation.ipynb index 4ae5360..0c06ef8 100644 --- a/doc/source/examples/advanced_concatenation.ipynb +++ b/doc/source/examples/advanced_concatenation.ipynb @@ -142,7 +142,7 @@ "metadata": {}, "source": [ "## Calculating the pulse correlation filter functions\n", - "Finally, we can set up the `PulseSequence`s and compute the filter functions. For the Hadamard pulse, we will set an additional flag during concatenation to calculate the 'pulse correlation filter function' which captures the cross-correlational effects of individual pulses on the total pulse sequence's susceptibility to noise." + "Finally, we can set up the `PulseSequence`s and compute the (fidelity) filter functions. For the Hadamard pulse, we will set an additional flag during concatenation to calculate the 'pulse correlation filter function' which captures the cross-correlational effects of individual pulses on the total pulse sequence's susceptibility to noise." ] }, { @@ -244,7 +244,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.1" + "version": "3.8.3" } }, "nbformat": 4, diff --git a/doc/source/examples/extending_pulses.ipynb b/doc/source/examples/extending_pulses.ipynb index 34aa856..6fb29dd 100644 --- a/doc/source/examples/extending_pulses.ipynb +++ b/doc/source/examples/extending_pulses.ipynb @@ -230,7 +230,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.1" + "version": "3.8.3" } }, "nbformat": 4, diff --git a/doc/source/examples/getting_started.ipynb b/doc/source/examples/getting_started.ipynb index 68b9fff..84fb472 100644 --- a/doc/source/examples/getting_started.ipynb +++ b/doc/source/examples/getting_started.ipynb @@ -20,6 +20,8 @@ "\n", "Therefore, the decoherence can be very intuitively understood and studied by examining the filter function in the frequency domain -- if it is large, the system will be very susceptible to noise at that frequency, while if it is small, noise of that frequency will be 'filtered out'. In this regard, the filter function is very similar to the transfer function of electrical circuits.\n", "\n", + "Other filter functions than the one *fidelity* filter function above which describes phase coherence may be defined in a similar way and obtained as linear combinations of *generalized* filter functions. These functions are defined with respect to an orthonormal operator basis $\\lbrace C_k\\rbrace_{k=0}^{d^2-1}$ and the fidelity filter function, for example, may be obtained from them by tracing out the basis indices, $F(\\omega) = \\mathrm{tr}(F_{kl}(\\omega))$.\n", + "\n", "### The `PulseSequence` class\n", "The central object of this package is the `PulseSequence` class. It is used to represent the control operation implementing a certain Hamiltonian ${H}_c$ and the sensitivities of the noise afflicting the system. More concisely, the total Hamiltonian of the system is modelled as\n", "\n", @@ -46,7 +48,7 @@ "metadata": {}, "source": [ "## A simple example\n", - "To illustrate the instantiation syntax, we would like to compute the filter function for a Free Induction Decay (FID) experiment, which is given by [[Cywiński et al. (2008)]]\n", + "To illustrate the instantiation syntax, we would like to compute the fidelity filter function for a Free Induction Decay (FID) experiment, which is given by [[Cywiński et al. (2008)]]\n", "\n", "$$\n", " F(\\omega\\tau) = \\frac{2\\sin^2\\frac{\\omega\\tau}{2}}{\\omega^2}\n", @@ -162,6 +164,27 @@ "ff.plot_bloch_vector_evolution(FID, psi0=psi_i, n_samples=101, figsize=(4, 4))" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Generalized filter functions\n", + "As mentioned above, we can of course also compute the fidelity filter function from the generalized filter functions. These may be obtained by specifying a corresponding flag when calling `PulseSequence.get_filter_function()`. We can verify that the fidelity filter function above is calculated correctly by taking the trace over the basis elements of the generalized filter functions. By convention, the first two axes of this quantity correspond to the noise operators $B_\\alpha,B_\\beta$, and the third and fourth correspond to the basis elements $C_k,C_l$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "F_fid = FID.get_filter_function(omega, which='fidelity')\n", + "F_gen = FID.get_filter_function(omega, which='generalized')\n", + "\n", + "print(F_fid.shape, F_gen.shape)\n", + "print(np.allclose(F_fid, F_gen.trace(axis1=2, axis2=3)))" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -272,7 +295,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.1" + "version": "3.8.3" } }, "nbformat": 4, diff --git a/filter_functions/numeric.py b/filter_functions/numeric.py index a4130f1..ec8f967 100644 --- a/filter_functions/numeric.py +++ b/filter_functions/numeric.py @@ -420,33 +420,59 @@ def calculate_error_vector_correlation_functions( return u_kl -def calculate_filter_function(R: ndarray) -> ndarray: - """ - Compute the filter function from the control matrix. +@util.parse_which_FF_parameter +def calculate_filter_function(R: ndarray, which: str = 'fidelity') -> ndarray: + r"""Compute the filter function from the control matrix. Parameters ---------- R: array_like, shape (n_nops, d**2, n_omega) The control matrix. + which : str, optional + Which filter function to return. Either 'fidelity' (default) or + 'generalized' (see :ref:`Notes `). Returns ------- - F: ndarray, shape (n_nops, n_nops, n_omega) + F: ndarray, shape (n_nops, n_nops, [d**2, d**2], n_omega) The filter functions for each noise operator correlation. The diagonal corresponds to the filter functions for uncorrelated noise sources. + Notes + ----- + The generalized filter function is given by + + .. math:: + + F_{\alpha\beta,kl}(\omega) = \mathcal{R}_{\alpha k}^\ast(\omega) + \mathcal{R}_{\beta l}(\omega), + + where :math:`\alpha,\beta` are indices counting the noise operators + :math:`B_\alpha` and :math:`k,l` indices counting the basis elements + :math:`C_k`. + + The fidelity filter function is obtained by tracing over the basis indices: + + .. math:: + + F_{\alpha\beta}(\omega) = \sum_{k} F_{\alpha\beta,kk}(\omega). + See Also -------- calculate_control_matrix_from_scratch: Control matrix from scratch. calculate_control_matrix_from_atomic: Control matrix from concatenation. calculate_pulse_correlation_filter_function: Pulse correlations. """ - return np.einsum('iko,jko->ijo', R.conj(), R) + if which == 'fidelity': + return np.einsum('ako,bko->abo', R.conj(), R) + elif which == 'generalized': + return np.einsum('ako,blo->abklo', R.conj(), R) -def calculate_pulse_correlation_filter_function(R: ndarray) -> ndarray: - r""" - Compute the pulse correlation filter function from the control matrix. +@util.parse_which_FF_parameter +def calculate_pulse_correlation_filter_function( + R: ndarray, which: str = 'fidelity') -> ndarray: + r"""Compute pulse correlation filter function from the control matrix. Parameters ---------- @@ -455,23 +481,34 @@ def calculate_pulse_correlation_filter_function(R: ndarray) -> ndarray: Returns ------- - F_pc: ndarray, shape (n_pulses, n_pulses, n_nops, n_nops, n_omega) + F_pc: ndarray, shape (n_pulses, n_pulses, n_nops, n_nops, [d**2, d**2], n_omega) # noqa The pulse correlation filter functions for each pulse and noise operator correlations. The first two axes hold the pulse correlations, the second two the noise correlations. + which : str, optional + Which filter function to return. Either 'fidelity' (default) or + 'generalized' (see :ref:`Notes `). Notes ----- - The pulse correlation filter function is given by + The generalized pulse correlation filter function is given by .. math:: - F_{\alpha\beta}^{(gg')}(\omega) = e^{i\omega(t_{g-1} - t_{g'-1})} - \mathcal{R}^{(g)}(\omega)\mathcal{Q}^{(g-1)} + F_{\alpha\beta,kl}^{(gg')}(\omega) = \bigl[ \mathcal{Q}^{(g'-1)\dagger}\mathcal{R}^{(g')\dagger}(\omega) + \bigr]_{k\alpha} \bigl[ + \mathcal{R}^{(g)}(\omega)\mathcal{Q}^{(g-1)} + \bigr]_{\beta l} e^{i\omega(t_{g-1} - t_{g'-1})}, with :math:`\mathcal{R}^{(g)}` the control matrix of the :math:`g`-th - pulse. + pulse. The fidelity pulse correlation function is obtained by tracing out + the basis indices, + + .. math:: + + F_{\alpha\beta}^{(gg')}(\omega) = + \sum_{k} F_{\alpha\beta,kk}^{(gg')}(\omega) See Also -------- @@ -479,16 +516,18 @@ def calculate_pulse_correlation_filter_function(R: ndarray) -> ndarray: calculate_control_matrix_from_atomic: Control matrix from concatenation. calculate_filter_function: Regular filter function. """ - try: - F_pc = np.einsum('gjko,hlko->ghjlo', R.conj(), R) - except ValueError: + if R.ndim != 4: raise ValueError('Expected R.ndim == 4.') - return F_pc + if which == 'fidelity': + return np.einsum('gako,hbko->ghabo', R.conj(), R) + elif which == 'generalized': + return np.einsum('gako,hblo->ghabklo', R.conj(), R) def diagonalize(H: ndarray, dt: Coefficients) -> Tuple[ndarray]: - r""" + r"""Diagonalize a Hamiltonian. + Diagonalize the Hamiltonian *H* which is piecewise constant during the times given by *dt* and return eigenvalues, eigenvectors, and the cumulative propagators :math:`Q_l`. Note that we calculate in units where diff --git a/filter_functions/plotting.py b/filter_functions/plotting.py index 80eab0b..1ace148 100644 --- a/filter_functions/plotting.py +++ b/filter_functions/plotting.py @@ -293,8 +293,8 @@ def plot_filter_function(pulse: 'PulseSequence', gridspec_kw: Optional[dict] = None, **figure_kw) -> FigureAxesLegend: r""" - Plot the filter function(s) of the given PulseSequence for positive - frequencies. As of now only the diagonal elements of + Plot the fidelity filter function(s) of the given PulseSequence for + positive frequencies. As of now only the diagonal elements of :math:`F_{\alpha\beta}` are implemented, i.e. the filter functions corresponding to uncorrelated noise sources. @@ -412,8 +412,9 @@ def plot_pulse_correlation_filter_function( gridspec_kw: Optional[dict] = None, **figure_kw) -> FigureAxesLegend: r""" - Plot the pulse correlation filter functions of the given PulseSequence if - they were computed during concatenation for positive frequencies. + Plot the fidelity pulse correlation filter functions of the given + PulseSequence if they were computed during concatenation for positive + frequencies. Returns a figure with *n* by *n* subplots where *n* is the number of pulses that were concatenated. As of now only the diagonal elements of diff --git a/filter_functions/pulse_sequence.py b/filter_functions/pulse_sequence.py index bc76963..f5864bc 100644 --- a/filter_functions/pulse_sequence.py +++ b/filter_functions/pulse_sequence.py @@ -18,9 +18,8 @@ # # Contact email: tobias.hangleiter@rwth-aachen.de # ============================================================================= -""" -This module defines the PulseSequence class, the central object of the -formalism. + +"""This module defines the PulseSequence class, the package's central object. Classes ------- @@ -212,9 +211,9 @@ class PulseSequence: eigenvalues and -vectors as well as cumulative propagators get_control_matrix(omega, show_progressbar=False) Calculate the control matrix for frequencies omega - get_filter_function(omega, show_progressbar=False) + get_filter_function(omega, which='fidelity', show_progressbar=False) Calculate the filter function for frequencies omega - get_pulse_correlation_filter_function() + get_pulse_correlation_filter_function(which='fidelity') Get the pulse correlation filter function (only possible if computed during concatenation) propagator_at_arb_t(t) @@ -268,13 +267,15 @@ def __init__(self, *args, **kwargs) -> None: self._HD = None self._HV = None self._Q = None - # Attributes that are cached during calculation of the filter function self._total_phases = None self._total_Q = None self._total_Q_liouville = None self._R = None + self._R_pc = None self._F = None + self._F_kl = None self._F_pc = None + self._F_pc_kl = None def __str__(self): """String method.""" @@ -395,9 +396,12 @@ def is_cached(self, attr: str) -> bool: 'total phases': '_total_phases', 'filter function': '_F', 'fidelity filter function': '_F', + 'generalized filter function': '_F_kl', 'pulse correlation filter function': '_F_pc', 'fidelity pulse correlation filter function': '_F_pc', - 'control matrix': '_R'} + 'generalized pulse correlation filter function': '_F_pc_kl', + 'control matrix': '_R', + 'pulse correlation control matrix': '_R_pc'} alias = attr.lower().replace('_', ' ') if alias in aliases: @@ -444,13 +448,14 @@ def get_control_matrix(self, omega: Coefficients, if np.array_equal(self.omega, omega): if self.is_cached('R'): return self._R + else: + if self.is_cached('R_pc'): + self._R = self._R_pc.sum(axis=0) + return self._R else: # Getting with different frequencies. Remove all cached attributes # that are frequency-dependent - self._R = None - self._F = None - self._F_pc = None - self._total_phases = None + self.cleanup('frequency dependent') # Make sure the Hamiltonian has been diagonalized self.diagonalize() @@ -475,7 +480,7 @@ def cache_control_matrix(self, omega: Coefficients, ---------- omega: array_like, shape (n_omega,) The frequencies for which to cache the filter function. - R: array_like, shape (n_nops [, n_nops], d**2, n_omega), optional + R: array_like, shape (n_nops, [n_nops,] d**2, n_omega), optional The control matrix for the frequencies *omega*. If ``None``, it is computed. show_progressbar: bool @@ -485,7 +490,11 @@ def cache_control_matrix(self, omega: Coefficients, R = self.get_control_matrix(omega, show_progressbar) self.omega = omega - self._R = R + if R.ndim == 4: + # Pulse correlation control matrix + self._R_pc = R + else: + self._R = R # Cache total phase and total transfer matrices as well self.cache_total_phases(omega) @@ -494,20 +503,22 @@ def cache_control_matrix(self, omega: Coefficients, self.total_Q, self.basis ) - def get_filter_function(self, omega: Coefficients, - show_progressbar: bool = False) -> ndarray: - r""" - Calculate the first-order filter function + def get_pulse_correlation_control_matrix(self) -> ndarray: + """Get the pulse correlation control matrix if it was cached.""" + if self.is_cached('R_pc'): + return self._R_pc - .. math:: - - F_{\alpha\beta}(\omega) = \left[\mathcal{R}\mathcal{R}^\dagger - \right]_{\alpha\beta}(\omega), + raise util.CalculationError( + "Could not get the pulse correlation control matrix since it " + + "was not computed during concatenation. Please run the " + + "concatenation again with 'calc_pulse_correlation_ff' set to " + + "True." + ) - where :math:`\alpha,\beta` are indices counting the noise operators - :math:`B_\alpha`. Thus, the filter function :math:`\alpha,\beta` - corresponds to the noise correlations between :math:`B_\alpha` and - :math:`B_{\beta}`. + @util.parse_which_FF_parameter + def get_filter_function(self, omega: Coefficients, which: str = 'fidelity', + show_progressbar: bool = False) -> ndarray: + r"""Get the first-order filter function. The filter function is cached so it doesn't need to be calculated twice for the same frequencies. @@ -516,37 +527,67 @@ def get_filter_function(self, omega: Coefficients, ---------- omega: array_like, shape (n_omega,) The frequencies at which to evaluate the filter function. - show_progressbar: bool + which: str, optional + Which filter function to return. Either 'fidelity' (default) or + 'generalized' (see :ref:`Notes `). + show_progressbar: bool, optional Show a progress bar for the calculation of the control matrix. Returns ------- - F: ndarray, shape (n_nops, n_nops, n_omega) + F: ndarray, shape (n_nops, n_nops, [d**2, d**2,] n_omega) The filter function for each combination of noise operators as a function of omega. + + Notes + ----- + The generalized filter function is given by + + .. math:: + + F_{\alpha\beta,kl}(\omega) = \mathcal{R}_{\alpha k}^\ast(\omega) + \mathcal{R}_{\beta l}(\omega), + + where :math:`\alpha,\beta` are indices counting the noise operators + :math:`B_\alpha` and :math:`k,l` indices counting the basis elements + :math:`C_k`. + + The fidelity filter function is obtained by tracing over the basis + indices: + + .. math:: + + F_{\alpha\beta}(\omega) = \sum_{k} F_{\alpha\beta,kk}(\omega). + """ # Only calculate if not calculated before for the same frequencies if np.array_equal(self.omega, omega): - if self.is_cached('F'): - return self._F + if which == 'fidelity': + if self.is_cached('F'): + return self._F + elif which == 'generalized': + if self.is_cached('F_kl'): + return self._F_kl else: # Getting with different frequencies. Remove all cached attributes # that are frequency-dependent - self._R = None - self._F = None - self._F_pc = None - self._total_phases = None + self.cleanup('frequency dependent') - # Cache filter function self.cache_filter_function( - omega, R=self.get_control_matrix(omega, show_progressbar) + omega, R=self.get_control_matrix(omega, show_progressbar), + which=which ) - return self._F + if which == 'fidelity': + return self._F + elif which == 'generalized': + return self._F_kl + @util.parse_which_FF_parameter def cache_filter_function(self, omega: Coefficients, R: Optional[ndarray] = None, F: Optional[ndarray] = None, + which: str = 'fidelity', show_progressbar: bool = False) -> None: r""" Cache the filter function. If R.ndim == 4, it is taken to be the 'pulse @@ -559,39 +600,52 @@ def cache_filter_function(self, omega: Coefficients, ---------- omega: array_like, shape (n_omega,) The frequencies for which to cache the filter function. - R: array_like, shape (n_nops [, n_nops], d**2, n_omega), optional + R: array_like, shape (n_nops, [n_nops,] d**2, n_omega), optional The control matrix for the frequencies *omega*. If ``None``, it is computed and the filter function derived from it. - F: array_like, shape (n_nops, n_nops, n_omega), optional + F: array_like, shape (n_nops, n_nops, [d**2, d**2,] n_omega), optional The filter function for the frequencies *omega*. If ``None``, it is computed from R. + which: str, optional + Which filter function to return. Either 'fidelity' (default) or + 'generalized'. show_progressbar: bool Show a progress bar for the calculation of the control matrix. + + See Also + -------- + PulseSequence.get_filter_function : Getter method """ if F is None: if R is None: R = self.get_control_matrix(omega, show_progressbar) + self.cache_control_matrix(omega, R) if R.ndim == 4: - # Cache regular control matrix - self.cache_control_matrix(omega, R.sum(axis=0), False) # Calculate pulse correlation FF and derive canonical FF from # it - self._F_pc = \ - numeric.calculate_pulse_correlation_filter_function(R) + F_pc = numeric.calculate_pulse_correlation_filter_function( + R, which) - F = np.einsum('ghjlo->jlo', self._F_pc) + if which == 'fidelity': + self._F_pc = F_pc + elif which == 'generalized': + self._F_pc_kl = F_pc + F = F_pc.sum(axis=(0, 1)) else: # Regular case - self.cache_control_matrix(omega, R, - show_progressbar=show_progressbar) - F = numeric.calculate_filter_function(R) + F = numeric.calculate_filter_function(R, which=which) self.omega = omega - self._F = F - - def get_pulse_correlation_filter_function(self) -> ndarray: + if which == 'fidelity': + self._F = F + elif which == 'generalized': + self._F_kl = F + + @util.parse_which_FF_parameter + def get_pulse_correlation_filter_function( + self, which: str = 'fidelity') -> ndarray: r""" Get the pulse correlation filter function given by @@ -632,13 +686,28 @@ def get_pulse_correlation_filter_function(self) -> ndarray: for :math:`g,g'\in\{A, B, C\}`. """ - if self.is_cached('F_pc'): - return self._F_pc + if which == 'fidelity': + if self.is_cached('F_pc'): + return self._F_pc + elif which == 'generalized': + if self.is_cached('F_pc_kl'): + return self._F_pc_kl + + if self.is_cached('R_pc'): + F_pc = numeric.calculate_pulse_correlation_filter_function( + self._R_pc, which) + + if which == 'fidelity': + self._F_pc = F_pc + elif which == 'generalized': + self._F_pc_kl = F_pc + + return F_pc raise util.CalculationError( - "Could not get the pulse correlation function since it was not " + - "computed during concatenation. Please run the concatenation " + - "again with 'calc_pulse_correlation_ff' set to True." + "Could not get the pulse correlation filter function since it " + + "was not computed during concatenation. Please run the " + + "concatenation again with 'calc_pulse_correlation_ff' set to True." ) def get_total_phases(self, omega: Coefficients) -> ndarray: @@ -650,10 +719,7 @@ def get_total_phases(self, omega: Coefficients) -> ndarray: else: # Getting with different frequencies. Remove all cached attributes # that are frequency-dependent - self._R = None - self._F = None - self._F_pc = None - self._total_phases = None + self.cleanup('frequency dependent') self.cache_total_phases(omega) return self._total_phases @@ -768,6 +834,8 @@ def nbytes(self) -> int: return sum(_nbytes) + @util.parse_optional_parameter( + 'method', ('conservative', 'greedy', 'frequency dependent', 'all')) def cleanup(self, method: str = 'conservative') -> None: """ Delete cached byproducts of the calculation of the filter function that @@ -775,7 +843,7 @@ def cleanup(self, method: str = 'conservative') -> None: Parameters ---------- - method: {'conservative', 'greedy', 'all'}, optional + method: optional If set to 'conservative' (the default), only the following attributes are deleted: @@ -786,28 +854,43 @@ def cleanup(self, method: str = 'conservative') -> None: If set to 'greedy', all of the above as well as the following attributes are deleted: - - _total_phases - _total_Q - _total_Q_liouville + - _total_phases - _R + - _R_pc If set to 'all', all of the above as well as the following attributes are deleted: - omega - _F + - _F_kl - _F_pc + - _F_pc_kl + + If set to 'frequency dependent' only attributes that are functions + of frequency are initalized to ``None``. Note that if this ``PulseSequence`` is concatenated with another one, some of the attributes might need to be calculated again, resulting in slower execution of the concatenation. """ - attrs = ['_HD', '_HV', '_Q'] - if method != 'conservative': - attrs.extend(['_R', '_total_phases', '_total_Q', - '_total_Q_liouville']) - if method != 'greedy': - attrs.extend(['omega', '_F', '_F_pc']) + default_attrs = {'_HD', '_HV', '_Q'} + concatenation_attrs = {'_total_Q', '_total_Q_liouville', '_R', '_R_pc', + '_total_phases'} + filter_function_attrs = {'omega', '_F', '_F_kl', '_F_pc', '_F_pc_kl'} + + if method == 'conservative': + attrs = default_attrs + elif method == 'greedy': + attrs = default_attrs.union(concatenation_attrs) + elif method == 'frequency dependent': + attrs = filter_function_attrs.union({'_R', '_R_pc', + '_total_phases'}) + else: + attrs = filter_function_attrs.union(default_attrs, + concatenation_attrs) for attr in attrs: setattr(self, attr, None) @@ -858,7 +941,7 @@ def _parse_args(H_c: Hamiltonian, H_n: Hamiltonian, dt: Coefficients, object. """ - if not hasattr(dt, '__getitem__'): + if not hasattr(dt, '__len__'): raise TypeError('Expected a sequence of time steps, not {}'.format( type(dt))) @@ -932,7 +1015,7 @@ def _parse_Hamiltonian(H: Hamiltonian, n_dt: int, raise TypeError('Expected operators in '.format(H_str) + 'to be NumPy arrays or QuTiP Qobjs!') - if not all(hasattr(coeff, '__getitem__') for coeff in coeffs): + if not all(hasattr(coeff, '__len__') for coeff in coeffs): raise TypeError('Expected coefficients in '.format(H_str) + 'to be a sequence') @@ -949,10 +1032,6 @@ def _parse_Hamiltonian(H: Hamiltonian, n_dt: int, raise ValueError('Expected all operators in {} '.format(H_str) + 'to be two-dimensional!') - if len(set(oper.shape for oper in opers)) != 1: - raise ValueError('Expected all operators in {} '.format(H_str) + - 'to have the same dimensions!') - if len(set(opers[0].shape)) != 1: raise ValueError('Expected operators in {} '.format(H_str) + 'to be square!') @@ -1346,9 +1425,11 @@ def concatenate_without_filter_function( return newpulse +@util.parse_which_FF_parameter def concatenate(pulses: Iterable[PulseSequence], calc_pulse_correlation_ff: bool = False, calc_filter_function: Optional[bool] = None, + which: str = 'fidelity', omega: Optional[Coefficients] = None, show_progressbar: bool = False) -> PulseSequence: r""" @@ -1372,14 +1453,19 @@ def concatenate(pulses: Iterable[PulseSequence], calculation of the composite filter function is forced. calc_pulse_correlation_ff: bool, optional Switch to control whether the pulse correlation filter function (see - :ref:`Notes `) is calculated. If *omega* is not given, the - cached frequencies of all *pulses* need to be equal. + :meth:`PulseSequence.get_pulse_correlation_filter_function`) is + calculated. If *omega* is not given, the cached frequencies of all + *pulses* need to be equal. calc_filter_function: bool, optional Switch to force the calculation of the filter function to be carried out or not. Overrides the automatic behavior of calculating it if at least one pulse has a cached control matrix. If ``True`` and no pulse has a cached control matrix, a list of frequencies must be supplied as *omega*. + which: str, optional + Which filter function to compute. Either 'fidelity' (default) or + 'generalized' (see :meth:`PulseSequence.get_filter_function` and + :meth:`PulseSequence.get_pulse_correlation_filter_function`). omega: array_like, optional Frequencies at which to evaluate the (pulse correlation) filter functions. If ``None``, an attempt is made to use cached frequencies. @@ -1391,22 +1477,6 @@ def concatenate(pulses: Iterable[PulseSequence], pulse: PulseSequence The concatenated pulse. - .. _notes: - - Notes - ----- - The pulse correlation filter function is given by - - .. math:: - - F_{\alpha\beta}^{(gg')}(\omega) = e^{i\omega(t_{g-1} - t_{g'-1})} - \left[\mathcal{Q}^{(g'-1)\dagger} - \mathcal{R}^{(g')\dagger}(\omega)\right]_{k\alpha} - \left[\mathcal{R}^{(g)}(\omega) - \mathcal{Q}^{(g-1)}\right]_{\beta l}, - - where :math:`g,g'` index the pulse in the sequence and :math:`\alpha,\beta` - index the noise operators. """ pulses = tuple(pulses) if len(pulses) == 1: @@ -1448,18 +1518,25 @@ def concatenate(pulses: Iterable[PulseSequence], equal_n_opers = (n_opers_present.sum(axis=0) > 1).any() if omega is None: cached_ctrl_mat = [pls.is_cached('R') for pls in pulses] - equal_omega = util.all_array_equal( - (pls.omega for pls in compress(pulses, cached_ctrl_mat)) - ) + if any(cached_ctrl_mat): + equal_omega = util.all_array_equal( + (pls.omega for pls in compress(pulses, cached_ctrl_mat)) + ) + else: + cached_omega = [pls.is_cached('omega') for pls in pulses] + equal_omega = util.all_array_equal( + (pls.omega for pls in compress(pulses, cached_omega)) + ) + if not equal_omega: if calc_filter_function: raise ValueError("Calculation of filter function forced " + - "but not all pulses have the same " + - "frequencies cached and none were supplied!") + "but not all pulses have the same " + + "frequencies cached and none were supplied!") if calc_pulse_correlation_ff: raise ValueError("Cannot compute the pulse correlation " + - "filter functions; do not have the " + - "frequencies at which to evaluate.") + "filter functions; do not have the " + + "frequencies at which to evaluate.") return newpulse @@ -1472,12 +1549,16 @@ def concatenate(pulses: Iterable[PulseSequence], # Can reuse cached filter functions or calculation explicitly asked # for; run calculation. Get the index of the first pulse with cached FF # to steal some attributes from. - ind = np.nonzero(cached_ctrl_mat)[0][0] + if any(cached_ctrl_mat): + ind = np.nonzero(cached_ctrl_mat)[0][0] + else: + ind = np.nonzero(cached_omega)[0][0] + omega = pulses[ind].omega if not equal_n_opers: # Cannot reuse atomic filter functions - newpulse.cache_filter_function(omega) + newpulse.cache_filter_function(omega, which=which) return newpulse # Get the phase factors at the correct times (the individual gate @@ -1531,7 +1612,7 @@ def concatenate(pulses: Iterable[PulseSequence], # Set the attribute and calculate filter function (if the pulse correlation # FF has been calculated, this is a little overhead but negligible) - newpulse.cache_filter_function(omega, R) + newpulse.cache_filter_function(omega, R, which=which) return newpulse diff --git a/filter_functions/util.py b/filter_functions/util.py index 3ce102b..ded367f 100644 --- a/filter_functions/util.py +++ b/filter_functions/util.py @@ -70,6 +70,7 @@ """ import functools +import inspect import io import json import operator @@ -1098,6 +1099,36 @@ def progressbar_range(*args, show_progressbar: Optional[bool] = True, return range(*args) +def parse_optional_parameter(name: str, allowed: Sequence) -> Callable: + """Decorator factory to parse optional parameter with certain legal values. + + If the parameter value corresponding to ``name`` (either in args or kwargs) + is not contained in ``allowed`` a ``ValueError`` is raised. + """ + def decorator(func): + @functools.wraps(func) + def wrapper(*args, **kwargs): + parameters = inspect.signature(func).parameters + idx = tuple(parameters).index(name) + try: + value = args[idx] + except IndexError: + value = kwargs.get(name, parameters[name].default) + + if value not in allowed: + raise ValueError( + "Invalid value for {}: {}. ".format(name, value) + + "Should be one of {}.".format(allowed) + ) + return func(*args, **kwargs) + return wrapper + return decorator + + +parse_which_FF_parameter = parse_optional_parameter( + 'which', ('fidelity', 'generalized')) + + class CalculationError(Exception): """Indicates a quantity could not be computed.""" diff --git a/tests/test_core.py b/tests/test_core.py index 2aa4cd0..19b1176 100644 --- a/tests/test_core.py +++ b/tests/test_core.py @@ -28,7 +28,7 @@ import numpy as np import filter_functions as ff -from filter_functions import numeric +from filter_functions import numeric, util from tests import testutil @@ -63,7 +63,7 @@ def test_pulse_sequence_constructor(self): dt[idx] *= -1 with self.assertRaises(ValueError): - # imagniary dt + # imaginary dt dt = dt.astype(complex) dt[idx] += 1j ff.PulseSequence(H_c, H_n, dt) @@ -86,6 +86,14 @@ def test_pulse_sequence_constructor(self): # Noise Hamiltonian not list or tuple ff.PulseSequence(H_c, np.array(H_n), dt) + with self.assertRaises(TypeError): + # Element of control Hamiltonian not list or tuple + ff.PulseSequence([np.array(H_c[0])], H_n, dt) + + with self.assertRaises(TypeError): + # Element of noise Hamiltonian not list or tuple + ff.PulseSequence(H_c, [np.array(H_n[0])], dt) + idx = testutil.rng.randint(0, 3) with self.assertRaises(TypeError): # Control Hamiltonian element not list or tuple @@ -171,6 +179,14 @@ def test_pulse_sequence_constructor(self): for hn in H_n: hn[0] = hn[0].reshape(2, 2) + with self.assertRaises(ValueError): + # Control and noise operators not same dimension + for hn in H_n: + hn[0] = np.block([[hn[0], hn[0]], [hn[0], hn[0]]]) + ff.PulseSequence(H_c, H_n, dt) + + for hn in H_n: + hn[0] = hn[0][:2, :2] with self.assertRaises(ValueError): # Control identifiers not unique identifier = H_c[idx][2] @@ -213,10 +229,10 @@ def test_pulse_sequence_constructor(self): # Fewer identifiers than opers pulse_2 = ff.PulseSequence( - [[ff.util.P_np[1], [1], 'X'], - [ff.util.P_np[2], [1]]], - [[ff.util.P_np[1], [1]], - [ff.util.P_np[2], [1], 'Y']], + [[util.P_np[1], [1], 'X'], + [util.P_np[2], [1]]], + [[util.P_np[1], [1]], + [util.P_np[2], [1], 'Y']], [1] ) self.assertArrayEqual(pulse_2.c_oper_identifiers, ('A_1', 'X')) @@ -224,7 +240,7 @@ def test_pulse_sequence_constructor(self): def test_pulse_sequence_attributes(self): """Test attributes of single instance""" - X, Y, Z = ff.util.P_np[1:] + X, Y, Z = util.P_np[1:] n_dt = testutil.rng.randint(1, 10) # trivial case @@ -329,10 +345,23 @@ def test_pulse_sequence_attributes(self): setattr(A, attr, 'foo') assertion = self.assertTrue else: + setattr(A, attr, None) assertion = self.assertFalse assertion(A.is_cached(attr)) - setattr(A, attr, None) + + # Diagonalization attributes + A.diagonalize() + self.assertIsNotNone(A.HD) + self.assertIsNotNone(A.HV) + self.assertIsNotNone(A.Q) + + A.cleanup('conservative') + self.assertIsNotNone(A.HD) + A.cleanup('conservative') + self.assertIsNotNone(A.HV) + A.cleanup('conservative') + self.assertIsNotNone(A.Q) aliases = {'eigenvalues': '_HD', 'eigenvectors': '_HV', @@ -343,9 +372,12 @@ def test_pulse_sequence_attributes(self): 'total phases': '_total_phases', 'filter function': '_F', 'fidelity filter function': '_F', + 'generalized filter function': '_F_kl', 'pulse correlation filter function': '_F_pc', 'fidelity pulse correlation filter function': '_F_pc', - 'control matrix': '_R'} + 'generalized pulse correlation filter function': '_F_pc_kl', + 'control matrix': '_R', + 'pulse correlation control matrix': '_R_pc'} for alias, attr in aliases.items(): # set mock attribute at random @@ -353,19 +385,20 @@ def test_pulse_sequence_attributes(self): setattr(A, attr, 'foo') assertion = self.assertTrue else: + setattr(A, attr, None) assertion = self.assertFalse assertion(A.is_cached(alias)) assertion(A.is_cached(alias.upper())) assertion(A.is_cached(alias.replace(' ', '_'))) - setattr(A, attr, None) + A.cleanup('all') # Test cleanup C = ff.concatenate((A, A), calc_pulse_correlation_ff=True, - omega=ff.util.get_sample_frequencies(A)) + which='generalized', + omega=util.get_sample_frequencies(A)) C.diagonalize() - C.cache_filter_function(ff.util.get_sample_frequencies(A)) attrs = ['_HD', '_HV', '_Q'] for attr in attrs: self.assertIsNotNone(getattr(C, attr)) @@ -375,6 +408,7 @@ def test_pulse_sequence_attributes(self): self.assertIsNone(getattr(C, attr)) C.diagonalize() + C.cache_control_matrix(A.omega) attrs.extend(['_R', '_total_phases', '_total_Q', '_total_Q_liouville']) for attr in attrs: self.assertIsNotNone(getattr(C, attr)) @@ -383,18 +417,35 @@ def test_pulse_sequence_attributes(self): for attr in attrs: self.assertIsNone(getattr(C, attr)) - C.cache_filter_function(ff.util.get_sample_frequencies(A)) + C.cache_filter_function(A.omega, which='generalized') + for attr in attrs + ['omega', '_F_kl', '_F_pc_kl']: + self.assertIsNotNone(getattr(C, attr)) + + C = ff.concatenate((A, A), calc_pulse_correlation_ff=True, + which='fidelity', omega=A.omega) + C.diagonalize() + C.cache_filter_function(A.omega, which='fidelity') attrs.extend(['omega', '_F', '_F_pc']) for attr in attrs: self.assertIsNotNone(getattr(C, attr)) C.cleanup('all') - for attr in attrs: + for attr in attrs + ['_F_kl', '_F_pc_kl']: self.assertIsNone(getattr(C, attr)) + C.cache_filter_function(A.omega, which='fidelity') + C.cleanup('frequency dependent') + freq_attrs = {'omega', '_R', '_F', '_F_kl', '_F_pc', '_F_pc_kl', + '_total_phases'} + for attr in freq_attrs: + self.assertIsNone(getattr(C, attr)) + + for attr in set(attrs).difference(freq_attrs): + self.assertIsNotNone(getattr(C, attr)) + def test_pulse_sequence_attributes_concat(self): """Test attributes of concatenated sequence.""" - X, Y, Z = ff.util.P_np[1:] + X, Y, Z = util.P_np[1:] n_dt_1 = testutil.rng.randint(5, 11) x_coeff_1 = testutil.rng.randn(n_dt_1) z_coeff_1 = testutil.rng.randn(n_dt_1) @@ -414,6 +465,13 @@ def test_pulse_sequence_attributes_concat(self): [[Z, np.abs(testutil.rng.randn(2))]], [1, 1]) + # Concatenate with different noise opers + pulses = [testutil.rand_pulse_sequence(2, 1) for _ in range(2)] + pulses[0].omega = np.arange(10) + pulses[1].omega = np.arange(10) + newpulse = ff.concatenate(pulses, calc_filter_function=True) + self.assertTrue(newpulse.is_cached('filter function')) + pulse_12 = pulse_1 @ pulse_2 pulse_21 = pulse_2 @ pulse_1 @@ -493,7 +551,7 @@ def test_filter_function(self): total_HD, total_HV, _ = numeric.diagonalize( np.einsum('il,ijk->ljk', c_coeffs, c_opers), total_pulse.dt ) - omega = ff.util.get_sample_frequencies(total_pulse, n_samples=100) + omega = util.get_sample_frequencies(total_pulse, n_samples=100) # Try the progress bar R = total_pulse.get_control_matrix(omega, show_progressbar=True) @@ -543,9 +601,45 @@ def test_filter_function(self): np.isreal(F[np.eye(len(n_opers), dtype=bool)]).all() ) + # Check switch between fidelity and generalized filter function + F_generalized = total_pulse.get_filter_function( + omega, which='generalized') + + F_fidelity = total_pulse.get_filter_function( + omega, which='fidelity') + + # Check that F_fidelity is correctly reduced from F_generalized + self.assertArrayAlmostEqual(F_fidelity, + F_generalized.trace(axis1=2, axis2=3)) + + # Hit getters again to check caching functionality + F_generalized = total_pulse.get_filter_function( + omega, which='generalized') + + F_fidelity = total_pulse.get_filter_function( + omega, which='fidelity') + + # Check that F_fidelity is correctly reduced from F_generalized + self.assertArrayAlmostEqual(F_fidelity, + F_generalized.trace(axis1=2, axis2=3)) + + # Different set of frequencies than cached + F_generalized = total_pulse.get_filter_function( + omega + 1, which='generalized') + + F_fidelity = total_pulse.get_filter_function( + omega + 1, which='fidelity') + + # Check that F_fidelity is correctly reduced from F_generalized + self.assertArrayAlmostEqual(F_fidelity, + F_generalized.trace(axis1=2, axis2=3)) + def test_pulse_correlation_filter_function(self): - """Test calculation of pulse correlation filter function""" - X, Y, Z = ff.util.P_np[1:] + """ + Test calculation of pulse correlation filter function and control + matrix. + """ + X, Y, Z = util.P_np[1:] T = 1 omega = np.linspace(-2e1, 2e1, 250) H_c, H_n, dt = dict(), dict(), dict() @@ -577,10 +671,16 @@ def test_pulse_correlation_filter_function(self): pulse_1 = pulses['X'] @ pulses['Y'] pulse_2 = ff.concatenate([pulses['X'], pulses['Y']], - calc_pulse_correlation_ff=True) + calc_pulse_correlation_ff=True, + which='fidelity') + pulse_3 = ff.concatenate([pulses['X'], pulses['Y']], + calc_pulse_correlation_ff=True, + which='generalized') - with self.assertRaises(ValueError): - numeric.calculate_pulse_correlation_filter_function(pulse_1._R) + self.assertTrue(pulse_2.is_cached('R_pc')) + self.assertTrue(pulse_2.is_cached('F_pc')) + self.assertTrue(pulse_3.is_cached('R_pc')) + self.assertTrue(pulse_3.is_cached('F_pc_kl')) # Check if the filter functions on the diagonals are real F = pulse_2.get_pulse_correlation_filter_function() @@ -598,8 +698,58 @@ def test_pulse_correlation_filter_function(self): self.assertArrayAlmostEqual(pulse_1.get_filter_function(omega), pulse_2._F) + # Test the behavior of the pulse correlation control matrix + with self.assertRaises(ValueError): + # R wrong dimension + numeric.calculate_pulse_correlation_filter_function(pulse_1._R) + + with self.assertRaises(util.CalculationError): + # not calculated + pulse_1.get_pulse_correlation_control_matrix() + + R_pc = pulse_2.get_pulse_correlation_control_matrix() + self.assertArrayEqual( + F, numeric.calculate_pulse_correlation_filter_function( + R_pc, 'fidelity' + ) + ) + + R_pc = pulse_3.get_pulse_correlation_control_matrix() + F = pulse_3.get_pulse_correlation_filter_function(which='fidelity') + self.assertArrayEqual( + F, numeric.calculate_pulse_correlation_filter_function( + R_pc, 'fidelity' + ) + ) + + F = pulse_3.get_pulse_correlation_filter_function(which='generalized') + self.assertArrayEqual( + F, numeric.calculate_pulse_correlation_filter_function( + R_pc, 'generalized' + ) + ) + + # If for some reason F_pc_xy is removed, check if recovered from R_pc + pulse_2._F_pc = None + pulse_3._F_pc_kl = None + + R_pc = pulse_3.get_pulse_correlation_control_matrix() + F = pulse_3.get_pulse_correlation_filter_function(which='fidelity') + self.assertArrayEqual( + F, numeric.calculate_pulse_correlation_filter_function( + R_pc, 'fidelity' + ) + ) + + F = pulse_3.get_pulse_correlation_filter_function(which='generalized') + self.assertArrayEqual( + F, numeric.calculate_pulse_correlation_filter_function( + R_pc, 'generalized' + ) + ) + S = omega**0*1e-2 - with self.assertRaises(ff.util.CalculationError): + with self.assertRaises(util.CalculationError): infid_1 = ff.infidelity(pulse_1, S, omega, which='correlations') with self.assertRaises(ValueError): diff --git a/tests/test_sequencing.py b/tests/test_sequencing.py index 1d89a78..a0d7a95 100644 --- a/tests/test_sequencing.py +++ b/tests/test_sequencing.py @@ -22,6 +22,7 @@ This module tests the concatenation functionality for PulseSequence's. """ +from copy import copy import string from itertools import product from random import sample @@ -29,12 +30,37 @@ import numpy as np import filter_functions as ff -from filter_functions import util +from filter_functions import pulse_sequence, util from tests import testutil class ConcatenationTest(testutil.TestCase): + def test_concatenate_base(self): + """Basic functionality.""" + pulse_1, pulse_2 = [testutil.rand_pulse_sequence(2, 1, 2, 3) + for _ in range(2)] + + # Trivial case, copy + c_pulse = ff.concatenate([pulse_1]) + self.assertEqual(pulse_1, c_pulse) + self.assertFalse(pulse_1 is c_pulse) + + # Don't cache filter function, expect same result as with + # concatenate_without_filter_function + c_pulse_1 = ff.concatenate([pulse_1, pulse_2], + calc_filter_function=False) + c_pulse_2 = pulse_sequence.concatenate_without_filter_function( + [pulse_1, pulse_2], return_identifier_mappings=False + ) + self.assertEqual(c_pulse_1, c_pulse_2) + + # Try concatenation with different frequencies but FF calc. forced + with self.assertRaises(ValueError): + pulse_1.omega = [1, 2] + pulse_2.omega = [3, 4] + ff.concatenate([pulse_1, pulse_2], calc_filter_function=True) + def test_concatenate_without_filter_function(self): """Concatenate two Spin Echos without filter functions.""" tau = 10 @@ -67,6 +93,35 @@ def test_concatenate_without_filter_function(self): CPMG_concat = ff.concatenate((SE_1, SE_2), omega=omega) self.assertIsNotNone(CPMG_concat._F) + pulse = testutil.rand_pulse_sequence(2, 1, 2, 3) + # Concatenate pulses without filter functions + with self.assertRaises(TypeError): + # Not all pulse sequence + pulse_sequence.concatenate_without_filter_function([pulse, 2]) + + with self.assertRaises(TypeError): + # Not iterable + pulse_sequence.concatenate_without_filter_function(pulse) + + with self.assertRaises(ValueError): + # Incompatible Hamiltonian shapes + pulse_sequence.concatenate_without_filter_function( + [testutil.rand_pulse_sequence(2, 1), + testutil.rand_pulse_sequence(3, 1)] + ) + + with self.assertRaises(ValueError): + # Incompatible bases + pulse = testutil.rand_pulse_sequence(4, 1, btype='GGM') + cpulse = copy(pulse) + cpulse.basis = ff.Basis.pauli(2) + pulse_sequence.concatenate_without_filter_function([pulse, cpulse]) + + pulse = pulse_sequence.concatenate_without_filter_function( + [pulse, pulse], return_identifier_mappings=False + ) + self.assertFalse(pulse.is_cached('filter function')) + def test_concatenate_with_filter_function_SE1(self): """ Concatenate two Spin Echos with the first having a filter function. @@ -312,7 +367,7 @@ def test_concatenate_split_cnot(self): rtol, atol) def test_different_n_opers(self): - """Test behavior when concatenating with different n_opers""" + """Test behavior when concatenating with different n_opers.""" for d, n_dt in zip(testutil.rng.randint(2, 5, 20), testutil.rng.randint(1, 11, 20)): opers = testutil.rand_herm_traceless(d, 10) diff --git a/tests/test_util.py b/tests/test_util.py index 1001dab..7396152 100644 --- a/tests/test_util.py +++ b/tests/test_util.py @@ -583,3 +583,27 @@ def test_progressbar_range(self): ii.append(i) self.assertEqual(ii, list(range(523, 123, -32))) + + def test_parse_optional_parameter(self): + + @util.parse_optional_parameter('foo', [1, 'bar', (2, 3)]) + def foobar(a, b, foo=None, x=2): + pass + + with self.assertRaises(ValueError) as err: + foobar(1, 1, 2) + self.assertEqual(str(err.exception), + "Invalid value for foo: {}.".format(2) + + " Should be one of {}".format([1, 'bar', (2, 3)])) + + with self.assertRaises(ValueError): + foobar(1, 1, 'x') + self.assertEqual(str(err.exception), + "Invalid value for foo: {}.".format('x') + + " Should be one of {}".format([1, 'bar', (2, 3)])) + + with self.assertRaises(ValueError): + foobar(1, 1, [1, 2]) + self.assertEqual(str(err.exception), + "Invalid value for foo: {}.".format([1, 2]) + + " Should be one of {}".format([1, 'bar', (2, 3)]))