# Source code for qutip.tensor

# This file is part of QuTiP: Quantum Toolbox in Python.
#
#    Copyright (c) 2011 and later, Paul D. Nation and Robert J. Johansson.
#
#    Redistribution and use in source and binary forms, with or without
#    modification, are permitted provided that the following conditions are
#    met:
#
#    1. Redistributions of source code must retain the above copyright notice,
#       this list of conditions and the following disclaimer.
#
#    2. Redistributions in binary form must reproduce the above copyright
#       notice, this list of conditions and the following disclaimer in the
#       documentation and/or other materials provided with the distribution.
#
#    3. Neither the name of the QuTiP: Quantum Toolbox in Python nor the names
#       of its contributors may be used to endorse or promote products derived
#       from this software without specific prior written permission.
#
#    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
#    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
#    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
#    PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
#    HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
#    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
#    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
#    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
#    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
#    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
#    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
###############################################################################
"""
Module for the creation of composite quantum objects via the tensor product.
"""

__all__ = ['tensor', 'super_tensor', 'composite', 'tensor_contract']

import numpy as np
import scipy.sparse as sp

from qutip.qobj import Qobj
from qutip.permute import reshuffle
from qutip.superoperator import operator_to_vector

import qutip.settings
import qutip.superop_reps  # Avoid circular dependency here.

[docs]def tensor(*args): """Calculates the tensor product of input operators. Parameters ---------- args : array_like list or array of quantum objects for tensor product. Returns ------- obj : qobj A composite quantum object. Examples -------- >>> tensor([sigmax(), sigmax()]) Quantum object: dims = [[2, 2], [2, 2]], \ shape = [4, 4], type = oper, isHerm = True Qobj data = [[ 0.+0.j 0.+0.j 0.+0.j 1.+0.j] [ 0.+0.j 0.+0.j 1.+0.j 0.+0.j] [ 0.+0.j 1.+0.j 0.+0.j 0.+0.j] [ 1.+0.j 0.+0.j 0.+0.j 0.+0.j]] """ if not args: raise TypeError("Requires at least one input argument") if len(args) == 1 and isinstance(args[0], (list, np.ndarray)): # this is the case when tensor is called on the form: # tensor([q1, q2, q3, ...]) qlist = args[0] elif len(args) == 1 and isinstance(args[0], Qobj): # tensor is called with a single Qobj as an argument, do nothing return args[0] else: # this is the case when tensor is called on the form: # tensor(q1, q2, q3, ...) qlist = args if not all([isinstance(q, Qobj) for q in qlist]): # raise error if one of the inputs is not a quantum object raise TypeError("One of inputs is not a quantum object") out = Qobj() if qlist[0].issuper: out.superrep = qlist[0].superrep if not all([q.superrep == out.superrep for q in qlist]): raise TypeError("In tensor products of superroperators, all must" + "have the same representation") out.isherm = True for n, q in enumerate(qlist): if n == 0: out.data = q.data out.dims = q.dims else: out.data = sp.kron(out.data, q.data, format='csr') out.dims = [out.dims[0] + q.dims[0], out.dims[1] + q.dims[1]] out.isherm = out.isherm and q.isherm if not out.isherm: out._isherm = None return out.tidyup() if qutip.settings.auto_tidyup else out
[docs]def super_tensor(*args): """Calculates the tensor product of input superoperators, by tensoring together the underlying Hilbert spaces on which each vectorized operator acts. Parameters ---------- args : array_like list or array of quantum objects with type="super". Returns ------- obj : qobj A composite quantum object. """ if isinstance(args[0], list): args = args[0] # Check if we're tensoring vectors or superoperators. if all(arg.issuper for arg in args): if not all(arg.superrep == "super" for arg in args): raise TypeError( "super_tensor on type='super' is only implemented for " "superrep='super'." ) # Reshuffle the superoperators. shuffled_ops = list(map(reshuffle, args)) # Tensor the result. shuffled_tensor = tensor(shuffled_ops) # Unshuffle and return. out = reshuffle(shuffled_tensor) out.superrep = args[0].superrep return out elif all(arg.isoperket for arg in args): # Reshuffle the superoperators. shuffled_ops = list(map(reshuffle, args)) # Tensor the result. shuffled_tensor = tensor(shuffled_ops) # Unshuffle and return. out = reshuffle(shuffled_tensor) return out elif all(arg.isoperbra for arg in args): return super_tensor(*(arg.dag() for arg in args)).dag() else: raise TypeError( "All arguments must be the same type, " "either super, operator-ket or operator-bra." )
def _isoperlike(q): return q.isoper or q.issuper def _isketlike(q): return q.isket or q.isoperket def _isbralike(q): return q.isbra or q.isoperbra
[docs]def composite(*args): """ Given two or more operators, kets or bras, returns the Qobj corresponding to a composite system over each argument. For ordinary operators and vectors, this is the tensor product, while for superoperators and vectorized operators, this is the column-reshuffled tensor product. If a mix of Qobjs supported on Hilbert and Liouville spaces are passed in, the former are promoted. Ordinary operators are assumed to be unitaries, and are promoted using to_super, while kets and bras are promoted by taking their projectors and using operator_to_vector(ket2dm(arg)). """ # First step will be to ensure everything is a Qobj at all. if not all(isinstance(arg, Qobj) for arg in args): raise TypeError("All arguments must be Qobjs.") # Next, figure out if we have something oper-like (isoper or issuper), # or something ket-like (isket or isoperket). Bra-like we'll deal with # by turning things into ket-likes and back. if all(map(_isoperlike, args)): # OK, we have oper/supers. if any(arg.issuper for arg in args): # Note that to_super does nothing to things # that are already type=super, while it will # promote unitaries to superunitaries. return super_tensor(*map(qutip.superop_reps.to_super, args)) else: # Everything's just an oper, so ordinary tensor products work. return tensor(*args) elif all(map(_isketlike, args)): # Ket-likes. if any(arg.isoperket for arg in args): # We have a vectorized operator, we we may need to promote # something. return super_tensor(*( arg if arg.isoperket else operator_to_vector(qutip.states.ket2dm(arg)) for arg in args )) else: # Everything's ordinary, so we can use the tensor product here. return tensor(*args) elif all(map(_isbralike, args)): # Turn into ket-likes and recurse. return composite(*(arg.dag() for arg in args)).dag() else: raise TypeError("Unsupported Qobj types [{}].".format( ", ".join(arg.type for arg in args) ))
def flatten(l): """Flattens a list of lists to the first level. Given a list containing a mix of scalars and lists, flattens down to a list of the scalars within the original list. Examples -------- >>> print(flatten([[[0], 1], 2])) [0, 1, 2] """ if not isinstance(l, list): return [l] else: return sum(map(flatten, l), []) def _enumerate_flat(l, idx=0): if not isinstance(l, list): # Found a scalar, so return and increment. return idx, idx + 1 else: # Found a list, so append all the scalars # from it and recurse to keep the increment # correct. acc = [] for elem in l: labels, idx = _enumerate_flat(elem, idx) acc.append(labels) return acc, idx def enumerate_flat(l): """Labels the indices at which scalars occur in a flattened list. Given a list containing a mix of scalars and lists, returns a list of the same structure, where each scalar has been replaced by an index into the flattened list. Examples -------- >>> print(enumerate_flat([[[10], [20, 30]], 40])) [[[0], [1, 2]], 3] """ return _enumerate_flat(l)[0] def deep_remove(l, *what): """Removes scalars from all levels of a nested list. Given a list containing a mix of scalars and lists, returns a list of the same structure, but where one or more scalars have been removed. Examples -------- >>> print(deep_remove([[[[0, 1, 2]], [3, 4], [5], [6, 7]]], 0, 5)) [[[[1, 2]], [3, 4], [], [6, 7]]] """ if isinstance(l, list): # Make a shallow copy at this level. l = l[:] for to_remove in what: if to_remove in l: l.remove(to_remove) else: l = list(map(lambda elem: deep_remove(elem, to_remove), l)) return l def unflatten(l, idxs): """Unflattens a list by a given structure. Given a list of scalars and a deep list of indices as produced by flatten, returns an "unflattened" form of the list. This perfectly inverts flatten. Examples -------- >>> l = [[[10, 20, 30], [40, 50, 60]], [[70, 80, 90], [100, 110, 120]]] >>> idxs = enumerate_flat(l) >>> print(unflatten(flatten(l)), idxs) == l True """ acc = [] for idx in idxs: if isinstance(idx, list): acc.append(unflatten(l, idx)) else: acc.append(l[idx]) return acc def _tensor_contract_single(arr, i, j): """ Contracts a dense tensor along a single index pair. """ if arr.shape[i] != arr.shape[j]: raise ValueError("Cannot contract over indices of different length.") idxs = np.arange(arr.shape[i]) sl = tuple(slice(None, None, None) if idx not in (i, j) else idxs for idx in range(arr.ndim)) return np.sum(arr[sl], axis=0) def _tensor_contract_dense(arr, *pairs): """ Contracts a dense tensor along one or more index pairs, keeping track of how the indices are relabeled by the removal of other indices. """ axis_idxs = list(range(arr.ndim)) for pair in pairs: # axis_idxs.index effectively evaluates the mapping from # original index labels to the labels after contraction. arr = _tensor_contract_single(arr, *map(axis_idxs.index, pair)) list(map(axis_idxs.remove, pair)) return arr
[docs]def tensor_contract(qobj, *pairs): """Contracts a qobj along one or more index pairs. Note that this uses dense representations and thus should *not* be used for very large Qobjs. Parameters ---------- pairs : tuple One or more tuples (i, j) indicating that the i and j dimensions of the original qobj should be contracted. Returns ------- cqobj : Qobj The original Qobj with all named index pairs contracted away. """ # Record and label the original dims. dims = qobj.dims dims_idxs = enumerate_flat(dims) flat_dims = flatten(dims) # Convert to dense first, since sparse won't support the reshaping we need. qtens = qobj.data.toarray() # Reshape by the flattened dims. qtens = qtens.reshape(flat_dims) # Contract out the indices from the flattened object. qtens = _tensor_contract_dense(qtens, *pairs) # Remove the contracted indexes from dims so we know how to # reshape back. contracted_idxs = deep_remove(dims_idxs, *flatten(list(map(list, pairs)))) contracted_dims = unflatten(flat_dims, contracted_idxs) l_mtx_dims, r_mtx_dims = map(np.product, contracted_dims) # Reshape back into a 2D matrix. qmtx = qtens.reshape((l_mtx_dims, r_mtx_dims)) # Return back as a qobj. return Qobj(qmtx, dims=contracted_dims, superrep=qobj.superrep)
import qutip.states