Implement bases whose elements have finite trace (#11)

* Implement non-traceless bases

When dealing with different subspaces, one might want to generate an
operator basis that completely separates into elements that live
exclusively on the one and elements that live exclusively on the other
subspace. This is not possible when forcing all operator bases to have
an identity element (and thus the remaining elements traceless). The
core feature affected by this change is the way a complete basis is
generated from partial orthogonal elements. A switch now controls
whether a traceless or non-traceless is generated. Additionally, the
checks at the end of _full_from_partial() have been removed. The n
elements from which the basis is constructed now are the first n of
the final basis.

Including non-traceless bases also shone light on a bug in
calculate_error_transfer_matrix for a single qubit and the Pauli basis,
which did not include contributions from identity elements in case of a
cross-correlated spectrum. This induces non-unitality.

* Cleanup and comments

* Add more tests

* Drop imaginary part of integrands when purely real

* Fix bug for U.ndim > 2

* Add comments

* Add utf-8 coding hint

* Increase atol to prevent test from sometimes failing

* Get n_oper_identifiers from pulse if available

* Check same dims of control and noise Hamiltonian

* Add function to more quickly compute exp(1j*x)

This is a major time sink in the calculation of the control matrix using
calculate_control_matrix_from_scratch(). Computing real and imaginary
part separately and writing to a preallocated array should make the
whole calculation faster by a factor of two.

* Update notebooks due to deprecation warnings

* Make QFT example work again

* Update RB example to work with the current version

* Add show_progressbar switch to error_transfer_matrix()

* Remove test that didn't make sense

* Project out subspace identity element in CNOT infid

* Simplify oper_equiv

* Replace :meth: with more accurate :func: sphinx refs

* Catch non-traceless bases in infidelity()

In this case, additional terms need to be considered.

* Set correct cutoff frequencies and fix scaling in test_infidelity_cnot()

* Sort imports

Author: thangleiter

doc/source/examples/advanced_concatenation.ipynb

@@ -244,7 +244,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.4"
+   "version": "3.7.6"
   }
  },
  "nbformat": 4,

doc/source/examples/calculating_error_transfer_matrices.ipynb

@@ -21,7 +21,7 @@
    "outputs": [],
    "source": [
     "import numpy as np\n",
-    "import qutip as qt\n",
+    "from qutip.qip import operations\n",
     "\n",
     "import filter_functions as ff\n",
     "from filter_functions import util\n",
@@ -139,7 +139,7 @@
    "outputs": [],
    "source": [
     "cnot_12 = hadamard_2 @ cphase_12 @ hadamard_2\n",
-    "util.oper_equiv(cnot_12.total_Q, qt.cnot(control=0, target=1))"
+    "util.oper_equiv(cnot_12.total_Q, operations.cnot(control=0, target=1))"
    ]
   },
   {
@@ -168,7 +168,7 @@
     "\n",
     "cnot_21 = ff.remap(cnot_12, (1, 0), oper_identifier_mapping=mapping)\n",
     "\n",
-    "util.oper_equiv(cnot_21.total_Q, qt.cnot(control=1, target=0))\n",
+    "util.oper_equiv(cnot_21.total_Q, operations.cnot(control=1, target=0))\n",
     "print(cnot_21.c_oper_identifiers, cnot_21.n_oper_identifiers, sep='\\t')"
    ]
   },
@@ -230,7 +230,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.4"
+   "version": "3.7.6"
   }
  },
  "nbformat": 4,

doc/source/examples/extending_pulses.ipynb

@@ -252,6 +252,13 @@
    "source": [
     "The discrepancy to the analytical solution is due to the finite pulse width we have chosen here. We can see that, contrary to the FID-like filter function along $\\sigma_x$, the dephasing filter function of the SE pulse drops to (almost) zero for DC noise, affirming that indeed the pulse enhances coherence as it filters out frequencies slower than the pulse duration."
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {
@@ -270,7 +277,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.4"
+   "version": "3.7.6"
   }
  },
  "nbformat": 4,

doc/source/examples/getting_started.ipynb

@@ -138,7 +138,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.4"
+   "version": "3.7.6"
   }
  },
  "nbformat": 4,

doc/source/examples/periodic_driving.ipynb

@@ -29,25 +29,13 @@
     Hamiltonians: Example with ion traps, 1(1), 1–7.
     Retrieved from http://arxiv.org/abs/1503.08806
 """
-import filter_functions as ff
 import numpy as np
-import qutip as qt
-# %% Define some functions
-
-
-def get_a_k(k, n, N: int = 4):
-    if k < n:
-        a_k = 0
-    elif k == n:
-        a_k = 1/np.sqrt(2**(N - n + 4))
-    else:
-        a_k = 2**(n - k - 1)/get_a_k(n, n, N)
-
-    return a_k
 
+import filter_functions as ff
+import qutip as qt
+from qutip import qip
 
-def get_J_kl(k, l, n, N: int = 4):
-    return -np.pi*get_a_k(k, n, N)*get_a_k(l, n, N)/2
+# %% Define some functions
 
 
 def R_k_pulse(k, theta, phi, N: int = 4, tau: float = 1):
@@ -63,10 +51,10 @@ def R_k_pulse(k, theta, phi, N: int = 4, tau: float = 1):
         X = qt.tensor(X)
         Y = qt.tensor(Y)
 
-    H_c = [[X*np.cos(phi), [theta/2/tau]],
-           [Y*np.sin(phi), [theta/2/tau]]]
-    H_n = [[X/np.sqrt(X.shape[0]), [np.sqrt(X.shape[0])]],
-           [Y/np.sqrt(Y.shape[0]), [np.sqrt(Y.shape[0])]]]
+    H_c = [[X, [theta/2/tau*np.cos(phi)], 'I'*k + 'X' + 'I'*(N - k - 1)],
+           [Y, [theta/2/tau*np.sin(phi)], 'I'*k + 'Y' + 'I'*(N - k - 1)]]
+    H_n = [[X/np.sqrt(X.shape[0]), [1], 'I'*k + 'X' + 'I'*(N - k - 1)],
+           [Y/np.sqrt(Y.shape[0]), [1], 'I'*k + 'Y' + 'I'*(N - k - 1)]]
     dt = [tau]
 
     return ff.PulseSequence(H_c, H_n, dt)
@@ -76,17 +64,19 @@ def T_I_pulse(N: int = 4, tau: float = 1):
     Id = qt.qeye(2)
     if N == 1:
         T = Id
-        H_c = [[Id, [0]]]
-        H_n = [[T/np.sqrt(T.shape[0]), [np.sqrt(T.shape[0])]]]
+        H_c = [[Id, [0], 'I']]
+        H_n = [[T/np.sqrt(T.shape[0]), [1], 'I']]
     else:
         H_c = []
         H_n = []
         T = [Id]*(N - 1)
         T.insert(0, qt.sigmaz())
         for k in range(1, N+1):
-            H_c.append([qt.tensor(np.roll(T, k-1)),
-                        [np.pi/4*(1 - 2**(1 - k))/tau]])
-            H_n.append([qt.tensor(T)/np.sqrt(2**len(T)), [np.sqrt(2**len(T))]])
+            H_c.append([qt.tensor(T[-k+1:] + T[:-k+1]),
+                        [np.pi/4*(1 - 2**(1 - k))/tau],
+                        'I'*(k - 1) + 'Z' + 'I'*(N - k)])
+            H_n.append([qt.tensor(T[-k+1:] + T[:-k+1])/np.sqrt(2**len(T)),
+                        [1], 'I'*(k - 1) + 'Z' + 'I'*(N - k)])
 
     return ff.PulseSequence(H_c, H_n, [tau])
 
@@ -95,38 +85,35 @@ def T_F_pulse(N: int = 4, tau: float = 1):
     Id = qt.qeye(2)
     if N == 1:
         T = Id
-        H_c = [[Id, [0]]]
-        H_n = [[T/np.sqrt(T.shape[0]), [np.sqrt(T.shape[0])]]]
+        H_c = [[Id, [0], 'I']]
+        H_n = [[T/np.sqrt(T.shape[0]), [1], 'I']]
     else:
         H_c = []
         H_n = []
         T = [Id]*(N - 1)
         T.insert(0, qt.sigmaz())
         for k in range(1, N+1):
-            H_c.append([qt.tensor(np.roll(T, k-1)),
-                        [np.pi/4*(1 - 2**(k - N))/tau]])
-            H_n.append([qt.tensor(T)/np.sqrt(2**len(T)), [np.sqrt(2**len(T))]])
+            H_c.append([qt.tensor(T[-k+1:] + T[:-k+1]),
+                        [np.pi/4*(1 - 2**(k - N))/tau],
+                        'I'*(k - 1) + 'Z' + 'I'*(N - k)])
+            H_n.append([qt.tensor(T[-k+1:] + T[:-k+1])/np.sqrt(2**len(T)),
+                        [1], 'I'*(k - 1) + 'Z' + 'I'*(N - k)])
 
     return ff.PulseSequence(H_c, H_n, [tau])
 
 
-def U_pulse(n, N: int = 4, tau: float = 1):
+def P_n_pulse(n, N: int = 4, tau: float = 1):
     Id = qt.qeye(2)
-    if N == 2:
-        Z = qt.tensor([qt.sigmaz()]*2)
-        H_c = [[Z, [get_J_kl(1, 2, 1, N)/tau]]]
-        H_n = [[Z/np.sqrt(Z.shape[0]), [np.sqrt(Z.shape[0])]]]
-    else:
-        H_c = []
-        H_n = []
-        for k in range(n, N+1):
-            for l in range(k+1, N+1):
-                Z = [Id]*(N - 2)
-                Z.insert(k-1, qt.sigmaz())
-                Z.insert(l-1, qt.sigmaz())
-                Z = qt.tensor(Z)
-                H_c.append([Z, [get_J_kl(k, l, n, N)/tau]])
-                H_n.append([Z/np.sqrt(Z.shape[0]), [np.sqrt(Z.shape[0])]])
+    H_c = []
+    H_n = []
+    for l in range(n+1, N+1):
+        Z = [Id]*(N - 2)
+        Z.insert(n-1, qt.sigmaz())
+        Z.insert(l-1, qt.sigmaz())
+        Z = qt.tensor(Z)
+        identifier = ('I'*(n-1) + 'Z' + 'I'*(l - n - 1) + 'Z' + 'I'*(N - l))
+        H_c.append([Z, [-np.pi/4*2**(n - l)/tau], identifier])
+        H_n.append([Z/np.sqrt(Z.shape[0]), [1], identifier])
 
     return ff.PulseSequence(H_c, H_n, [tau])
 
@@ -140,7 +127,7 @@ def QFT_pulse(N: int = 4, tau: float = 1):
     pulses = [T_I_pulse(N, tau)]
     for n in range(N-1):
         pulses.append(H_k_pulse(n, N, tau))
-        pulses.append(U_pulse(n+1, N, tau))
+        pulses.append(P_n_pulse(n+1, N, tau))
 
     pulses.append(H_k_pulse(N-1, N, tau))
     pulses.append(T_F_pulse(N, tau))
@@ -150,7 +137,18 @@ def QFT_pulse(N: int = 4, tau: float = 1):
 
 # %% Get a 4-qubit QFT and plot the filter function
 omega = np.logspace(-2, 2, 500)
-QFT = QFT_pulse(4)
+
+N = 4
+QFT = QFT_pulse(N)
+
+# Check the pulse produces the correct propagator after swapping the qubits
+swaps = [qip.operations.swap(N, [i, j]).full()
+         for i, j in zip(range(N//2), range(N-1, N//2-1, -1))]
+prop = ff.util.mdot(swaps) @ QFT.total_Q
+qt.matrix_histogram_complex(prop)
+print('Correct action: ',
+      ff.util.oper_equiv(prop, qip.algorithms.qft.qft(N), eps=1e-14))
+
 fig, ax, _ = ff.plot_filter_function(QFT, omega)
 # Move the legend to the side because of many entries
-ax.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
+ax.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)

doc/source/examples/quantum_fourier_transform.ipynb

@@ -97,17 +97,17 @@ def run_randomized_benchmarking(N_G: int, N_l: int, min_l: int, max_l: int,
 H_n = {}
 H_c = {}
 dt = {}
-H_c['Id'] = [[qt.sigmax().full(), [0]]]
-H_n['Id'] = [[qt.sigmax().full(), [1]],
-             [qt.sigmay().full(), [1]]]
+H_c['Id'] = [[qt.sigmax().full(), [0], 'X']]
+H_n['Id'] = [[qt.sigmax().full(), [1], 'X'],
+             [qt.sigmay().full(), [1], 'Y']]
 dt['Id'] = [T]
-H_c['X2'] = [[qt.sigmax().full(), [np.pi/4/T]]]
-H_n['X2'] = [[qt.sigmax().full(), [1]],
-             [qt.sigmay().full(), [1]]]
+H_c['X2'] = [[qt.sigmax().full(), [np.pi/4/T], 'X']]
+H_n['X2'] = [[qt.sigmax().full(), [1], 'X'],
+             [qt.sigmay().full(), [1], 'Y']]
 dt['X2'] = [T]
-H_c['Y2'] = [[qt.sigmay().full(), [np.pi/4/T]]]
-H_n['Y2'] = [[qt.sigmax().full(), [1]],
-             [qt.sigmay().full(), [1]]]
+H_c['Y2'] = [[qt.sigmay().full(), [np.pi/4/T], 'Y']]
+H_n['Y2'] = [[qt.sigmax().full(), [1], 'X'],
+             [qt.sigmay().full(), [1], 'Y']]
 dt['Y2'] = [T]
 
 # %% Set up PulseSequences
@@ -169,7 +169,7 @@ def run_randomized_benchmarking(N_G: int, N_l: int, min_l: int, max_l: int,
 # Scaling factor for the noise so that alpha = 0 and alpha = 0.7 have the same
 # power
 noise_scaling_factor = {
-    0.0: 0.39651222009884685,
+    0.0: 0.4415924985735799,
     0.7: 1
 }
 
@@ -200,16 +200,16 @@ def run_randomized_benchmarking(N_G: int, N_l: int, min_l: int, max_l: int,
 for i, alpha in enumerate((0.0, 0.7)):
 
     means = np.mean(fidelities[alpha], axis=1)
-    errs = np.std(fidelities[alpha], axis=1)/np.sqrt(N_G)
+    stds = np.std(fidelities[alpha], axis=1)
 
-    popt, pcov = optimize.curve_fit(fitfun, lengths, means, [0, 1], errs,
+    popt, pcov = optimize.curve_fit(fitfun, lengths, means, [0, 1], stds,
                                     absolute_sigma=True)
 
     for j in range(N_G):
         fid = ax[i].plot(lengths, fidelities[alpha][:, j], 'k.', alpha=0.1,
                          zorder=2)
 
-    mean = ax[i].errorbar(lengths, means, errs, fmt='.', zorder=3,
+    mean = ax[i].errorbar(lengths, means, stds, fmt='.', zorder=3,
                           color='tab:red')
     fit = ax[i].plot(lengths, fitfun(lengths, *popt), '--', zorder=4,
                      color='tab:red')

doc/source/examples/qutip_integration.ipynb

@@ -1,3 +1,4 @@
+# -*- coding: utf-8 -*-
 # =============================================================================
 #     filter_functions
 #     Copyright (C) 2019 Quantum Technology Group, RWTH Aachen University
@@ -21,21 +22,21 @@
 
 from . import analytic, basis, numeric, plotting, pulse_sequence, util
 from .basis import Basis
-from .numeric import (infidelity, liouville_representation,
-                      error_transfer_matrix)
-from .plotting import (plot_bloch_vector_evolution, plot_filter_function,
-                       plot_pulse_correlation_filter_function,
-                       plot_pulse_train, plot_error_transfer_matrix)
+from .numeric import (error_transfer_matrix, infidelity,
+                      liouville_representation)
+from .plotting import (
+    plot_bloch_vector_evolution, plot_error_transfer_matrix,
+    plot_filter_function, plot_pulse_correlation_filter_function,
+    plot_pulse_train)
 from .pulse_sequence import (PulseSequence, concatenate, concatenate_periodic,
                              extend, remap)
 
-__all__ = ['PulseSequence', 'concatenate', 'extend', 'Basis', 'numeric',
-           'plot_filter_function', 'plot_bloch_vector_evolution',
-           'plot_pulse_train', 'util', 'plotting', 'analytic', 'basis',
-           'infidelity', 'plot_pulse_correlation_filter_function',
-           'concatenate_periodic', 'pulse_sequence', 'remap',
-           'error_transfer_matrix', 'plot_error_transfer_matrix',
-           'liouville_representation']
+__all__ = ['Basis', 'PulseSequence', 'analytic', 'basis', 'concatenate',
+           'concatenate_periodic', 'error_transfer_matrix', 'extend',
+           'infidelity', 'liouville_representation', 'numeric',
+           'plot_bloch_vector_evolution', 'plot_error_transfer_matrix',
+           'plot_filter_function', 'plot_pulse_correlation_filter_function',
+           'plot_pulse_train', 'plotting', 'pulse_sequence', 'remap', 'util']
 
 __version__ = '0.1.0'
 __license__ = 'GNU GPLv3+'