Commit 2094bdc1 authored by Tobias Hangleiter's avatar Tobias Hangleiter
Browse files

Make linalg.tensor more flexible and add UI module

The UI module for now just has a progress bar for loops.
parent c7013555
from . import linalg, matlab, plotting, ui, qi
__version__ = '0.1'
__all__ = ['linalg', 'matlab', 'ui', 'plotting', 'qi']
......@@ -25,9 +25,9 @@ Functions
:meth:`dot_HS`
Hilbert-Schmidt inner product
"""
import operator
import string
from functools import reduce
from math import copysign
from typing import List, Sequence, Tuple, Union
import numpy as 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,26 +77,95 @@ if numba:
)(abs2)
def tensor(*args):
def tensor(*args, rank: int = 2, optimize: Union[bool, str] = False):
r"""
Fast tensor product using einsum. The arguments may be of arbitrary
dimension but are required to be square or pairwise transposed on its last
two axes and the shapes of the remaining axes must be the same. For
example, the following shapes are compatible::
Fast 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 be of the same shape.
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, c, b) -> (a, b*c, c*b)
(a, b, c), (a, d, e) -> (a, b*d, c*e)
- for vectors (``rank == 1``)::
(a, b, c, d), (a, b, c, e) -> (a, b, c, d*e)
(a, b), (a, c) -> (a, b*c)
(a, b, 1), (a, c) -> (a, b, c)
Parameters
----------
args : array_like
The elements of the tensor product
rank : int, optional
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
Optimize the tensor contraction order. Passed through to
:meth:`oe.contract`.
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)
>>> try:
... result = tensor(A, B, rank=1)
... except ValueError as err:
... print(err)
Require all args to have the same shape except on the last 1 axes.
>>> result = tensor(A, B, rank=2)
>>> result.shape == (1*3, 3*4)
True
"""
# Squeeze out extra dimension of vectors
args = [arg[None, :] if arg.ndim == 1 else arg for arg in args]
if all(1 in arg.shape[-2:] for arg in args) and rank == 1:
# Vectors, but numpy arrays are still two-dimensional
rank += 1
if len(set(arg.shape[:-2] for arg in args)) != 1:
if len(set(arg.shape[:-rank] for arg in args)) != 1:
raise ValueError('Require all args to have the same shape except ' +
'on the last two axes.')
'on the last {} axes.'.format(rank))
chars = string.ascii_lowercase + string.ascii_uppercase
# 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 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)
"""Function to compute the Kronecker product of two tensors"""
if optimize:
path, _ = np.einsum_path(subscripts, A, B, optimize=optimize)
else:
path = False
outshape = A.shape[:-rank] + tuple(
reduce(operator.mul, dimensions)
for dimensions in zip(*[arg.shape[-rank:] for arg in (A, B)])
)
return np.einsum(subscripts, A, B, optimize=path).reshape(outshape)
return reduce(binary_tensor, args)
......@@ -249,31 +322,34 @@ def check_phase_eq(psi: Union[qt.Qobj, Sequence],
if eps is None:
# Tolerance the floating point eps times the # of flops for the matrix
# multiplication, i.e. for psi and phi n x m matrices n**2*m
# 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]
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],
......@@ -309,9 +385,9 @@ def dot_HS(U: Union[ndarray, qt.Qobj],
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
......
......@@ -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 Iterable, List, Union
import requests
from requests.compat import urljoin
try:
import jupyter_client
from notebook.utils import check_pid
from jupyter_core.paths import jupyter_runtime_dir
_has_notebook = True
except ImportError:
_has_notebook = False
def _list_running_servers(runtime_dir: str = None) -> List[str]:
"""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() -> Union[None, str]:
"""
Return the full path of the jupyter notebook.
See https://github.com/jupyter/notebook/issues/1000
"""
try:
connection_file = jupyter_client.find_connection_file()
except OSError:
return
kernel_id = re.search('kernel-(.*).json', connection_file).group(1)
servers = _list_running_servers()
for ss in servers:
response = requests.get(urljoin(ss['url'], 'api/sessions'),
params={'token': ss.get('token', '')})
for nn in json.loads(response.text):
if nn['kernel']['id'] == kernel_id:
relative_path = nn['notebook']['path']
return os.path.join(ss['notebook_dir'], relative_path)
try:
if _has_notebook and _get_notebook_name():
from tqdm import tqdm_notebook as tqdm
else:
from tqdm import tqdm
_has_tqdm = True
except ImportError:
_has_tqdm = False
__all__ = ['progressbar']
def _progressbar(iterable: Iterable, prefix: str = "Computing: ",
size: int = 25, file=sys.stdout):
"""
Primitive base progress bar stolen from
https://stackoverflow.com/a/34482761
"""
count = len(iterable)
def show(j):
x = int(size*j/count)
file.write("\r{}[{}{}] {} %".format(prefix, "#"*x, "."*(size - x),
int(100*j/count)))
file.flush()
show(0)
for i, item in enumerate(iterable):
yield item
show(i + 1)
file.write("\n")
file.flush()
def progressbar(iterable: Iterable, *args, **kwargs):
"""
Progress bar for loops. Uses tqdm if available or a quick-and-dirty
implementation from stackoverflow.
Usage::
for i in progressbar(range(10)):
do_something()
"""
if _has_tqdm:
return tqdm(iterable, *args, **kwargs)
else:
return _progressbar(iterable, *args, **kwargs)
......@@ -34,6 +34,7 @@ setup(name='qutil',
package_dir={'qutil': 'qutil'},
python_requires='>=3.5',
install_requires=['numpy', 'hdf5storage', 'pandas'],
extras_require={'fancy_progressbar': ['tqdm']},
classifiers=[
"Programming Language :: Python :: 3",
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment