Merge 580eb01 into 52230d6

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,

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,

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,

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 <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,40 +481,53 @@ 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>`).
 
     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
     --------
     calculate_control_matrix_from_scratch: Control matrix from scratch.
     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

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