Commit 62c2dbcb authored by Qutech lab's avatar Qutech lab
Browse files

Merge remote-tracking branch 'origin/master' into selene/master

parents e52f65df c28fddc6
......@@ -39,6 +39,9 @@ for i in qutil.ui.progressbar(range(n)):
## qutil.qi
In this module there are some quantities and functions related to quantum information, like the Pauli matrices in different data types.
## qutil.random
Here we collect functions for random numbers like `random_hermitian` to generate random Hermitian matrices.
## qutil.itertools
This module contains a everything from `itertools`, `more_itertools` and custom functions.
......
from . import const, linalg, matlab, plotting, ui, qi, io
from . import const, linalg, matlab, plotting, ui, qi, random, io
__version__ = '0.1'
__all__ = ['const', 'linalg', 'matlab', 'ui', 'plotting', 'qi', 'io']
__all__ = ['const', 'linalg', 'matlab', 'ui', 'plotting', 'qi', 'io', 'random']
import sys
import csv
import os.path
import pickle
def show_saved_figure(path):
"""Unpickle and show a pickled matplotlib figure.
Parameters
----------
path: str | pathlib.Path
The path where the figure is saved.
Returns
-------
fig: matplotlib.figure.Figure
The matplotlib figure instance.
"""
fig = pickle.load(open(path, 'rb'))
import matplotlib.pyplot as plt
dummy = plt.figure()
new_manager = dummy.canvas.manager
new_manager.canvas.figure = fig
fig.set_canvas(new_manager.canvas)
plt.show()
return fig
def save_figure(fig, path):
"""Pickle a matplotlib figure.
Parameters
----------
fig: matplotlib.figure.Figure
The matplotlib figure instance to be saved.
path: str | pathlib.Path
The path where the figure is saved.
"""
pickle.dump(fig, open(path, 'wb'))
def query_yes_no(question, default="yes"):
"""Ask a yes/no question via input() and return their answer.
......@@ -52,11 +94,11 @@ class CsvLogger:
if reader.fieldnames != fieldnames:
raise RuntimeError("Existing file has differing fieldnames", reader.fieldnames, fieldnames)
else:
with open(self.filename, 'x') as file:
csv.DictWriter(file, fieldnames=self.fieldnames, dialect=self.dialect).writeheader()
def write(self, *args):
with open(self.filename, 'a+') as file:
csv.DictWriter(file, fieldnames=self.fieldnames, dialect=self.dialect).writerow(dict(zip(self.fieldnames, args)))
......@@ -5,6 +5,10 @@ Functions
---------
:meth:`abs2` :
Fast absolute value squared
:meth:`cexp` :
Fast complex exponential
:meth:`closest_unitary` :
Get the closest unitary on a subblock of a matrix
:meth:`tensor`
Fast, flexible tensor product of an arbitrary number of inputs using
:meth:`~numpy.einsum` that supports broadcasting
......@@ -33,6 +37,7 @@ from typing import List, Sequence, Tuple, Union
import numpy as np
from numpy import linalg, ndarray
from scipy.linalg import polar
from .qi import paulis
......@@ -42,8 +47,8 @@ try:
except ImportError:
numba = False
__all__ = ['abs2', 'check_phase_eq', 'density', 'dot_HS', 'max_abs_diff',
'max_rel_diff', 'mdot', 'pauli_expm', 'ptrace',
__all__ = ['abs2', 'cexp', 'check_phase_eq', 'density', 'dot_HS',
'max_abs_diff', 'max_rel_diff', 'mdot', 'pauli_expm', 'ptrace',
'remove_float_errors', 'sparsity', 'tensor']
......@@ -76,23 +81,68 @@ if numba:
)(abs2)
def cexp(x: ndarray) -> ndarray:
r"""Fast complex exponential.
Parameters
----------
x : ndarray
Argument of the complex exponential :math:`\exp(i x)`.
Returns
-------
y : ndarray
Complex exponential :math:`y = \exp(i x)`.
References
----------
https://software.intel.com/en-us/forums/intel-distribution-for-python/topic/758148 # noqa
"""
df_exp = np.empty(x.shape, dtype=np.complex128)
trig_buf = np.cos(x)
df_exp.real[:] = trig_buf
np.sin(x, out=trig_buf)
df_exp.imag[:] = trig_buf
return df_exp
def closest_unitary(Q: ndarray, subspace: Sequence[int] = None) -> ndarray:
"""
Compute the closest unitary to Q[np.ix_(*subspace)] on a given subspace
using left polar decomposition.
"""
if subspace is None:
subspace = (range(Q.shape[-2]), range(Q.shape[-1]))
idx = np.ix_(*subspace)
V = Q.copy()
if Q.ndim == 2:
V[idx], _ = polar(Q[idx], side='left')
else:
for i, q in enumerate(Q):
V[i][idx], _ = polar(q[idx], side='left')
return V
def tensor(*args, rank: int = 2,
optimize: Union[bool, str] = False) -> ndarray:
"""
Fast, flexible tensor product using einsum. The product is taken over the
last *rank* axes (the exception being ``rank == 1`` since a column vector
is also represented as a 2d-array) and broadcast over the remaining axes
which thus need to follow numpy broadcasting rules.
last *rank* axes and broadcast over the remaining axes which thus need to
follow numpy broadcasting rules. Note that vectors are treated as rank 2
tensors with shape (1, x) or (x, 1).
For example, the following shapes are compatible:
For example, the following shapes are compatible.
- for matrices (``rank == 2``)::
- ``rank == 2`` (e.g. matrices or vectors)::
(a, b, c, d, d), (a, b, c, e, e) -> (a, b, c, d*e, d*e)
(a, b, c), (a, d, e) -> (a, b*d, c*e)
(a, b), (c, d, e) -> (c, a*d, b*e)
- for vectors (``rank == 1``)::
(1, a), (b, 1, c) -> (b, 1, a*c)
- ``rank == 1``::
(a, b, 1, d), (a, 1, c, e) -> (a, b, c, d*e)
(a, b), (a, c) -> (a, b*c)
(a, b, 1), (a, c) -> (a, b, c)
......@@ -102,8 +152,7 @@ def tensor(*args, rank: int = 2,
The elements of the tensor product
rank : int, optional (default: 2)
The rank of the tensors. E.g., for a Kronecker product between two
vectors, ``rank == 1``, and between two matrices ``rank == 2``. The
remaining axes are broadcast over.
matrices ``rank == 2``. The remaining axes are broadcast over.
optimize : bool|str, optional (default: False)
Optimize the tensor contraction order. Passed through to
:meth:`numpy.einsum`.
......@@ -145,13 +194,6 @@ def tensor(*args, rank: int = 2,
--------
:meth:`numpy.kron`
"""
# Squeeze out extra dimension of vectors
args = [arg[None, :] if arg.ndim == 1 else arg for arg in args]
if rank == 1 and all(1 in arg.shape[-2:] for arg in args):
# Vectors, but numpy arrays are still two-dimensional
rank += 1
chars = string.ascii_letters
# All the subscripts we need
A_chars = chars[:rank]
......@@ -160,8 +202,8 @@ def tensor(*args, rank: int = 2,
A_chars, B_chars, ''.join(i + j for i, j in zip(A_chars, B_chars))
)
def tensor_product_shape(shape_A: Sequence[int], shape_B: Sequence[int],
rank: int):
def _tensor_product_shape(shape_A: Sequence[int], shape_B: Sequence[int],
rank: int):
"""Get shape of the tensor product between A and B of rank rank"""
broadcast_shape = ()
# Loop over dimensions from last to first, filling the 'shorter' shape
......@@ -198,7 +240,7 @@ def tensor(*args, rank: int = 2,
while B.ndim < rank:
B = B[None, :]
outshape = tensor_product_shape(A.shape, B.shape, rank)
outshape = _tensor_product_shape(A.shape, B.shape, rank)
return np.einsum(subscripts, A, B, optimize=optimize).reshape(outshape)
# Compute the tensor products in a binary tree-like structure, calculating
......@@ -208,8 +250,8 @@ def tensor(*args, rank: int = 2,
n = len(args)
bit = n % 2
while n > 1:
args = args[:bit] + [binary_tensor(*args[i:i+2])
for i in range(bit, n, 2)]
args = args[:bit] + tuple(binary_tensor(*args[i:i+2])
for i in range(bit, n, 2))
n = len(args)
bit = n % 2
......@@ -285,7 +327,7 @@ def pauli_expm(a: Sequence, nonzero_ind: List[bool]) -> ndarray:
Examples
--------
A Hadamard gate, i.e. a rotation by pi about X + Z (dividing by sqrt(2)
normalizes the rotation axis vector)
normalizes the rotation axis vector):
>>> a = np.full((2,), np.pi/2/np.sqrt(2))
>>> pauli_expm(a, [True, False, True]) # a[0]=a[2]=pi/2/sqrt(2), a[1]=0
......@@ -304,7 +346,8 @@ def pauli_expm(a: Sequence, nonzero_ind: List[bool]) -> ndarray:
# Replace nans originating from divide by zero
n[np.isnan(n)] = 0
real = np.einsum('ij,...->ij...', paulis[0], np.cos(theta))
imag = np.einsum('i...,ijk,...->jk...', n, paulis[1:][idx], np.sin(theta))
imag = np.einsum('i...,ijk,...->jk...', n, np.asarray(paulis[1:])[idx],
np.sin(theta))
result = real - 1j*imag
return result.transpose([2, 0, 1]) if result.ndim == 3 else result
......
......@@ -12,19 +12,50 @@ __all_ = ['load_special_measure_scan', 'cached_load_mat_file', 'special_measure_
def special_measure_to_dataframe(loaded_scan_data: dict,
squeeze_constant_setchan: bool = True) -> pandas.DataFrame:
loops = loaded_scan_data['scan']['loops'].squeeze(axis=(0, 1))
n_loops = loops.shape[-1]
n_points = loops['npoints'].squeeze(axis=2).astype(numpy.int64)
rngs = loops['rng']
setchan = loops['setchan'].squeeze(axis=2)
"""Try to interpret the data returned from hdf5storage.loadmat(filename)
as a pandas.DataFrame.
Not handled/tested yet:
- Buffered measurements
- trafofn
- procfn
"""
scan = loaded_scan_data['scan']
assert scan.shape == (1, 1)
scan = scan[0, 0]
loops = scan['loops']
assert len(loops.shape) == 2
assert loops.shape[0] == 1
loops = loops[0, :]
n_loops = loops.size
# fails if a loop has more than one npoints
npoints = list(loops['npoints'].astype(numpy.int64))
setchan = list(loops['setchan'])
for loop_idx, chan in enumerate(setchan):
assert len(chan.shape) == 2
assert chan.shape[0] == 1
setchan[loop_idx] = tuple(ch[0] for ch in chan.ravel())
# rngs can be per channel or
rngs = list(loops['rng'])
for loop_idx, rng in enumerate(rngs):
assert rng.shape == (1, 2)
rngs[loop_idx] = tuple(rng.ravel())
# This code needs to be adapted if it is possible to use different ranges
# for different setchannels in the same loop. This might require
# interpreting the trafofn
sweeps = []
for loop_idx in range(n_loops):
loop_sweeps = {}
for ch, rng, n in zip(setchan[loop_idx], rngs[loop_idx], n_points[loop_idx]):
ch = ch[0, 0]
span = numpy.linspace(*rng, num=n)
span = numpy.linspace(*rngs[loop_idx], num=npoints[loop_idx])
for ch in setchan[loop_idx]:
loop_sweeps[ch] = span
sweeps.append(loop_sweeps)
......@@ -37,12 +68,20 @@ def special_measure_to_dataframe(loaded_scan_data: dict,
idx = pandas.MultiIndex.from_product(values, names=labels)
# buffered measurements no handled?
getchan = loops['getchan']
for loop_idx, chan in enumerate(getchan):
assert len(chan.shape) == 2
assert chan.shape[0] == 1
getchan[loop_idx] = tuple(ch[0] for ch in chan.ravel())
measured = [ch[0, 0][0, 0] for ch in getchan if all(ch.shape)]
measured = list(itertools.chain.from_iterable(getchan))
# always vector cell
data = loaded_scan_data['data'].flatten()
data = loaded_scan_data['data']
assert data.shape == (1, len(measured))
data = data[0, :]
result = pandas.DataFrame(index=idx)
assert len(measured) == len(data)
......@@ -56,6 +95,11 @@ def special_measure_to_dataframe(loaded_scan_data: dict,
for lvl_idx, lvl_dim in enumerate(idx.levshape)
if lvl_dim == 1]
result = result.droplevel(to_squeeze)
try:
result.attrs['consts'] = loaded_scan_data['scan']['consts']
except ValueError:
pass
return result
......
......@@ -123,7 +123,11 @@ def plot_2d_dataframe(df: pd.DataFrame,
for mesh in ax.meshes.values():
mesh.set_clim(vmin, vmax)
ax.custom_colorbar.set_clim(vmin, vmax)
try:
ax.custom_colorbar.set_clim(vmin, vmax)
except AttributeError:
ax.custom_colorbar.mappable.set_clim(vmin, vmax)
else:
# TODO: fix
warnings.warn("for update_mode='overwrite' the colorbar code is stupid")
......
"""This module contains utilities pertaining to quantum information"""
from collections import namedtuple
from typing import Callable
import numpy as np
from . import linalg
__all__ = ['paulis', 'hadamard', 'cnot', 'cphase', 'crot', 'swap']
......@@ -20,14 +24,14 @@ paulis = namedtuple('PauliMatrices', ['sigma_0',
'sigma_x',
'sigma_y',
'sigma_z'])(
np.array([[1, 0],
[0, 1]]),
np.array([[0, 1],
[1, 0]]),
np.array([[0, -1j],
[1j, 0]]),
np.array([[1, 0],
[0, -1]]),
np.array([[1, 0],
[0, 1]]),
np.array([[0, 1],
[1, 0]]),
np.array([[0, -1j],
[1j, 0]]),
np.array([[1, 0],
[0, -1]]),
)
type(paulis).__doc__ = "Collection of pauli matrices that can be addressed by 'name' or by index"
type(paulis).__repr__ = lambda self: ('PauliMatrices(sigma_0=array([[1, 0],\n'
......@@ -64,3 +68,217 @@ swap = np.array(
[0, 1, 0, 0],
[0, 0, 0, 1]]
)
def pauli_basis(n: int):
"""n-qubit normalized Pauli basis.
Parameters
----------
n: int
Number of qubits.
Returns
-------
sigma: ndarray
n-qubit Pauli basis.
"""
normalization = np.sqrt(2**n)
combinations = np.indices((4,)*n).reshape(n, 4**n)
sigma = linalg.tensor(*np.array(paulis)[combinations], rank=2)
sigma /= normalization
return sigma
def apply_channel_kraus(rho, K):
"""Apply the Kraus representation of a channel to the operator rho.
Parameters
----------
rho: ndarray, shape (d, d)
The density operator going through the channel.
K: ndarray, shape (..., d, d)
The Kraus operators of the channel.
Returns
-------
rhop: ndarray, shape (d, d)
The output density operator of the channel.
"""
rhop = (K @ rho @ K.conj().transpose(0, 2, 1)).sum(axis=0)
return rhop
def apply_channel_chi(rho, chi):
"""Apply the Chi matrix representation of a channel to the operator rho.
Parameters
----------
rho: ndarray, shape (d, d)
The density operator going through the channel.
chi: ndarray, shape (d**2, d**2)
The chi matrix of the channel.
Returns
-------
rhop: ndarray, shape (d, d)
The output density operator of the channel.
"""
n = int(np.log2(rho.shape[0]))
sigma = pauli_basis(n)
rhop = np.einsum('jk,jab,bc,kcd', chi, sigma, rho, sigma)
return rhop
def apply_channel_liouville(rho, U, basis):
"""Apply the Liouville representation of a channel U to the operator rho.
Parameters
----------
rho: ndarray, shape (d, d)
The density operator going through the channel.
U: ndarray, shape (d**2, d**2)
The Liouville representation of the channel.
basis: ndarray, shape (d**2, d, d)
The operator basis defining the Liouville representation.
Returns
-------
rhop: ndarray, shape (d, d)
The output density operator of the channel.
"""
rhoket = np.tensordot(rho, basis, axes=[(-2, -1), (-1, -2)])
rhoketp = U @ rhoket
rhop = np.tensordot(rhoketp, basis, axes=[0, 0])
return rhop
def kraus_to_chi(K):
"""Convert from Kraus to chi matrix representation.
Parameters
----------
K: ndarray, shape (..., d, d)
The Kraus operators.
Returns
-------
chi: ndarray, shape (d**2, d**2)
The chi (process) matrix.
"""
n = int(np.log2(K.shape[-1]))
sigma = pauli_basis(n)
a_ij = np.einsum('iab,jba', K, sigma)
chi = (a_ij.conj().T @ a_ij).T
return chi
def chi_to_liouville(chi, basis):
"""Convert from chi to Liouville representation.
Parameters
----------
chi: ndarray, shape (d**2, d**2)
The chi (process) matrix.
basis: ndarray, shape (d**2, d, d)
The operator basis defining the Liouville representation.
Returns
-------
U: ndarray, shape (d**2, d**2)
The Liouville representation.
"""
traces = np.einsum('iab,jbc,kcd,lda', *[basis]*4)
U = np.einsum('kl,ikjl', chi, traces)
return U
def kraus_to_liouville(K, basis):
"""Convert from Kraus to Liouville representation.
Parameters
----------
K: ndarray, shape (..., d, d)
The Kraus operators.
basis: ndarray, shape (d**2, d, d)
The operator basis defining the Liouville representation.
Returns
-------
U: ndarray, shape (d**2, d**2)
The Liouville representation.
"""
U = np.einsum('iab,kbc,jcd,kad', basis, K, basis, K.conj())
return U
def liouville_to_choi(U, basis):
"""Convert from Liouville to Choi matrix representation.
Parameters
----------
U: ndarray, shape (d**2, d**2)
The Liouville representation.
basis: ndarray, shape (d**2, d, d)
The operator basis defining the Liouville representation.
Returns
-------
choi: ndarray, shape (d**2, d**2)
The Choi matrix of the channel U.
"""
d2 = U.shape[0]
try:
basis_kroned = np.einsum(
'iab,jcd->ijacbd', basis.transpose(0, 2, 1), basis
).reshape(d2, d2, d2, d2)
choi = np.tensordot(U, basis_kroned, axes=[(0, 1), (1, 0)])
except MemoryError:
basis_kroned = np.empty((d2, d2, d2), dtype=complex)
choi = np.empty((d2, d2), dtype=complex)
for i in range(d2):
basis_kroned = np.einsum(
'ab,jcd->ijacbd', basis[i].transpose(0, 2, 1), basis,
out=basis_kroned
).reshape(d2, d2, d2)
choi[i] = np.tensordot(U, basis_kroned, axes=[(0, 1), (1, 0)])
choi /= np.sqrt(d2)
return choi
def channel_to_choi(U: Callable, d: int, *args, **kwargs):
"""Convert a channel U to the Choi representation.
Parameters
----------
U: Callable
Callable that takes an operator as input and returns another operator.
d: int
Dimension of the Hilbert space.
*args, **kwargs: Additional arguments passed through to U.
Returns
-------
choi: ndarray, shape (d**2, d**2)
The Choi matrix of the channel U.
"""
eye = np.eye(d)
choi = np.zeros((d**2, d**2), dtype=complex)
for i in range(d):
for j in range(d):
ij = linalg.tensor(eye[:, i:i+1], eye[j:j+1])
choi += linalg.tensor(ij, U(ij, *args, **kwargs))
choi /= d
return choi
"""
This module contains random functions
"""
import numpy as np
from numpy import ndarray
def random_hermitian</