Merge 2cee556 into 2578a34

README.md

@@ -4,11 +4,59 @@
 [![Documentation Status](https://readthedocs.org/projects/filter-functions/badge/?version=latest)](https://filter-functions.readthedocs.io/en/latest/?badge=latest)
 [![PyPI version](https://badge.fury.io/py/filter-functions.svg)](https://badge.fury.io/py/filter-functions)
 
-Simply put, filter functions characterize a pulse's susceptibility to noise at a given frequency and can thus be used to gain insight into the decoherence of the system. The formalism allows for efficient calculation of several quantities of interest such as average gate fidelity. Moreover, the filter function of a composite pulse can be easily derived from those of the constituent pulses, allowing for efficient assembly and characterization of pulse sequences.
+## Introduction
+Simply put, filter functions characterize a quantum system's susceptibility to noise at a given frequency during a control operation and can thus be used to gain insight into its decoherence. The formalism allows for efficient calculation of several quantities of interest such as average gate fidelity and even the entire quantum process up to a unitary rotation. Moreover, the filter function of a composite pulse can be easily derived from those of the constituent pulses, allowing for efficient assembly and characterization of pulse sequences.
 
-Previously, filter functions have only been computed analytically for select pulses such as dynamical decoupling sequences [1], [2]. With this project we aim to provide a toolkit for calculating and inspecting filter functions for arbitrary pulses including pulses without analytic form such as one might get from numerical pulse optimization algorithms. 
+Previously, filter functions have only been computed analytically for select pulses such as dynamical decoupling sequences [1, 2]. With this project we aim to provide a toolkit for calculating and inspecting filter functions for arbitrary pulses including pulses without analytic form such as one might get from numerical pulse optimization algorithms.
 
-The package is built to interface with [QuTiP](http://qutip.org/), a widely-used quantum toolbox for Python, and comes with extensive documentation and a test suite.
+The `filter_functions` package is built to interface with [QuTiP](http://qutip.org/), a widely-used quantum toolbox for Python, and comes with extensive documentation and a test suite. Note that the project is still in pre-release and thus liable to breaking API changes.
+
+As a very brief introduction, consider a Hadamard gate implemented by a pi/2 Y-gate followed by a NOT-gate using simple square pulses. We can calculate and plot the dephasing filter function of the gate with the following code:
+
+```python
+import filter_functions as ff
+import qutip as qt
+from math import pi
+
+H_c = [[qt.sigmax()/2,   [0, pi], 'X'],
+       [qt.sigmay()/2, [pi/2, 0], 'Y']]     # control Hamiltonian
+H_n = [[qt.sigmaz()/2,    [1, 1], 'Z']]     # constant coupling to z noise
+dt = [1, 1]                                 # time steps
+
+hadamard = ff.PulseSequence(H_c, H_n, dt)   # Central object representing a control pulse
+omega = ff.util.get_sample_frequencies(hadamard)
+F = hadamard.get_filter_function(omega)
+
+ff.plot_filter_function(hadamard)           # Filter function cached from before
+```
+
+![Hadamard dephasing filter function](./doc/source/_static/hadamard.png)
+
+An alternative way of obtaining the Hadamard `PulseSequence` is by concatenating the composing pulses:
+
+```python
+Y2 = ff.PulseSequence([[qt.sigmay()/2, [pi/2], 'Y']],
+                      [[qt.sigmaz()/2,    [1], 'Z']],
+                      [1])
+X = ff.PulseSequence([[qt.sigmax()/2, [pi], 'X']],
+                     [[qt.sigmaz()/2,  [1], 'Z']],
+                     [1])
+
+Y2.cache_filter_function(omega)
+X.cache_filter_function(omega)
+
+hadamard = Y2 @ X           # equivalent: ff.concatenate((Y2, X))
+hadamard.is_cached('F')
+# True  (filter function cached during concatenation)
+```
+
+To compute, for example, the infidelity of the gate in the presence of an arbitrary classical noise spectrum, we can simply call `infidelity()`:
+
+```python
+spectrum = 1e-2/abs(omega)  # omega is symmetric about zero
+infidelity = ff.infidelity(hadamard, spectrum, omega)
+# array([0.006037])  (one contribution per noise operator)
+```
 
 ## Installation
 To install the package from PyPI, run `pip install filter_functions`. It is recommended to install QuTiP before by following the [instructions on their website](http://qutip.org/docs/latest/installation.html) rather than installing it through `pip`. To install the package from source run `python setup.py develop` to install using symlinks or `python setup.py install` without.
@@ -25,4 +73,4 @@ Building the documentation requires the following additional dependencies: `nbsp
 ## References
 [1]: Cywinski, L., Lutchyn, R. M., Nave, C. P., & Das Sarma, S. (2008). How to enhance dephasing time in superconducting qubits. Physical Review B - Condensed Matter and Materials Physics, 77(17), 1–11. [https://doi.org/10.1103/PhysRevB.77.174509](https://doi.org/10.1103/PhysRevB.77.174509)
 
-[2]: Green, T. J., Sastrawan, J., Uys, H., & Biercuk, M. J. (2013). Arbitrary quantum control of qubits in the presence of universal noise. New Journal of Physics, 15(9), 095004. [https://doi.org/10.1088/1367-2630/15/9/095004](https://doi.org/10.1088/1367-2630/15/9/095004)
+[2]: Green, T. J., Sastrawan, J., Uys, H., & Biercuk, M. J. (2013). Arbitrary quantum control of qubits in the presence of universal noise. New Journal of Physics, 15(9), 095004. [https://doi.org/10.1088/1367-2630/15/9/095004](https://doi.org/10.1088/1367-2630/15/9/095004)

doc/source/conf.py

@@ -49,6 +49,7 @@
     'sphinx.ext.viewcode',
     'sphinx.ext.ifconfig',
     'sphinx.ext.napoleon',
+    'sphinx.ext.intersphinx',
     #'sphinxcontrib.apidoc',
     #'IPython.sphinxext.ipython_console_highlighting',
     #'IPython.sphinxext.ipython_directive',
@@ -274,11 +275,11 @@
 
 # Example configuration for intersphinx: refer to the Python standard library.
 intersphinx_mapping = {
-    'python': ('https://docs.python.org/dev', None),
-    'numpy': ('https://docs.scipy.org/doc/numpy', None),
-    'scipy': ('https://docs.scipy.org/doc/scipy/reference', None),
-    'matplotlib': ('https://matplotlib.org', None),
-    'qutip': ('https://qutip.org/docs/latest', None)
+    'python': ('https://docs.python.org/3/', None),
+    'numpy': ('https://numpy.org/doc/stable/', None),
+    'scipy': ('https://docs.scipy.org/doc/scipy/reference/', None),
+    'matplotlib': ('https://matplotlib.org/', None),
+    'qutip': ('http://qutip.org/docs/latest/', None)
 }
 
 # -- Options for todo extension ----------------------------------------------

filter_functions/basis.py

@@ -89,7 +89,6 @@ class Basis(ndarray):
 
     Parameters
     ----------
-
     basis_array : array_like, shape (n, d, d)
         An array or list of square matrices that are elements of an operator
         basis spanning :math:`\mathbb{C}^{d\times d}`. *n* should be smaller
@@ -142,10 +141,10 @@ class Basis(ndarray):
     Other than the methods inherited from ``ndarray``, a ``Basis`` instance has
     the following methods:
 
-    :meth:`normalize`
+    normalize(b)
         Normalizes the basis in-place (used internally when creating a basis
         from elements)
-    :meth:`tidyup`
+    tidyup(eps_scale=None)
         Cleans up floating point errors in-place to make zeros actual zeros.
         ``eps_scale`` is an optional argument multiplied to the data type's
         ``eps`` to get the absolute tolerance.
@@ -614,7 +613,8 @@ def expand(M: Union[ndarray, Basis], basis: Union[ndarray, Basis],
 
     .. math::
         M &= \sum_j c_j C_j, \\
-        c_j &= \frac{\mathrm{tr}\big(M C_j\big)}{\mathrm{tr}\big(C_j^2\big)}.
+        c_j &= \frac{\mathrm{tr}\big(M C_j\big)}
+                    {\mathrm{tr}\big(C_j^\dagger C_j\big)}.
 
     """
     coefficients = np.einsum('...ij,bji->...b', np.asarray(M), basis)

filter_functions/numeric.py

@@ -74,7 +74,7 @@
 
 
 def calculate_control_matrix_from_atomic(
-        phases: ndarray, R_l: ndarray, Q_liouville: ndarray,
+        phases: ndarray, R_g: ndarray, Q_liouville: ndarray,
         show_progressbar: Optional[bool] = None) -> ndarray:
     r"""
     Calculate the control matrix from the control matrices of atomic segments.
@@ -83,7 +83,7 @@ def calculate_control_matrix_from_atomic(
     ----------
     phases : array_like, shape (n_dt, n_omega)
         The phase factors for :math:`l\in\{0, 1, \dots, n-1\}`.
-    R_l : array_like, shape (n_dt, n_nops, d**2, n_omega)
+    R_g : array_like, shape (n_dt, n_nops, d**2, n_omega)
         The pulse control matrices for :math:`l\in\{1, 2, \dots, n\}`.
     Q_liouville : array_like, shape (n_dt, n_nops, d**2, d**2)
         The transfer matrices of the cumulative propagators for
@@ -107,24 +107,22 @@ def calculate_control_matrix_from_atomic(
 
     See Also
     --------
-    :func:`calculate_control_matrix_from_scratch`
-
-    :func:`liouville_representation`
+    calculate_control_matrix_from_scratch : Control matrix from scratch.
+    liouville_representation : Liouville representation for a given basis.
     """
-    n = len(R_l)
+    n = len(R_g)
     # Allocate memory
-    R = np.zeros(R_l[0].shape, dtype=complex)
+    R = np.zeros(R_g.shape[1:], dtype=complex)
 
     # Set up a reusable contraction expression. In some cases it is faster to
     # also contract the time dimension in the same expression instead of
     # looping over it, but we don't distinguish here for readability.
-    R_expr = contract_expression('o,ijo,jk->iko', phases[0].shape,
-                                 R_l[0].shape, Q_liouville[0].shape,
-                                 optimize=[(0, 1), (0, 1)])
+    R_expr = contract_expression('ijo,jk->iko',
+                                 R_g.shape[1:], Q_liouville.shape[1:])
 
-    for l in progressbar_range(n, show_progressbar=show_progressbar,
+    for g in progressbar_range(n, show_progressbar=show_progressbar,
                                desc='Calculating control matrix'):
-        R += R_expr(phases[l], R_l[l], Q_liouville[l])
+        R += R_expr(phases[g]*R_g[g], Q_liouville[g])
 
     return R
 
@@ -210,20 +208,16 @@ def calculate_control_matrix_from_scratch(
 
     See Also
     --------
-    :func:`calculate_control_matrix_from_atomic`
+    calculate_control_matrix_from_atomic : Control matrix from concatenation.
     """
     if t is None:
         t = np.concatenate(([0], np.asarray(dt).cumsum()))
 
     d = HV.shape[-1]
-    n = len(dt)
     # We're lazy
     E = omega
     n_coeffs = np.asarray(n_coeffs)
 
-    # Allocate memory
-    R = np.zeros((len(n_opers), len(basis), len(E)), dtype=complex)
-
     # Precompute noise opers transformed to eigenbasis of each pulse
     # segment and Q^\dagger @ HV
     if d < 4:
@@ -233,39 +227,35 @@ def calculate_control_matrix_from_scratch(
                       optimize=['einsum_path', (0, 1), (0, 1)])
     else:
         QdagV = Q[:-1].transpose(0, 2, 1).conj() @ HV
-        B = np.empty((len(n_opers), n, d, d), dtype=complex)
+        B = np.empty((len(n_opers), len(dt), d, d), dtype=complex)
         for j, n_oper in enumerate(n_opers):
             B[j] = HV.conj().transpose(0, 2, 1) @ n_oper @ HV
 
-    # Allocate array for the integral
-    integral = np.empty((len(E), d, d), dtype=complex)
+    # Allocate result and buffers for intermediate arrays
+    R = np.zeros((len(n_opers), len(basis), len(E)), dtype=complex)
+    exp_buf = np.empty((len(E), d, d), dtype=complex)
+    int_buf = np.empty((len(E), d, d), dtype=complex)
     R_path = ['einsum_path', (0, 3), (0, 1), (0, 2), (0, 1)]
 
-    for l in progressbar_range(n, show_progressbar=show_progressbar,
+    for l in progressbar_range(len(dt), show_progressbar=show_progressbar,
                                desc='Calculating control matrix'):
-        # Create a (n_E, d, d)-shaped array containing the energy
-        # differences in its last two dimensions
+
         dE = np.subtract.outer(HD[l], HD[l])
-        # Add the frequencies to get EdE_nm = omega + omega_n - omega_m
-        EdE = np.add.outer(E, dE)
-
-        # Mask the integral to avoid convergence problems
-        mask_small = np.abs(EdE*dt[l]) <= 1e-7
-        integral[mask_small] = dt[l]
-        integral[~mask_small] = (cexp(EdE[~mask_small]*dt[l]) - 1) /\
-                                (1j*(EdE[~mask_small]))
-        """
-        # Test convergence of integral as argument of exponential -> 0
-        fn = lambda a, dt: (np.exp(a*dt) - 1)/a - dt
-        fn1 = lambda a: fn(a, 1)
-        a = np.logspace(-10, 0, 100)
-        plt.loglog(a, fn1(a))
-        """
+        # iEdE_nm = 1j*(omega + omega_n - omega_m)
+        int_buf.real = 0
+        int_buf.imag = np.add.outer(E, dE, out=int_buf.imag)
+
+        # Use expm1 for better convergence with small arguments
+        exp_buf = np.expm1(int_buf*dt[l], out=exp_buf)
+        # Catch zero-division warnings
+        mask = (int_buf != 0)
+        int_buf = np.divide(exp_buf, int_buf, out=int_buf, where=mask)
+        int_buf[~mask] = dt[l]
 
         # Faster for d = 2 to also contract over the time dimension instead of
         # loop, but for readability we don't distinguish.
         R += np.einsum('o,j,jmn,omn,knm->jko',
-                       cexp(E*t[l]), n_coeffs[:, l], B[:, l], integral,
+                       cexp(E*t[l]), n_coeffs[:, l], B[:, l], int_buf,
                        QdagV[l].conj().T @ basis @ QdagV[l],
                        optimize=R_path)
 
@@ -446,11 +436,9 @@ def calculate_filter_function(R: ndarray) -> ndarray:
 
     See Also
     --------
-    :func:`calculate_control_matrix_from_scratch`
-
-    :func:`calculate_control_matrix_from_atomic`
-
-    :func:`calculate_pulse_correlation_filter_function`
+    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)
 
@@ -486,11 +474,9 @@ def calculate_pulse_correlation_filter_function(R: ndarray) -> ndarray:
 
     See Also
     --------
-    :func:`calculate_control_matrix_from_scratch`
-
-    :func:`calculate_control_matrix_from_atomic`
-
-    :func:`calculate_filter_function`
+    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)
@@ -682,9 +668,8 @@ def error_transfer_matrix(
 
     See Also
     --------
-    :func:`calculate_error_vector_correlation_functions`
-
-    :func:`infidelity`
+    calculate_error_vector_correlation_functions
+    infidelity : Calculate only infidelity of a pulse.
     """
     N, d = pulse.basis.shape[:2]
     u_kl = calculate_error_vector_correlation_functions(pulse, S, omega,
@@ -773,12 +758,16 @@ def infidelity(pulse: 'PulseSequence',
         into account.
     which : str, optional
         Which infidelities should be calculated, may be either 'total'
-        (default) or 'correlations'. In the former case, only the total
-        infidelities for each noise operator are returned, in the latter all of
-        the individual infidelity contributions including the pulse
-        correlations (note that in this case no checks are performed if the
-        frequencies are compliant). See :func:`~pulse_sequence.concatenate`
-        for more details.
+        (default) or 'correlations'. In the former case, one value is returned
+        for each noise operator, corresponding to the total infidelity of the
+        pulse (or pulse sequence). In the latter, an array of infidelities is
+        returned where element (i,j) corresponds to the infidelity contribution
+        of the correlations between pulses i and j (see :ref:`Notes <notes>`).
+        Note that this option is only available if the pulse correlation filter
+        functions have been computed during concatenation (see
+        :func:`calculate_pulse_correlation_filter_function` and
+        :func:`~filter_functions.pulse_sequence.concatenate`) and that in this
+        case no checks are performed if the frequencies are compliant.
     return_smallness : bool, optional
         Return the smallness parameter :math:`\xi` for the given spectrum.
     test_convergence : bool, optional
@@ -805,26 +794,30 @@ def infidelity(pulse: 'PulseSequence',
         The matplotlib figure and axis instances used for plotting.
         Only if *test_convergence* is ``True``.
 
+    .. _notes:
+
     Notes
     -----
     The infidelity is given by
 
     .. math::
 
-        \big\langle\mathcal{I}_\mathrm{e}\big\rangle &=
-                \frac{1}{d^2}\left\langle
-                \mathrm{tr}\big\lvert\tilde{U}(\tau)\big\rvert^2
-                \right\rangle \\
-            &= \frac{1}{2\pi d}\int_{-\infty}^{\infty}\mathrm{d}\omega
-                \,\mathrm{tr}\bigl(S(\omega)F(\omega)\bigr) +
-                \mathcal{O}\big(\xi^4\big)
+        \big\langle\mathcal{I}_\mathrm{e}\big\rangle_{\alpha\beta} &=
+                \frac{1}{2\pi d}\int_{-\infty}^{\infty}\mathrm{d}\omega\,
+                S_{\alpha\beta}(\omega)F_{\alpha\beta}(\omega)
+                +\mathcal{O}\big(\xi^4\big) \\
+            &= \sum_{g,g'=1}^G \big\langle
+                \mathcal{I}_\mathrm{e}\big\rangle_{\alpha\beta}^{(gg')}
 
     with :math:`S_{\alpha\beta}(\omega)` the two-sided noise spectral density
     and :math:`F_{\alpha\beta}(\omega)` the first-order filter function for
     noise sources :math:`\alpha,\beta`. The noise spectrum may include
     correlated noise sources, that is, its entry at :math:`(\alpha,\beta)`
     corresponds to the correlations between sources :math:`\alpha` and
     :math:`\beta`.
+    :math:`\big\langle\mathcal{I}_\mathrm{e}\big\rangle_{\alpha\beta}^{(gg')}`
+    are the correlation infidelities that can be computed by setting
+    ``which='correlations'``.
 
     To convert to the average gate infidelity, use the
     following relation given by Horodecki et al. [Hor99]_ and

filter_functions/plotting.py

@@ -145,7 +145,7 @@ def plot_bloch_vector_evolution(pulse: 'PulseSequence',
     show** : bool, optional
         Whether to show the sphere (by calling :code:`b.make_sphere()`).
     return_Bloch : bool, optional
-        Whether to return the ``qutip.Bloch`` instance
+        Whether to return the :class:`qutip.bloch.Bloch` instance
     bloch_kwargs : dict, optional
         A dictionary with keyword arguments to be fed into the Bloch
         constructor (if *b* not given).
@@ -162,7 +162,8 @@ def plot_bloch_vector_evolution(pulse: 'PulseSequence',
 
     See Also
     --------
-    :class:`qutip.Bloch`
+    qutip.bloch.Bloch : Qutip's Bloch sphere implementation.
+    scipy.linalg.expm : asdf
     """
     # Raise an exception if not a one-qubit pulse
     if not pulse.d == 2: