Commit e945e70b authored by Tobias Hangleiter's avatar Tobias Hangleiter
Browse files

Merge branch 'tobias'

parents 6275fc89 2f7743e7
......@@ -18,11 +18,21 @@ This will link the files into your environment instead of copying them. If you a
## qutil.matlab
In this module there are functions that are helpful for reading `.mat` files, especially those created with special measure. If you simply want to open a random `.mat` file you can use `hdf5storage.loadmat`.
## qutil.const
This module defines all the constants you could wish for as well as functions to convert temperatures (`convert_temperature`) or between wavelengths and frequencies (`lambda2nu`, `nu2lambda`). For an overview, see the module docstring.
## qutil.linalg
This module provides several handy linear algebra functions. While some are implemented elsewhere, the implementation here is typically speedier for large arrays. For example, `pauli_expm` exploits the fact that a matrix exponential of Pauli matrices can be written as a cosine times the identity matrix plus a sine times the Paulis to speed up the calculation.
For an overview of the included functions, see the module docstring.
## qutil.ui
This module collects UI helpers, such as a progress bar for loops that can be used like so:
```python
for i in qutil.ui.progressbar(range(n)):
do_something()
```
## qutil.qi
In this module there are some quantities and functions related to quantum information, like the Pauli matrices in different data types.
......
from . import const, linalg, matlab, plotting, ui, qi
__version__ = '0.1'
__all__ = ['const', 'linalg', 'matlab', 'ui', 'plotting', 'qi']
r"""
This module defines constants, mostly wrapped from scipy. Many constants are in
the top-level namespace. Many more are in the dictionary ``physical_constants``
that can be searched using :meth:`find` and accessed using :meth:`value`.
On top of the scipy constants, the following shorthands are defined:
``e2`` : :math:`e^2`
``pi2`` : :math:`\pi/2`
``pi4`` : :math:`\pi/4`
``two_pi`` : :math:`2\pi`
``four_pi`` : :math:`4\pi`
Functions
---------
:meth:`convert_temperature` :
Convert temperature between different scales
:meth:`find` :
Find constants in the ``physical_constants`` dictionary matching a string
:meth:`lambda2nu` :
Convert wavelength to frequency
:meth:`nu2lambda` :
Convert frequency to wavelength
:meth:`unit` :
Get the unit of a constant in the ``physical_constants`` dictionary.
:meth:`value` :
Get the value of a constant in the ``physical_constants`` dictionary.
:meth:`precision` :
Get the precision of a constant in the ``physical_constants`` dictionary.
"""
from scipy.constants import (
Avogadro,
Boltzmann,
Btu,
Btu_IT,
Btu_th,
G,
Julian_year,
N_A,
Planck,
R,
Rydberg,
Stefan_Boltzmann,
Wien,
acre,
alpha,
angstrom,
arcmin,
arcminute,
arcsec,
arcsecond,
astronomical_unit,
atm,
atmosphere,
atomic_mass,
atto,
au,
bar,
barrel,
bbl,
blob,
c,
calorie,
calorie_IT,
calorie_th,
carat,
centi,
convert_temperature,
day,
deci,
degree,
degree_Fahrenheit,
deka,
dyn,
dyne,
e,
eV,
electron_mass,
electron_volt,
elementary_charge,
epsilon_0,
erg,
exa,
exbi,
femto,
fermi,
find,
fine_structure,
fluid_ounce,
fluid_ounce_US,
fluid_ounce_imp,
foot,
g,
gallon,
gallon_US,
gallon_imp,
gas_constant,
gibi,
giga,
golden,
golden_ratio,
grain,
gram,
gravitational_constant,
h,
hbar,
hectare,
hecto,
horsepower,
hour,
hp,
inch,
k,
kgf,
kibi,
kilo,
kilogram_force,
kmh,
knot,
lambda2nu,
lb,
lbf,
light_year,
liter,
litre,
long_ton,
m_e,
m_n,
m_p,
m_u,
mach,
mebi,
mega,
metric_ton,
micro,
micron,
mil,
mile,
milli,
minute,
mmHg,
mph,
mu_0,
nano,
nautical_mile,
neutron_mass,
nu2lambda,
ounce,
oz,
parsec,
pebi,
peta,
physical_constants,
pi,
pico,
point,
pound,
pound_force,
precision,
proton_mass,
psi,
pt,
short_ton,
sigma,
slinch,
slug,
speed_of_light,
speed_of_sound,
stone,
survey_foot,
survey_mile,
tebi,
tera,
ton_TNT,
torr,
troy_ounce,
troy_pound,
u,
unit,
value,
week,
yard,
year,
yobi,
yotta,
zebi,
zepto,
zero_Celsius,
zetta
)
e2 = e**2
pi2 = pi/2
pi4 = pi/4
two_pi = pi*2
four_pi = pi*4
......@@ -6,8 +6,8 @@ Functions
:meth:`abs2` :
Fast absolute value squared
:meth:`tensor`
Fast tensor product of an arbitrary number of inputs using
:meth:`~numpy.einsum`
Fast, flexible tensor product of an arbitrary number of inputs using
:meth:`~numpy.einsum` that supports broadcasting
:meth:`mdot`
Multiple matrix product along an axis
:meth:`ptrace`
......@@ -25,15 +25,15 @@ Functions
:meth:`dot_HS`
Hilbert-Schmidt inner product
"""
import operator
import string
from functools import reduce
from math import copysign
from itertools import zip_longest
from typing import List, Sequence, Tuple, Union
import numpy as np
import qutip as qt
from numpy import ndarray
from numpy.linalg import norm
from numpy import linalg, ndarray
from .qi import P_np
......@@ -43,6 +43,10 @@ try:
except ImportError:
numba = False
__all__ = ['abs2', 'check_phase_eq', 'density', 'dot_HS', 'max_abs_diff',
'max_rel_diff', 'mdot', 'pauli_expm', 'ptrace',
'remove_float_errors', 'sparsity', 'tensor']
def abs2(arr):
r"""
......@@ -73,29 +77,145 @@ if numba:
)(abs2)
def tensor(*args):
r"""
Fast tensor product using einsum. The arguments may be of arbitrary
dimension but are required to be square on its last two axes and the shapes
of the remaining axes must be the same. For example, the following shapes
are compatible::
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.
For example, the following shapes are compatible.
- for matrices (``rank == 2``)::
(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``)::
"""
(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)
if len(set(arg.shape[:-2] for arg in args)) != 1:
raise ValueError('Require all args to have the same shape except ' +
'on the last two axes.')
if not all(arg.shape[-2] == arg.shape[-1] for arg in args):
raise ValueError('Require all args to be square in its last two axes.')
Parameters
----------
args : array_like
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.
optimize : bool|str, optional (default: False)
Optimize the tensor contraction order. Passed through to
:meth:`numpy.einsum`.
Examples
--------
>>> Z = np.diag([1, -1])
>>> np.array_equal(tensor(Z, Z), np.kron(Z, Z))
True
>>> A, B = np.arange(2), np.arange(2, 5)
>>> tensor(A, B, rank=1)
array([[0, 0, 0, 2, 3, 4]])
>>> args = np.random.randn(4, 10, 3, 2)
>>> result = tensor(*args, rank=1)
>>> result.shape == (10, 3, 2**4)
True
>>> result = tensor(*args, rank=2)
>>> result.shape == (10, 3**4, 2**4)
True
>>> A, B = np.random.randn(1, 3), np.random.randn(3, 4)
>>> result = tensor(A, B)
>>> result.shape == (1*3, 3*4)
True
>>> A, B = np.random.randn(3, 1, 2), np.random.randn(2, 2, 2)
>>> try:
... result = tensor(A, B, rank=2)
... except ValueError as err: # cannot broadcast over axis 0
... print(err)
Incompatible shapes (3, 1, 2) and (2, 2, 2) for tensor product of rank 2.
>>> result = tensor(A, B, rank=3)
>>> result.shape == (3*2, 1*2, 2*2)
True
See Also
--------
: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]
B_chars = chars[rank:2*rank]
subscripts = '...{},...{}->...{}'.format(
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):
"""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
# with 1's once it is exhausted
for dims in zip_longest(shape_A[-rank-1::-1], shape_B[-rank-1::-1],
fillvalue=1):
if 1 in dims:
# Broadcast 1-d of argument to dimension of other
broadcast_shape = (max(dims),) + broadcast_shape
elif len(set(dims)) == 1:
# Both arguments have same dimension on axis.
broadcast_shape = dims[:1] + broadcast_shape
else:
raise ValueError('Incompatible shapes ' +
'{} and {} '.format(shape_A, shape_B) +
'for tensor product of rank {}.'.format(rank))
# Shape of the actual tensor product is product of each dimension,
# again broadcasting if need be
tensor_shape = tuple(
reduce(operator.mul, dimensions) for dimensions in zip_longest(
shape_A[:-rank-1:-1], shape_B[:-rank-1:-1], fillvalue=1
)
)[::-1]
return broadcast_shape + tensor_shape
def binary_tensor(A, B):
d = A.shape[-1]*B.shape[-1]
return np.einsum('...ij,...kl->...ikjl', A, B).reshape(*A.shape[:-2],
d, d)
"""Compute the Kronecker product of two tensors"""
# Add dimensions so that each arg has at least ndim == rank
while A.ndim < rank:
A = A[None, :]
while B.ndim < rank:
B = B[None, :]
outshape = tensor_product_shape(A.shape, B.shape, rank)
return np.einsum(subscripts, A, B, optimize=optimize).reshape(outshape)
return reduce(binary_tensor, args)
# Compute the tensor products in a binary tree-like structure, calculating
# the product of two leaves and working up. This is more memory-efficient
# than reduce(binary_tensor, args) which computes the products
# left-to-right.
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)]
n = len(args)
bit = n % 2
return args[0]
def mdot(arr: Sequence, axis: int = 0) -> ndarray:
......@@ -179,7 +299,7 @@ def pauli_expm(a: Sequence, nonzero_ind: List[bool]) -> ndarray:
P = np.asarray(P_np)
idx = np.asarray(nonzero_ind)
# Twice the rotation angle
theta = norm(a, axis=0)
theta = linalg.norm(a, axis=0)
# Rotation axis
n = a/theta
# Replace nans originating from divide by zero
......@@ -213,8 +333,8 @@ def remove_float_errors(arr: ndarray, eps_scale: float = None):
return arr
def check_phase_eq(psi: Union[qt.Qobj, ndarray],
phi: Union[qt.Qobj, ndarray],
def check_phase_eq(psi: Union[qt.Qobj, Sequence],
phi: Union[qt.Qobj, Sequence],
eps: float = None,
normalized: bool = False) -> Tuple[bool, float]:
r"""
......@@ -229,7 +349,7 @@ def check_phase_eq(psi: Union[qt.Qobj, ndarray],
Parameters
----------
psi, phi : Qobj or ndarray
psi, phi : Qobj or array_like
Vectors or operators to be compared
eps : float
The tolerance below which the two objects are treated as equal, i.e.,
......@@ -245,36 +365,39 @@ def check_phase_eq(psi: Union[qt.Qobj, ndarray],
>>> check_phase_eq(psi, phi)
(True, 1.2345)
"""
d = max(psi.shape)
psi, phi = [obj.full() if isinstance(obj, qt.Qobj) else obj
for obj in (psi, phi)]
if eps is None:
# Tolerance the floating point eps times the dimension squared
eps = max(np.finfo(psi.dtype).eps, np.finfo(phi.dtype).eps)*d**2
# Tolerance the floating point eps times the # of flops for the matrix
# multiplication, i.e. for psi and phi n x m matrices 2*n**2*m
eps = max(np.finfo(psi.dtype).eps, np.finfo(phi.dtype).eps) *\
np.prod(psi.shape)*phi.shape[-1]*2
if not normalized:
# normalization introduces more floating point error
eps *= (np.prod(psi.shape)*phi.shape[-1]*2)**2
if psi.ndim - psi.shape.count(1) == 1:
# Vector
inner_product = (psi.T.conj() @ phi).squeeze()
if normalized:
inner_product = (psi.T.conj() @ phi).squeeze()
norm = 1
else:
inner_product = ((psi.T.conj() @ phi) /
(norm(psi)*norm(phi))).squeeze()
norm = (linalg.norm(psi)*linalg.norm(phi)).squeeze()
elif psi.ndim == 2:
inner_product = dot_HS(psi, phi, eps)
# Matrix
if normalized:
inner_product = dot_HS(psi, phi, eps)
norm = 1
else:
inner_product = dot_HS(psi, phi, eps) /\
np.sqrt(dot_HS(psi, psi, eps)*dot_HS(phi, phi, eps))
norm = np.sqrt(dot_HS(psi, psi, eps)*dot_HS(phi, phi, eps))
else:
raise ValueError('Invalid dimension')
phase = np.angle(inner_product)
modulus = abs(inner_product)
return abs(1 - modulus) <= eps, phase
return abs(norm - modulus) <= eps, phase
def dot_HS(U: Union[ndarray, qt.Qobj],
......@@ -290,6 +413,11 @@ def dot_HS(U: Union[ndarray, qt.Qobj],
U, V : Qobj or ndarray
Objects to compute the inner product of.
Returns
-------
result : float, complex
The result rounded to precision eps.
Examples
--------
>>> U, V = qt.sigmax(), qt.sigmay()
......@@ -303,33 +431,22 @@ def dot_HS(U: Union[ndarray, qt.Qobj],
if isinstance(V, qt.Qobj):
V = V.full()
res = np.einsum('ij,ij', U.conj(), V)
if eps is None:
# Tolerance is the dtype precision times the number of flops for the
# matrix multiplication
# matrix multiplication times two to be on the safe side
try:
eps = np.finfo(U.dtype).eps*np.prod(U.shape)*V.shape[-1]
eps = np.finfo(U.dtype).eps*np.prod(U.shape)*V.shape[-1]*2
except ValueError:
# dtype is int and therefore exact
eps = 0
# Deal with real and imaginary part separately
if abs(res.real) <= eps:
real = 0.0
elif abs(1 - abs(res.real)) <= eps:
real = copysign(1, res.real)
else:
real = res.real
if abs(res.imag) <= eps:
imag = 0
elif abs(1 - abs(res.imag)) <= eps:
imag = copysign(1, res.imag)
if eps == 0:
decimals = 0
else:
imag = res.imag
decimals = abs(int(np.log10(eps)))
return real + 1j*imag if imag != 0 else real
res = np.round(np.einsum('ij,ij', U.conj(), V), decimals)
return res if res.imag else res.real
def sparsity(arr: ndarray, eps: float = None) -> float:
......
......@@ -16,3 +16,5 @@ P_sp = [P.data for P in P_qt]
# Clifford group generators
clifford_generators = [qt.snot(), qt.phasegate(np.pi/2).tidyup(), qt.cnot()]
__all__ = ['clifford_generators', 'P_np', 'P_sp', 'P_qt']
"""
This module contains UI utility functions
Functions
---------
:meth:`progressbar` :
A progress bar for loops
"""
import io
import json
import os
import re
import sys
from typing import Generator, Iterable
try:
import jupyter_client
import requests
from jupyter_core.paths import jupyter_runtime_dir
from notebook.utils import check_pid
from requests.compat import urljoin
def _list_running_servers(runtime_dir: str = None) -> Generator:
"""Iterate over the server info files of running notebook servers.
Given a runtime directory, find nbserver-* files in the security
directory, and yield dicts of their information, each one pertaining to
a currently running notebook server instance.
Copied from notebook.notebookapp.list_running_servers() (version 5.7.8)
since the highest version compatible with Python 3.5 (version 5.6.0)
has a bug.
"""
if runtime_dir is None:
runtime_dir = jupyter_runtime_dir()
# The runtime dir might not exist
if not os.path.isdir(runtime_dir):
return
for file_name in os.listdir(runtime_dir):
if re.match('nbserver-(.+).json', file_name):
with io.open(os.path.join(runtime_dir, file_name),
encoding='utf-8') as f:
info = json.load(f)
# Simple check whether that process is really still running
# Also remove leftover files from IPython 2.x without a pid
# field
if ('pid' in info) and check_pid(info['pid']):
yield info
else:
# If the process has died, try to delete its info file
try:
os.unlink(os.path.join(runtime_dir, file_name))
except OSError:
pass # TODO: This should warn or log or something
def _get_notebook_name() -> str: