def test_raise_error(self):
amat = matrix()
bmat = matrix()
bvec = vector()
# Test invalid length for axes
self.assertRaises(ValueError, tensordot, amat, bmat, (0, 1, 2))
# Test axes of uneven length
self.assertRaises(ValueError, tensordot, amat, bmat, ((0, 1), (0)))
# Test invalid len(axes) given inputs are matrices
self.assertRaises(ValueError, tensordot, amat, bmat, ((0, 1, 2), (0, 1, 2)))
# Test invalid axes[1] given that y is a vector
self.assertRaises(ValueError, tensordot, amat, bvec, (0, 1))
# Test invalid scalar axes given inputs are matrices
self.assertRaises(ValueError, tensordot, amat, bvec, 2)
python类tensordot()的实例源码
def test_weird_valid_axes(self):
# Test matrix-matrix
amat = matrix()
bmat = matrix()
for axes in [0,
(1, 0),
[1, 0],
(1, (0, )),
((1, ), 0),
([1], [0]),
([], [])]:
c = tensordot(amat, bmat, axes)
f3 = inplace_func([amat, bmat], c)
aval = rand(4, 7)
bval = rand(7, 9)
self.assertTrue(numpy.allclose(numpy.tensordot(aval, bval, axes),
f3(aval, bval)))
utt.verify_grad(self.TensorDot(axes), [aval, bval])
def test_scalar_axes(self):
# Test matrix-matrix
amat = fmatrix()
bmat = dmatrix()
# We let at float64 to test mix of float32 and float64.
axes = 1
aval = rand(4, 5).astype('float32')
bval = rand(5, 3)
c = tensordot(amat, bmat, axes)
f3 = inplace_func([amat, bmat], c)
self.assertTrue(numpy.allclose(numpy.tensordot(aval, bval, axes),
f3(aval, bval)))
utt.verify_grad(self.TensorDot(axes), [aval, bval])
# Test tensor-tensor
amat = tensor3()
bmat = tensor3()
axes = 2
aval = rand(3, 4, 5)
bval = rand(4, 5, 3)
c = tensordot(amat, bmat, axes)
f3 = inplace_func([amat, bmat], c)
self.assertTrue(numpy.allclose(numpy.tensordot(aval, bval, axes),
f3(aval, bval)))
utt.verify_grad(self.TensorDot(axes), [aval, bval])
def _develop(self, basis):
"""Compute the AR models and gains at instants fixed by newcols
returns:
ARcols : array containing the AR parts
Gcols : array containing the gains
The size of ARcols is (1 + ordar_, n_epochs, n_points)
The size of Gcols is (1, n_epochs, n_points)
"""
n_basis, n_epochs, n_points = basis.shape
ordar_ = self.ordar_
# -------- expand on the basis
AR_cols_ones = np.ones((1, n_epochs, n_points))
if ordar_ > 0:
AR_cols = np.tensordot(self.AR_, basis, axes=([1], [0]))
AR_cols = np.vstack((AR_cols_ones, AR_cols))
else:
AR_cols = AR_cols_ones
G_cols = self._develop_gain(basis, squared=False, log=False)
return (AR_cols, G_cols)
def _develop_parcor(self, LAR, basis):
single_dim = LAR.ndim == 1
LAR = np.atleast_2d(LAR)
# n_basis, n_epochs, n_points = basis.shape
# ordar, n_basis = LAR.shape
# ordar, n_epochs, n_points = lar_list.shape
lar_list = np.tensordot(LAR, basis, axes=([1], [0]))
if single_dim:
# n_epochs, n_points = lar_list.shape
lar_list = lar_list[0, :, :]
return lar_list
# ------------------------------------------------ #
# Functions that should be overloaded #
# ------------------------------------------------ #
def tensordotcontract(tensors):
'''
Use np.tensordot to implement the contraction of tensors.
Parameters
----------
tensors : list of Tensor
The tensors to be contracted.
Returns
-------
Tensor
The contracted tensor.
'''
def _contract_(a,b):
common=set(a.labels)&set(b.labels)
aaxes=[a.axis(label) for label in common]
baxes=[b.axis(label) for label in common]
labels=[label for label in it.chain(a.labels,b.labels) if label not in common]
axes=(aaxes,baxes) if len(common)>0 else 0
return Tensor(np.tensordot(a,b,axes=axes),labels=labels)
result=tensors[0]
for i in xrange(1,len(tensors)):
result=_contract_(result,tensors[i])
return result
def mgf_kmesh(self,omega,kmesh):
'''
Returns the mesh of the Green's functions in the mixed representation with respect to momentums.
Parameters
----------
omega : np.complex128/np.complex64
The frequency of the mixed Green's functions.
kmesh : (n+1)d ndarray like
The kmesh of the mixed Green's functions.
And n is the spatial dimension of the system.
Returns
-------
3d ndarray
The mesh of the mixed Green's functions.
'''
cgf=self.cgf(omega)
return np.tensordot(cgf,inv(np.identity(cgf.shape[0],dtype=cgf.dtype)-self.pt_kmesh(kmesh).dot(cgf)),axes=([1],[1])).transpose((1,0,2))
def VCAEB(engine,app):
'''
This method calculates the single particle spectrum along a path in the Brillouin zone.
'''
engine.rundependences(app.name)
engine.cache.pop('pt_kmesh',None)
erange,kmesh,nk=np.linspace(app.emin,app.emax,app.ne),app.path.mesh('k'),app.path.rank('k')
result=np.zeros((nk,app.ne,3))
result[:,:,0]=np.tensordot(np.array(xrange(nk)),np.ones(app.ne),axes=0)
result[:,:,1]=np.tensordot(np.ones(nk),erange,axes=0)
for i,omega in enumerate(erange):
result[:,i,2]=-(np.trace(engine.gf_kmesh(omega+app.mu+app.eta*1j,kmesh),axis1=1,axis2=2)).imag/engine.nopt/np.pi
name='%s_%s'%(engine,app.name)
if app.savedata: np.savetxt('%s/%s.dat'%(engine.dout,name),result.reshape((nk*app.ne,3)))
if app.plot: app.figure('P',result,'%s/%s'%(engine.dout,name))
if app.returndata: return result
def update_grad_input(self, x, grad_output):
self.grad_input = np.zeros_like(x)
if x.dims == 3:
C, H, W = x.shape
x_flat = x.view(C, H * W)
self.buffer = np.dot(grad_output, x_flat)
self.buffer += np.dot(grad_output.T, x_flat)
self.grad_input = self.buffer.view(C, H, W)
if x.dims == 4:
N, C, H, W = x.shape
x_flat = x.view(N, C, H * W)
self.buffer = np.tensordot(grad_output, x_flat, 2)
self.buffer += np.tensordot(grad_output.transpose(2, 3), x_flat, 2)
self.grad_input = self.buffer.view(N, C, H, W)
if self.normalize:
self.buffer /= C * H * W
return self.grad_input
def test_against_numpy_tensordot(self):
""" Test against numpy.tensordot in 2D case """
stream = tuple(np.random.random((8, 8)) for _ in range(2))
for axis in (0, 1, 2):
with self.subTest('axis = {}'.format(axis)):
from_numpy = np.tensordot(*stream)
from_stream = last(itensordot(stream))
self.assertSequenceEqual(from_numpy.shape, from_stream.shape)
self.assertTrue(np.allclose(from_numpy, from_stream))
def test_against_numpy_inner(self):
""" Test against numpy.tensordot in 2D case """
stream = tuple(np.random.random((8, 8)) for _ in range(2))
for axis in (0, 1, 2):
with self.subTest('axis = {}'.format(axis)):
from_numpy = np.inner(*stream)
from_stream = last(iinner(stream))
self.assertSequenceEqual(from_numpy.shape, from_stream.shape)
self.assertTrue(np.allclose(from_numpy, from_stream))
def _eig_rightvec_add(rightvec, mpo_lten, mps_lten):
"""Add one column to the right vector.
:param rightvec: existing right vector
It has three indices: mps bond, mpo bond, complex conjugate mps bond
:param op_lten: Local tensor of the MPO
:param mps_lten: Local tensor of the current MPS eigenstate
This does the same thing as _eig_leftvec_add(), except that
'left' and 'right' are exchanged in the contractions (but not in
the axis names of the input tensors).
"""
rightvec_names = ('mps_bond', 'mpo_bond', 'cc_mps_bond')
mpo_names = ('left_mpo_bond', 'phys_row', 'phys_col', 'right_mpo_bond')
mps_names = ('left_mps_bond', 'phys', 'right_mps_bond')
rightvec = named_ndarray(rightvec, rightvec_names)
mpo_lten = named_ndarray(mpo_lten, mpo_names)
mps_lten = named_ndarray(mps_lten, mps_names)
contract_mps = (('mps_bond', 'right_mps_bond'),)
rightvec = rightvec.tensordot(mps_lten, contract_mps)
rename_mps = (('left_mps_bond', 'mps_bond'),)
rightvec = rightvec.rename(rename_mps)
contract_mpo = (
('mpo_bond', 'right_mpo_bond'),
('phys', 'phys_col'))
rightvec = rightvec.tensordot(mpo_lten, contract_mpo)
contract_cc_mps = (
('cc_mps_bond', 'right_mps_bond'),
('phys_row', 'phys'))
rightvec = rightvec.tensordot(mps_lten.conj(), contract_cc_mps)
rename_mps_mpo = (
('left_mpo_bond', 'mpo_bond'),
('left_mps_bond', 'cc_mps_bond'))
rightvec = rightvec.rename(rename_mps_mpo)
rightvec = rightvec.to_array(rightvec_names)
return rightvec
def _eig_leftvec_add_mps(lv, lt1, lt2):
"""Add one column to the left vector (MPS version)"""
# MPS 1: Interpreted as |psiXpsi| part of the operator
# MPS 2: The current eigvectector candidate
# NB: It would be more efficient to store lt1.conj() instead of lt1.
# lv axes: 0: mps1 bond, 1: mps2 bond
lv = np.tensordot(lv, lt1.conj(), axes=(0, 0))
# lv axes: 0: mps2 bond, 1: physical leg, 2: mps1 bond
lv = np.tensordot(lv, lt2, axes=((0, 1), (0, 1)))
# lv axes: 0: mps1 bond, 1: mps2 bond
return lv
def _eig_rightvec_add_mps(rv, lt1, lt2):
"""Add one column to the right vector (MPS version)"""
# rv axes: 0: mps1 bond, 1: mps2 bond
rv = np.tensordot(rv, lt1.conj(), axes=(0, 2))
# rv axes: 0: mps2 bond, 1: mps1 bond, 2: physical leg
rv = np.tensordot(rv, lt2, axes=((0, 2), (2, 1)))
# rv axes: 0: mps1 bond, 1: mps2 bond
return rv
def dot(mpa1, mpa2, axes=(-1, 0), astype=None):
"""Compute the matrix product representation of the contraction of ``a``
and ``b`` over the given axes. [:ref:`Sch11 <Sch11>`, Sec. 4.2]
:param mpa1, mpa2: Factors as MPArrays
:param axes: Tuple ``(ax1, ax2)`` where ``ax1`` (``ax2``) is a single
physical leg number or sequence of physical leg numbers
referring to ``mpa1`` (``mpa2``). The first (second, etc) entries
of ``ax1`` and ``ax2`` will be contracted. Very similar to the
``axes`` argument for :func:`numpy.tensordot()`.
(default: ``(-1, 0)``)
.. note:: Note that the default value of ``axes`` is different compared to
:func:`numpy.tensordot`.
:param astype: Return type. If ``None``, use the type of ``mpa1``
:returns: Dot product of the physical arrays
"""
assert len(mpa1) == len(mpa2), \
"Length is not equal: {} != {}".format(len(mpa1), len(mpa2))
# adapt the axes from physical to true legs
if isinstance(axes[0], collections.Sequence):
axes = tuple(tuple(ax + 1 if ax >= 0 else ax - 1 for ax in axes2)
for axes2 in axes)
else:
axes = tuple(ax + 1 if ax >= 0 else ax - 1 for ax in axes)
ltens = [_local_dot(l, r, axes) for l, r in zip(mpa1.lt, mpa2.lt)]
if astype is None:
astype = type(mpa1)
return astype(ltens)
def _local_dot(ltens_l, ltens_r, axes):
"""Computes the local tensors of a dot product `dot(l, r)`.
Besides computing the normal dot product, this function rearranges the
virtual legs in such a way that the result is a valid local tensor again.
:param ltens_l: Array with `ndim > 1`
:param ltens_r: Array with `ndim > 1`
:param axes: Axes to compute dot product using the convention of
:func:`numpy.tensordot()`. Note that these correspond to the true
(and not the just the physical) legs of the local tensors
:returns: Correct local tensor representation
"""
# number of contracted legs need to be the same
clegs_l = len(axes[0]) if isinstance(axes[0], collections.Sequence) else 1
clegs_r = len(axes[1]) if isinstance(axes[0], collections.Sequence) else 1
assert clegs_l == clegs_r, \
"Number of contracted legs differ: {} != {}".format(clegs_l, clegs_r)
res = np.tensordot(ltens_l, ltens_r, axes=axes)
# Rearrange the virtual-dimension legs
res = np.rollaxis(res, ltens_l.ndim - clegs_l, 1)
res = np.rollaxis(res, ltens_l.ndim - clegs_l,
ltens_l.ndim + ltens_r.ndim - clegs_l - clegs_r - 1)
return res.reshape((ltens_l.shape[0] * ltens_r.shape[0], ) +
res.shape[2:-2] +
(ltens_l.shape[-1] * ltens_r.shape[-1],))
def _adapt_to_add_r(rightvec, compr_lten, tgt_lten):
"""Add one column to the right vector.
:param rightvec: existing right vector
It has two indices: `compr_mps_bond` and `tgt_mps_bond`
:param compr_lten: Local tensor of the compressed MPS
:param tgt_lten: Local tensor of the target MPS
Construct R from [:ref:`Sch11 <Sch11>`, Fig. 27, p. 48]. See comments in
:func:`_adapt_to_add_l()` for further details.
.. todo:: Adapt tensor leg names.
"""
rightvec_names = ('compr_bond', 'tgt_bond')
compr_names = ('compr_left_bond', 'compr_phys', 'compr_right_bond')
tgt_names = ('tgt_left_bond', 'tgt_phys', 'tgt_right_bond')
rightvec = named_ndarray(rightvec, rightvec_names)
compr_lten = named_ndarray(compr_lten, compr_names)
tgt_lten = named_ndarray(tgt_lten, tgt_names)
contract_compr_mps = (('compr_bond', 'compr_right_bond'),)
rightvec = rightvec.tensordot(compr_lten, contract_compr_mps)
contract_tgt_mps = (
('compr_phys', 'tgt_phys'),
('tgt_bond', 'tgt_right_bond'))
rightvec = rightvec.tensordot(tgt_lten.conj(), contract_tgt_mps)
rename = (
('compr_left_bond', 'compr_bond'),
('tgt_left_bond', 'tgt_bond'))
rightvec = rightvec.rename(rename)
rightvec = rightvec.to_array(rightvec_names)
return rightvec
def matdot(A, B, axes=((-1,), (0,))):
"""np.tensordot with sane defaults for matrix multiplication"""
return np.tensordot(A, B, axes=axes)
def pmps_dm_to_array(pmps, global_=False):
"""Convert PMPS to full array representation of the density matrix
The runtime of this method scales with D**3 instead of D**6 where
D is the rank and D**6 is the scaling of using :func:`pmps_to_mpo`
and :func:`to_array`. This is useful for obtaining reduced states
of a PMPS on non-consecutive sites, as normalizing before using
:func:`pmps_to_mpo` may not be sufficient to reduce the rank
in that case.
.. note:: The resulting array will have dimension-1 physical legs removed.
"""
out = np.ones((1, 1, 1))
# Axes: 0 phys, 1 upper rank, 2 lower rank
for lt in pmps.lt:
out = np.tensordot(out, lt, axes=(1, 0))
# Axes: 0 phys, 1 lower rank, 2 phys, 3 anc, 4 upper rank
out = np.tensordot(out, lt.conj(), axes=((1, 3), (0, 2)))
# Axes: 0 phys, 1 phys, 2 upper rank, 3 phys, 4 lower rank
out = np.rollaxis(out, 3, 2)
# Axes: 0 phys, 1 phys, 2 phys, 3 upper bound, 4 lower rank
out = out.reshape((-1, out.shape[3], out.shape[4]))
# Axes: 0 phys, 1 upper rank, 2 lower rank
out_shape = [dim for dim, _ in pmps.shape for rep in (1, 2) if dim > 1]
out = out.reshape(out_shape)
if global_:
assert len(set(out_shape)) == 1
out = local_to_global(out, sites=len(out_shape) // 2)
return out
def test_chain(nr_sites, local_dim, rank, rgen, dtype):
# This test produces at most `nr_sites` by tensoring two
# MPOs. This doesn't work for :code:`nr_sites = 1`.
if nr_sites < 2:
return
# NOTE: Everything here is in local form!!!
mpo = factory.random_mpa(nr_sites // 2, (local_dim, local_dim), rank,
randstate=rgen, dtype=dtype)
op = mpo.to_array()
# Test with 2-factors with full form
mpo_double = mp.chain((mpo, mpo))
op_double = np.tensordot(op, op, axes=(tuple(), ) * 2)
assert len(mpo_double) == 2 * len(mpo)
assert_array_almost_equal(op_double, mpo_double.to_array())
assert_array_equal(mpo_double.ranks, mpo.ranks + (1,) + mpo.ranks)
assert mpo.dtype == dtype
# Test 3-factors iteratively (since full form would be too large!!
diff = mp.chain((mpo, mpo, mpo)) - mp.chain((mpo, mp.chain((mpo, mpo))))
diff.canonicalize()
assert len(diff) == 3 * len(mpo)
assert mp.norm(diff) < 1e-6
# local_dim, rank
