Enhance documentation

filter_functions/basis.py

@@ -318,11 +318,8 @@ def H(self) -> 'Basis':
     @property
     def T(self) -> 'Basis':
         """Return the basis transposed element-wise."""
-        if self.ndim == 3:
-            return self.transpose(0, 2, 1)
-
-        if self.ndim == 2:
-            return self.transpose(1, 0)
+        if self.ndim >= 2:
+            return self.swapaxes(-1, -2)
 
         return self
 

filter_functions/numeric.py

@@ -1108,23 +1108,51 @@ def calculate_cumulant_function(
             cumulant_function[..., 1:, 1:] -= frequency_shifts[..., 1:, 1:]
             cumulant_function[..., 1:, 1:] += frequency_shifts[..., 1:, 1:].swapaxes(-1, -2)
     else:
-        # Multi qubit case. Use general expression. Drop imaginary part since
-        # result is guaranteed to be real (if we didn't do anything wrong)
+        # Multi qubit case. Use general expression.
         traces = pulse.basis.four_element_traces
+        # g_iklj = (
+        #     + traces.transpose(0, 1, 2, 3)
+        #     - traces.transpose(0, 1, 3, 2)
+        #     - traces.transpose(0, 3, 1, 2)
+        #     + traces.transpose(0, 3, 2, 1)
+        # )
+        # g_jikl = (
+        #     + traces.transpose(3, 0, 1, 2)
+        #     - traces.transpose(2, 0, 1, 3)
+        #     - traces.transpose(2, 0, 3, 1)
+        #     + traces.transpose(1, 0, 3, 2)
+        # )
         cumulant_function = - (
             + oe.contract('...kl,klji->...ij', decay_amplitudes, traces, backend='sparse')
             - oe.contract('...kl,kjli->...ij', decay_amplitudes, traces, backend='sparse')
             - oe.contract('...kl,kilj->...ij', decay_amplitudes, traces, backend='sparse')
             + oe.contract('...kl,kijl->...ij', decay_amplitudes, traces, backend='sparse')
         ) / 2
         if second_order:
+            # f_iklj = (
+            #     + traces.transpose(0, 1, 2, 3)
+            #     - traces.transpose(0, 2, 1, 3)
+            #     - traces.transpose(0, 2, 3, 1)
+            #     + traces.transpose(0, 3, 2, 1)
+            # )
+            #
+            # f_iklj = -g_iljk:
+            # f = -g.transpose(0, 2, 3, 1)
+            #
+            # f_jikl = (
+            #     + traces.transpose(3, 0, 1, 2)
+            #     - traces.transpose(3, 0, 2, 1)
+            #     - traces.transpose(1, 0, 2, 3)
+            #     + traces.transpose(1, 0, 3, 2)
+            # )
             cumulant_function -= (
                 + oe.contract('...kl,klji->...ij', frequency_shifts, traces, backend='sparse')
                 - oe.contract('...kl,lkji->...ij', frequency_shifts, traces, backend='sparse')
                 - oe.contract('...kl,klij->...ij', frequency_shifts, traces, backend='sparse')
                 + oe.contract('...kl,lkij->...ij', frequency_shifts, traces, backend='sparse')
             ) / 2
 
+    # Drop imaginary part since result is guaranteed to be real (if we didn't do anything wrong)
     return cumulant_function.real
 
 
@@ -1851,17 +1879,20 @@ def infidelity(
         The ``PulseSequence`` instance for which to calculate the
         infidelity for.
     spectrum: array_like, shape ([[n_nops,] n_nops,] omega) or callable
-        The two-sided noise power spectral density in units of inverse
-        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.
+        The noise power spectral density in units of inverse 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)``.
 
+        See :ref:`Notes <notes>` for a discussion on one- and two-sided
+        power spectral densities.
+
         If *test_convergence* is ``True``, a function handle to
         compute the power spectral density from a sequence of
         frequencies is expected.
@@ -1946,15 +1977,34 @@ def infidelity(
     infidelities that can be computed by setting
     ``which='correlations'``.
 
+    **One- and two-sided spectral densities**
 
-    To convert to the average gate infidelity, use the
-    following relation given by Horodecki et al. [Hor99]_ and
-    Nielsen [Nie02]_:
+    Since the real (imaginary) part of filter function :math:`F(\omega)`
+    is even (odd), it does not matter for integral whether
+    :math:`S(\omega)` is taken to be the one- or two-sided spectral
+    density. However, care should be taken that, if it is one or the
+    other, the frequencies :math:`\omega` are positive or symmetric
+    about zero, respectively.
+
+    To convert between one- and two-sided PSDs, use the following
+    relationship:
+
+    .. math::
+
+        S_\mathrm{onesided}(\omega) = 2 S_\mathrm{twosided}(\omega).
+
+    **Conversion to the Average Gate Infidelity (AGI)**
+
+    To convert the entanglement infidelity to the average gate
+    infidelity, use the following relation given by Horodecki et al.
+    [Hor99]_ and Nielsen [Nie02]_:
 
     .. math::
 
         \mathcal{I}_\mathrm{avg} = \frac{d}{d+1}\mathcal{I}.
 
+    **Goodness of approximation**
+
     The smallness parameter is given by
 
     .. math::

tests/test_basis.py

@@ -156,6 +156,25 @@ def test_basis_properties(self):
         self.assertTrue(basis.isorthonorm)
         self.assertArrayEqual(basis.T, basis.view(np.ndarray).T)
 
+    def test_transpose(self):
+        arr = rng.normal(size=(2, 3, 3))
+        b = arr.view(ff.Basis)
+        self.assertArrayEqual(b.T, arr.transpose(0, 2, 1))
+
+        arr = rng.normal(size=(2, 2))
+        b = arr.view(ff.Basis)
+        self.assertArrayEqual(b.T, arr.transpose(1, 0))
+
+        # Edge cases for not officially intended usage
+        arr = rng.normal(size=2)
+        b = arr.view(ff.Basis)
+        self.assertArrayEqual(b.T, arr)
+
+        arr = rng.normal(size=(2, 3, 4, 4))
+        b = arr.view(ff.Basis)
+        self.assertArrayEqual(b.T, arr.transpose(0, 1, 3, 2))
+
+
     def test_basis_expansion_and_normalization(self):
         """Correct expansion of operators and normalization of bases"""
         # dtype