Merge ba3b685 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

@@ -45,13 +45,13 @@
 
 import numpy as np
 import opt_einsum as oe
+from numpy import linalg as nla
 from numpy.core import ndarray
-from numpy.linalg import norm
 from qutip import Qobj
-from scipy.linalg import null_space
+from scipy import linalg as sla
 from sparse import COO
 
-from .util import P_np, remove_float_errors, tensor
+from . import util
 
 __all__ = ['Basis', 'expand', 'ggm_expand', 'normalize']
 
@@ -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.
@@ -285,7 +284,7 @@ def istraceless(self) -> bool:
         """
         if self._istraceless is None:
             trace = np.einsum('...jj', self)
-            trace = remove_float_errors(trace, self.d)
+            trace = util.remove_float_errors(trace, self.d)
             nonzero = trace.nonzero()
             if nonzero[0].size == 0:
                 self._istraceless = True
@@ -376,9 +375,9 @@ def four_element_traces(self, traces):
     def normalize(self) -> None:
         """Normalize the basis in-place"""
         if self.ndim == 2:
-            self /= norm(self)
+            self /= nla.norm(self)
         elif self.ndim == 3:
-            np.einsum('ijk,i->ijk', self, 1/norm(self, axis=(1, 2)),
+            np.einsum('ijk,i->ijk', self, 1/nla.norm(self, axis=(1, 2)),
                       out=self)
 
     def tidyup(self, eps_scale: Optional[float] = None) -> None:
@@ -425,7 +424,8 @@ def pauli(cls, n: int) -> 'Basis':
         """
         normalization = np.sqrt(2**n)
         combinations = np.indices((4,)*n).reshape(n, 4**n)
-        sigma = tensor(*np.array(P_np)[combinations], rank=2)/normalization
+        sigma = util.tensor(*np.array(util.P_np)[combinations], rank=2)
+        sigma /= normalization
         return cls(sigma, btype='Pauli', skip_check=True)
 
     @classmethod
@@ -535,7 +535,7 @@ def _full_from_partial(elems: Sequence, traceless: Union[None, bool]) -> Basis:
         # Those together with coeffs span the whole space, and therefore also
         # the linear combinations of GGMs weighted with the coefficients will
         # span the whole matrix space
-        coeffs = np.concatenate((coeffs, null_space(coeffs).T))
+        coeffs = np.concatenate((coeffs, sla.null_space(coeffs).T))
         # Our new basis is given by linear combinations of GGMs with coeffs
         basis = np.einsum('ij,jkl', coeffs, ggm)
     else:
@@ -573,9 +573,10 @@ def normalize(b: Sequence) -> Basis:
     """
     b = np.asanyarray(b)
     if b.ndim == 2:
-        return (b/norm(b)).view(Basis)
+        return (b/nla.norm(b)).view(Basis)
     if b.ndim == 3:
-        return np.einsum('ijk,i->ijk', b, 1/norm(b, axis=(1, 2))).view(Basis)
+        return np.einsum('ijk,i->ijk',
+                         b, 1/nla.norm(b, axis=(1, 2))).view(Basis)
 
     raise ValueError('Expected b.ndim to be either 2 or 3, not ' +
                      '{}.'.format(b.ndim))
@@ -614,15 +615,16 @@ 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)
 
     if not normalized:
         coefficients /= np.einsum('bij,bji->b', basis, basis).real
 
-    return remove_float_errors(coefficients) if tidyup else coefficients
+    return util.remove_float_errors(coefficients) if tidyup else coefficients
 
 
 def ggm_expand(M: Union[ndarray, Basis], traceless: bool = False) -> ndarray: