r"""This module contains a base class for a Matrix Product State (MPS).

An MPS looks roughly like this::

|   -- B[0] -- B[1] -- B[2] -- ...
|       |       |      |

We use the following label convention for the B (where arrows indicate qconj)::

|  vL ->- B ->- vR
|         |
|         ^
|         p

We store one 3-leg tensor _B[i] with labels 'vL', 'vR', 'p' for each of the L sites
0 <= i < L.
Additionally, we store L+1 singular value arrays _S[ib] on each bond 0 <= ib <= L,
independent of the boundary conditions.
_S[ib] gives the singlur values on the bond i-1, i.
However, be aware that e.g. :attr:~tenpy.networks.mps.MPS.chi returns only the dimensions of the
:attr:~tenpy.networks.mps.MPS.nontrivial_bonds depending on the boundary conditions.

The matrices and singular values always represent a normalized state
(i.e. np.linalg.norm(psi._S[ib]) == 1 up to roundoff errors),
but (for finite MPS) we keep track of the norm in :attr:~tenpy.networks.mps.MPS.norm
(which is respected by :meth:~tenpy.networks.mps.MPS.overlap, ...).

Valid MPS boundary conditions (not to confuse with bc_coupling of
:class:tenpy.models.model.CouplingModel)  are the following:

==========  ===================================================================================
bc        description
==========  ===================================================================================
'finite'    Finite MPS, G0 s1 G1 ... s{L-1} G{l-1}. This is acchieved
by using a trivial left and right bond s[0] = s[-1] = np.array([1.]).
'segment'   Generalization of 'finite', describes an MPS embedded in left and right
environments. The left environment is described by chi[0] *orthonormal* states
which are weighted by the singular values s[0]. Similar, s[L] weight some
right orthonormal states. You can think of the left and right states to be
generated by additional MPS, such that the overall structure is something like
... s L s L [s0 G0 s1 G1 ... s{L-1} G{L-1} s{L}] R s R s R ...
(where we save the part in the brackets [ ... ] ).
'infinite'  infinite MPS (iMPS): we save a 'MPS unit cell' [s0 G0 s1 G1 ... s{L-1} G{L-1}]
which is repeated periodically, identifying all indices modulo self.L.
In particular, the last bond L is identified with 0.
(The MPS unit cell can differ from a lattice unit cell).
bond is identified with the first one.
==========  ===================================================================================

An MPS can be in different 'canonical forms' (see [Vidal2004]_, [Schollwoeck2011]_).
To take care of the different canonical forms, algorithms should use functions like
:meth:~tenpy.networks.mps.MPS.get_theta, :meth:~tenpy.networks.mps.MPS.get_B
and :meth:~tenpy.networks.mps.MPS.set_B instead of accessing them directly,
as they return the B in the desired form (which can be chosen as an argument).
The values of the tuples for the form correspond to the exponent of the singular values
on the left and right.
To keep track of a "mixed" canonical form A A A s B B, we save the tuples for each
site of the MPS in :attr:MPS.form.

======== ========== ==========================================================================
form   tuple      description
======== ========== ==========================================================================
'B'  (0, 1)     right canonical: _B[i] = -- Gamma[i] -- s[i+1]--
The default form, which algorithms asssume.
'C'  (0.5, 0.5) symmetric form: _B[i] = -- s[i]**0.5 -- Gamma[i] -- s[i+1]**0.5--
'A'  (1, 0)     left canonical: _B[i] = -- s[i] -- Gamma[i] --.
'G'  (0, 0)     Save only _B[i] = -- Gamma[i] --.
'Th' (1, 1)     Form of a local wave function theta with singular value on both sides.
psi.get_B(i, 'Th') is equivalent to psi.get_theta(i, n=1).
None None   General non-canoncial form.
Valid form for initialization, but you need to call
:meth:~tenpy.networks.mps.MPS.canonical_form (or similar)
before using algorithms.
======== ========== ==========================================================================
"""
# Copyright 2018-2020 TeNPy Developers, GNU GPLv3

import numpy as np
import warnings
import random
from functools import reduce
import scipy.sparse.linalg.eigen.arpack

from ..linalg import np_conserved as npc
from ..linalg import sparse
from .site import GroupedSite, group_sites
from ..tools.misc import to_iterable, argsort, to_array
from ..tools.math import lcm, speigs, entropy
from ..tools.params import asConfig
from ..algorithms.truncation import TruncationError, svd_theta

__all__ = ['MPS', 'MPSEnvironment', 'TransferMatrix', 'build_initial_state']

class MPS:
r"""A Matrix Product State, finite (MPS) or infinite (iMPS).

Parameters
----------
sites : list of :class:~tenpy.networks.site.Site
Defines the local Hilbert space for each site.
Bs : list of :class:~tenpy.linalg.np_conserved.Array
The 'matrices' of the MPS. Labels are vL, vR, p (in any order).
SVs : list of 1D array
The singular values on *each* bond. Should always have length L+1.
Entries out of :attr:nontrivial_bonds are ignored.
bc : 'finite' | 'segment' | 'infinite'
Boundary conditions as described in the tabel of the module doc-string.
form : (list of) {'B' | 'A' | 'C' | 'G' | 'Th' | None | tuple(float, float)}
The form of the stored 'matrices', see table in module doc-string.
A single choice holds for all of the entries.

Attributes
----------
sites : list of :class:~tenpy.networks.site.Site
Defines the local Hilbert space for each site.
bc : {'finite', 'segment', 'infinite'}
Boundary conditions as described in above table.
form : list of {None | tuple(float, float)}
Describes the canonical form on each site.
None means non-canonical form.
For form = (nuL, nuR), the stored _B[i] are
s**form[0] -- Gamma -- s**form[1] (in Vidal's notation).
chinfo : :class:~tenpy.linalg.np_conserved.ChargeInfo
The nature of the charge.
dtype : type
The data type of the _B.
norm : float
The norm of the state, i.e. sqrt(<psi|psi>).
Ignored for (normalized) :meth:expectation_value, but important for :meth:overlap.
grouped : int
Number of sites grouped together, see :meth:group_sites.
_B : list of :class:npc.Array
The 'matrices' of the MPS. Labels are vL, vR, p (in any order).
We recommend using :meth:get_B and :meth:set_B, which will take care of the different
canonical forms.
_S : list of (None | 1D array)
The singular values on each virtual bond, length L+1.
May be None if the MPS is not in canonical form.
Otherwise, _S[i] is to the left of _B[i].
We recommend using :meth:get_SL, :meth:get_SR, :meth:set_SL, :meth:set_SR, which
takes proper care of the boundary conditions.
_valid_forms : dict
Class attribute.
Mapping for canonical forms to a tuple (nuL, nuR) indicating that
self._Bs[i] = s[i]**nuL -- Gamma[i] -- s[i]**nuR is saved.
_valid_bc : tuple of str
Class attribute. Possible valid boundary conditions.
_transfermatrix_keep : int
How many states to keep at least when diagonalizing a :class:TransferMatrix.
Important if the state develops a near-degeneracy.
_p_label, _B_labels : list of str
Class attribute. _p_label defines the physical legs of the B-tensors, _B_labels lists
all the labels of the B tensors. Used by methods like :meth:get_theta to avoid
the necessity of re-implementations for derived classes like the
:class:~tenpy.networks.purification_mps.Purification_MPS if just the number of physical
legs changed.
"""
# Canonical form conventions: the saved B = s**nu[0]--Gamma--s**nu[1].
# For the canonical forms, nu[0] + nu[1] = 1
_valid_forms = {
'A': (1., 0.),
'C': (0.5, 0.5),
'B': (0., 1.),
'G': (0., 0.),  # like Vidal's Gamma.
'Th': (1., 1.),
None: None,  # means 'not in any canonical form'
}

# valid boundary conditions. Don't overwrite this!
_valid_bc = ('finite', 'segment', 'infinite')
# the "physical" labels for each B
_p_label = ['p']
# All labels of each tensor in _B (order is used!)
_B_labels = ['vL', 'p', 'vR']

def __init__(self, sites, Bs, SVs, bc='finite', form='B', norm=1.):
self.sites = list(sites)
self.chinfo = self.sites[0].leg.chinfo
self.dtype = dtype = np.find_common_type([B.dtype for B in Bs], [])
self.form = self._parse_form(form)
self.bc = bc  # one of 'finite', 'periodic', 'segment'.
self.norm = norm
self.grouped = 1

# make copies of Bs and SVs
self._B = [B.astype(dtype, copy=True).itranspose(self._B_labels) for B in Bs]
self._S = [None] * (self.L + 1)
for i in range(self.L + 1)[self.nontrivial_bonds]:
self._S[i] = np.array(SVs[i], dtype=np.float)
if self.bc == 'infinite':
self._S[-1] = self._S[0]
elif self.bc == 'finite':
self._S[0] = self._S[-1] = np.ones([1])
self._transfermatrix_keep = 1
self.test_sanity()

def test_sanity(self):
"""Sanity check, raises ValueErrors, if something is wrong."""
if self.bc not in self._valid_bc:
raise ValueError("invalid boundary condition: " + repr(self.bc))
if len(self._B) != self.L:
raise ValueError("wrong len of self._B")
if len(self._S) != self.L + 1:
raise ValueError("wrong len of self._S")
for i, B in enumerate(self._B):
if B.get_leg_labels() != self._B_labels:
raise ValueError("B has wrong labels {0!r}, expected {1!r}".format(
B.get_leg_labels(), self._B_labels))
if self._S[i].shape[-1] != B.get_leg('vL').ind_len or \
self._S[i+1].shape[0] != B.get_leg('vR').ind_len:
raise ValueError("shape of B incompatible with len of singular values")
if not self.finite or i + 1 < self.L:
B2 = self._B[(i + 1) % self.L]
B.get_leg('vR').test_contractible(B2.get_leg('vL'))
if self.bc == 'finite':
if len(self._S[0]) != 1 or len(self._S[-1]) != 1:
raise ValueError("non-trivial outer bonds for finite MPS")
elif self.bc == 'infinite':
if np.any(self._S[self.L] != self._S[0]):
raise ValueError("iMPS with S[0] != S[L]")
assert len(self.form) == self.L
for f in self.form:
if f is not None:
assert isinstance(f, tuple)
assert len(f) == 2

def copy(self):
"""Returns a copy of self.

The copy still shares the sites, chinfo, and LegCharges of the B tensors, but the values of
B and S are deeply copied.
"""
# __init__ makes deep copies of B, S
cp = self.__class__(self.sites, self._B, self._S, self.bc, self.form, self.norm)
cp.grouped = self.grouped
cp._transfermatrix_keep = self._transfermatrix_keep
return cp

def save_hdf5(self, hdf5_saver, h5gr, subpath):
"""Export self into a HDF5 file.

This method saves all the data it needs to reconstruct self with :meth:from_hdf5.

Specifically, it saves
:attr:sites,
:attr:chinfo (under these names),
:attr:_B as "tensors",
:attr:_S as "singular_values",
:attr:bc as "boundary_condition", and
:attr:form converted to a single array of shape (L, 2) as "canonical_form",
Moreover, it saves :attr:norm, :attr:L, :attr:grouped and
:attr:_transfermatrix_keep (as "transfermatrix_keep") as HDF5 attributes, as well as
the maximum of :attr:chi under the name "max_bond_dimension".

Parameters
----------
hdf5_saver : :class:~tenpy.tools.hdf5_io.Hdf5Saver
Instance of the saving engine.
h5gr : :classGroup
HDF5 group which is supposed to represent self.
subpath : str
The name of h5gr with a '/' in the end.
"""
hdf5_saver.save(self.sites, subpath + "sites")
hdf5_saver.save(self._B, subpath + "tensors")
hdf5_saver.save(self._S, subpath + "singular_values")
hdf5_saver.save(self.bc, subpath + "boundary_condition")
hdf5_saver.save(np.array(self.form), subpath + "canonical_form")
hdf5_saver.save(self.chinfo, subpath + "chinfo")
h5gr.attrs["norm"] = self.norm
h5gr.attrs["grouped"] = self.grouped
h5gr.attrs["transfermatrix_keep"] = self._transfermatrix_keep
h5gr.attrs["max_bond_dimension"] = np.max(self.chi)  # same

@classmethod
"""Load instance from a HDF5 file.

This method reconstructs a class instance from the data saved with :meth:save_hdf5.

Parameters
----------
hdf5_loader : :class:~tenpy.tools.hdf5_io.Hdf5Loader
h5gr : :class:Group
HDF5 group which is represent the object to be constructed.
subpath : str
The name of h5gr with a '/' in the end.

Returns
-------
obj : cls
Newly generated class instance containing the required data.
"""
obj = cls.__new__(cls)  # create class instance, no __init__() call

obj.form = [tuple(f) for f in form]

obj.dtype = np.find_common_type([B.dtype for B in obj._B], [])
obj.test_sanity()
return obj

@classmethod
def from_lat_product_state(cls, lat, p_state, **kwargs):
"""Construct an MPS from a product state given in lattice coordinates.

This is a wrapper around :meth:from_product_state.
The purpuse is to make the p_state argument independent of the order of the Lattice,
and specify it in terms of lattice indices instead.

Parameters
----------
lat : :class:~tenpy.models.lattice.Lattice
The underlying lattice defining the geometry and Hilbert Space.
p_state : array_like of {int | str | 1D array}
Defines the product state to be represented.
Should be of dimension lat.dim+1, entries are indexed by lattice indices.
Entries of the array as for the p_state argument of :meth:from_product_state.
It gets tiled to the shape lat.shape, if it is smaller.
**kwargs :
Other keyword arguments as definied in :meth:from_product_state.
bc is set by default from lat.bc_MPS.

Returns
-------
product_mps : :class:MPS
An MPS representing the specified product state.

Examples
--------
Let's first consider a :class:~tenpy.models.lattice.Ladder composed of a
:class:~tenpy.networks.site.SpinHalfSite and a
:class:~tenpy.networks.site.FermionSite.

>>> spin_half = tenpy.networks.site.SpinHalfSite()
>>> fermion = tenpy.networks.site.FermionSite()

To initialize a state of up-spins on the spin sites and half-filled ferions, you can use:

>>> p_state = [["up", "empty"], ["up", "full"]]

Note that the same p_state works for a finite lattice of even length, say L=10, as
well. We then just "tile" in x-direction, i.e., repeat the specified state 5 times:

You can also easily half-fill a :class:~tenpy.models.lattice.Honeycomb, for example
with only the A sites occupied, or as stripe parallel to the x-direction (stripe_x,
alternating along y axis),
or as stripes parallel to the y-direction (stripe_y, alternating along x axis).

>>> honeycomb = tenpy.models.lattice.Honeycomb([4, 4], [fermion, fermion], bc_MPS="finite")
>>> p_state_only_A = [[["empty", "full"]]]
>>> psi_only_A = tenpy.networks.MPS.from_lat_product_state(honeycomb, p_state_only_A)
>>> p_state_stripe_x = [[["empty", "empty"],
...                      ["full", "full"]]]
>>> psi_stripe_x = tenpy.networks.MPS.from_lat_product_state(honeycomb, p_state_stripe_x)
>>> p_state_stripe_y = [[["empty", "empty"]],
...                      [["full", "full"]]]
>>> psi_stripe_y = tenpy.networks.MPS.from_lat_product_state(honeycomb, p_state_stripe_y)
"""
kwargs.setdefault("bc", lat.bc_MPS)
p_state = np.array(p_state, dtype=object)
if p_state.ndim == len(lat.shape):  # == lat.dim + 1
p_state = to_array(p_state, shape=lat.shape)  # tile to lattice shape
p_state_flat = p_state[tuple(lat.order.T)]  # "advanced" numpy indexing
elif p_state.ndim == len(lat.shape) + 1:
# extra dimension could be from purely 1D array entries
# make sure this is the case by converting to float
p_state = np.array(p_state, kwargs.get("dtype", np.float64))
# tile to lattice shape, ignore last dimension
p_state = to_array(p_state, shape=lat.shape + (None, ))
inds = tuple(lat.order.T) + (slice(None), )
p_state_flat = p_state[inds]  # "advanced" numpy indexing
else:
raise ValueError("wrong dimension of p_state. Expected {0:d}-dimensional array of "
"(string, int, or 1D array)".format(d=lat.dim + 1))
return cls.from_product_state(lat.mps_sites(), p_state_flat, **kwargs)

@classmethod
def from_product_state(cls,
sites,
p_state,
bc='finite',
dtype=np.float64,
permute=True,
form='B',
chargeL=None):
"""Construct a matrix product state from a given product state.

Parameters
----------
sites : list of :class:~tenpy.networks.site.Site
The sites defining the local Hilbert space.
p_state : list of {int | str | 1D array}
Defines the product state to be represented; one entry for each site of the MPS.
An entry of str type is translated to an int with the help of
:meth:~tenpy.networks.site.Site.state_labels.
An entry of int type represents the physical index of the state to be used.
An entry which is a 1D array defines the complete wavefunction on that site; this
allows to make a (local) superposition.
bc : {'infinite', 'finite', 'segmemt'}
MPS boundary conditions. See docstring of :class:MPS.
dtype : type or string
The data type of the array entries.
permute : bool
The :class:~tenpy.networks.Site might permute the local basis states if charge
conservation gets enabled.
If permute is True (default), we permute the given p_state locally according to
each site's :attr:~tenpy.networks.Site.perm.
The p_state entries should then always be given as if conserve=None in the Site.
form : (list of) {'B' | 'A' | 'C' | 'G' | None | tuple(float, float)}
Defines the canonical form. See module doc-string.
A single choice holds for all of the entries.
chargeL : charges
Leg charges at bond 0, which are purely conventional.

Returns
-------
product_mps : :class:MPS
An MPS representing the specified product state.

Examples
--------
Example to get a Neel state for a :class:~tenpy.models.tf_ising.TIChain:

>>> M = TFIChain({'L': 10})
>>> p_state = ["up", "down"] * (L//2)  # repeats entries L/2 times
>>> psi = MPS.from_product_state(M.lat.mps_sites(), p_state, bc=M.lat.bc_MPS)

The meaning of the labels "up","down" is defined by the :class:~tenpy.networks.Site,
in this example a :class:~tenpy.networks.site.SpinHalfSite.

Extending the example, we can replace the spin in the center with one with arbitrary
angles theta, phi in the bloch sphere:

>>> M = TFIChain({'L': 8, 'conserve': None})
>>> p_state = ["up", "down"] * (L//2)  # repeats entries L/2 times
>>> bloch_sphere_state = np.array([np.cos(theta/2), np.exp(1.j*phi)*np.sin(theta/2)])
>>> p_state[L//2] = bloch_sphere_state   # replace one spin in center
>>> psi = MPS.from_product_state(M.lat.mps_sites(), p_state, bc=M.lat.bc_MPS, dtype=np.complex)

Note that for the more general :class:~tenpy.models.spins.SpinChain,
the order of the two entries for the bloch_sphere_state would be *exactly the opposite*
(when we keep the the north-pole of the bloch sphere being the up-state).
The reason is that the SpinChain uses the general :class:~tenpy.networks.site.SpinSite,
where the states are orderd ascending from 'down' to 'up'.
The :class:~tenpy.networks.site.SpinHalfSite on the other hand uses the order
'up', 'down' where that the Pauli matrices look as usual.

Moreover, note that you can not write this bloch state (for theta != 0, pi) when
conserving symmetries, as the two physical basis states correspond to different symmetry
sectors.
"""
sites = list(sites)
L = len(sites)
p_state = list(p_state)
if len(p_state) != L:
raise ValueError("Length of p_state does not match number of sites.")
ci = sites[0].leg.chinfo
Bs = []
chargeL = ci.make_valid(chargeL)  # sets to zero if None
legL = npc.LegCharge.from_qflat(ci, [chargeL])  # (no need to bunch)
for p_st, site in zip(p_state, sites):
perm = permute
if isinstance(p_st, str):
p_st = site.state_labels[p_st]  # translate labels into "int"
perm = False
try:
iter(p_st)
except TypeError:
# just an int for p_st
B = np.zeros((site.dim, 1, 1), dtype)
B[p_st, 0, 0] = 1.0
else:  # iter works
if len(p_st) != site.dim:
raise ValueError("p_state incompatible with local dim:" + repr(p_st))
B = np.array(p_st, dtype).reshape((site.dim, 1, 1))
if perm:
B = B[site.perm, :, :]
Bs.append(B)
SVs = [[1.]] * (L + 1)
return cls.from_Bflat(sites, Bs, SVs, bc, dtype, False, form, legL)

@classmethod
def from_Bflat(cls,
sites,
Bflat,
SVs=None,
bc='finite',
dtype=None,
permute=True,
form='B',
legL=None):
"""Construct a matrix product state from a set of numpy arrays Bflat and singular vals.

Parameters
----------
sites : list of :class:~tenpy.networks.site.Site
The sites defining the local Hilbert space.
Bflat : iterable of numpy ndarrays
The matrix defining the MPS on each site, with legs 'p', 'vL', 'vR'
(physical, virtual left/right).
SVs : list of 1D array | None
The singular values on *each* bond. Should always have length L+1.
By default (None), set all singular values to the same value.
Entries out of :attr:nontrivial_bonds are ignored.
bc : {'infinite', 'finite', 'segmemt'}
MPS boundary conditions. See docstring of :class:MPS.
dtype : type or string
The data type of the array entries. Defaults to the common dtype of Bflat.
permute : bool
The :class:~tenpy.networks.Site might permute the local basis states if charge
conservation gets enabled.
If permute is True (default), we permute the given Bflat locally according to
each site's :attr:~tenpy.networks.Site.perm.
The p_state argument should then always be given as if conserve=None in the Site.
form : (list of) {'B' | 'A' | 'C' | 'G' | None | tuple(float, float)}
Defines the canonical form of Bflat. See module doc-string.
A single choice holds for all of the entries.
leg_L : LegCharge | None
Leg charges at bond 0, which are purely conventional.
If None, use trivial charges.

Returns
-------
mps : :class:MPS
An MPS with the matrices Bflat converted to npc arrays.
"""
sites = list(sites)
L = len(sites)
Bflat = list(Bflat)
if len(Bflat) != L:
raise ValueError("Length of Bflat does not match number of sites.")
ci = sites[0].leg.chinfo
if legL is None:
legL = npc.LegCharge.from_qflat(ci, [ci.make_valid(None)] * Bflat[0].shape[1])
legL = legL.bunch()[1]
if SVs is None:
SVs = [np.ones(B.shape[1]) / np.sqrt(B.shape[1]) for B in Bflat]
SVs.append(np.ones(Bflat[-1].shape[2]) / np.sqrt(Bflat[-1].shape[2]))
Bs = []
if dtype is None:
dtype = np.dtype(np.common_type(*Bflat))
for i, site in enumerate(sites):
B = np.array(Bflat[i], dtype)
if permute:
B = B[site.perm, :, :]
# calculate the LegCharge of the right leg
legs = [site.leg, legL, None]  # other legs are known
legs = npc.detect_legcharge(B, ci, legs, None, qconj=-1)
B = npc.Array.from_ndarray(B, legs, dtype)
B.iset_leg_labels(['p', 'vL', 'vR'])
Bs.append(B)
legL = legs[-1].conj()  # prepare for next i
if bc == 'infinite':
# for an iMPS, the last leg has to match the first one.
# so we need to gauge qtotal of the last B such that the right leg matches.
chdiff = Bs[-1].get_leg('vR').charges[0] - Bs[0].get_leg('vL').charges[0]
Bs[-1] = Bs[-1].gauge_total_charge('vR', ci.make_valid(chdiff))
return cls(sites, Bs, SVs, form=form, bc=bc)

@classmethod
def from_full(cls,
sites,
psi,
form=None,
cutoff=1.e-16,
normalize=True,
bc='finite',
outer_S=None):
"""Construct an MPS from a single tensor psi with one leg per physical site.

Performs a sequence of SVDs of psi to split off the B matrices and obtain the singular
values, the result will be in canonical form.
Obviously, this is only well-defined for finite or segment boundary conditions.

Parameters
----------
sites : list of :class:~tenpy.networks.site.Site
The sites defining the local Hilbert space.
psi : :class:~tenpy.linalg.np_conserved.Array
The full wave function to be represented as an MPS.
Should have labels 'p0', 'p1', ...,  'p{L-1}'.
Additionally, it may have (or must have for 'segment' bc) the legs 'vL', 'vR',
which are trivial for 'finite' bc.
form  : 'B' | 'A' | 'C' | 'G' | None
The canonical form of the resulting MPS, see module doc-string.
None defaults to 'A' form on the first site and 'B' form on all following sites.
cutoff : float
Cutoff of singular values used in the SVDs.
normalize : bool
Whether the resulting MPS should have 'norm' 1.
bc : 'finite' | 'segment'
Boundary conditions.
outer_S : None | (array, array)
For 'semgent' bc the singular values on the left and right of the considered segment,
None for 'finite' boundary conditions.

Returns
-------
psi_mps : :class:MPS
MPS representation of psi, in canonical form and possibly normalized.
"""
if form is not None and form not in ['B', 'A', 'C', 'G']:
raise ValueError("Invalid form: " + repr(form))
if bc != 'finite' and bc != 'segment':
raise ValueError("Wrong boundary conditions: " + repr(bc))
# perform SVDs to bring it into 'B' form, afterwards change the form.
L = len(sites)
assert (L >= 2)
B_list = [None] * L
S_list = [None] * (L + 1)
norm = 1. if normalize else npc.norm(psi)
if not psi.has_label('vL'):
elif bc == 'finite' and psi.get_leg('vL').ind_len != 1:
raise ValueError("non-trivial left leg for 'finite' bc!")
if not psi.has_label('vR'):
elif bc == 'finite' and psi.get_leg('vR').ind_len != 1:
raise ValueError("non-trivial left leg for 'finite' bc!")
labels = ['vL'] + ['p' + str(i) for i in range(L)] + ['vR']
psi.itranspose(labels)
# combine legs from left
for i in range(0, L - 1):
psi = psi.combine_legs([0, 1])  # combines the legs until i
# now psi has only three legs: '(((vL.p0).p1)...p{L-2})', 'p{L-1}', 'vR'
for i in range(L - 1, 0, -1):
# split off B[i]
psi = psi.combine_legs([labels[i + 1], 'vR'])
psi, S, B = npc.svd(psi, inner_labels=['vR', 'vL'], cutoff=cutoff)
S /= np.linalg.norm(S)  # normalize
if i > 1:
psi.iscale_axis(S, 1)
B_list[i] = B.split_legs(1).replace_label(labels[i + 1], 'p')
S_list[i] = S
psi = psi.split_legs(0)
# psi is now the first B in 'A' form
B_list[0] = psi.replace_label(labels[1], 'p')
B_form = ['A'] + ['B'] * (L - 1)
if bc == 'finite':
S_list[0] = S_list[-1] = np.ones([1], dtype=np.float)
elif outer_S is not None:
S_list[0], S_list[-1] = outer_S
res = cls(sites, B_list, S_list, bc=bc, form=B_form, norm=norm)
if form is not None:
res.convert_form(form)
return res

@classmethod
def from_singlets(cls,
site,
L,
pairs,
up='up',
down='down',
lonely=[],
lonely_state='up',
bc='finite'):
"""Create an MPS of entangled singlets.

Parameters
----------
site : :class:~tenpy.networks.site.Site
The site defining the local Hilbert space, taken uniformly for all sites.
L : int
The number of sites.
pairs : list of (int, int)
Pairs of sites to be entangled; the returned MPS will have a singlet
for each pair in pairs.
up, down : int | str
A singlet is defined as (|up down> - |down up>)/2**0.5,
up and down give state indices or labels defined on the corresponding site.
lonely : list of int
Sites which are not included into a singlet pair.
lonely_state : int | str
The state for the lonely sites.
bc : {'infinite', 'finite', 'segmemt'}
MPS boundary conditions. See docstring of :class:MPS.

Returns
-------
singlet_mps : :class:MPS
An MPS representing singlets on the specified pairs of sites.
"""
# sort each pair s.t. i < j
pairs = [((i, j) if i < j else (j, i)) for (i, j) in pairs]
# sort by smaller site of the pair
pairs.sort(key=lambda x: x[0])
pairs.append((L, L))
lonely = sorted(lonely) + [L]
# generate building block tensors
up = site.state_index(up)
down = site.state_index(down)
lonely_state = site.state_index(lonely_state)
Close = np.zeros([site.dim, site.dim], dtype=np.float_)
Close[up, down] = 1.
Close[down, up] = -1.
Close = npc.Array.from_ndarray(Close, [site.leg, site.leg])  # no conj() !
Id = npc.eye_like(Close, 0)
Lonely = np.zeros(site.dim, dtype=np.float_)
Lonely[lonely_state] = 1
Lonely = npc.Array.from_ndarray(Lonely, [site.leg])
Bs = []
Ss = [np.ones(1)]
forms = []
open_singlets = []  # the k-th open singlet should be closed at site open_singlets[k]
Ts = []  # the tensors on the current site
labels_L = []
for i in range(L):
labels_R = labels_L[:]
next_Ts = Ts[:]
if i == pairs[0][0]:  # open a new singlet
j = pairs[0][1]
lbl = 's{0:d}-{1:d}'.format(i, j)
pairs.pop(0)
open_singlets.append(j)
next_Ts.append(Id.copy().iset_leg_labels([lbl + 'L', lbl]))
Open.iset_leg_labels(['p', lbl])
Ts.append(Open.copy(deep=False))
labels_R.append(lbl)
forms.append('A')
elif i == lonely[0]:  # just a lonely state
Ts.append(Lonely)
lonely.pop(0)
forms.append('B')
else:  # close a singlet
k = open_singlets.index(i)
Close.iset_leg_labels([labels_L[k] + 'L', 'p'])
Ts[k] = Close
next_Ts.pop(k)
open_singlets.pop(k)
labels_R.pop(k)
forms.append('B')
# generate B from Ts
B = reduce(npc.outer, Ts)
labels_L = [lbl_ + 'L' for lbl_ in labels_L]
if len(labels_L) > 0 and len(labels_R) > 0:
B = B.combine_legs([labels_L, labels_R], new_axes=[0, 2], qconj=[+1, -1])
B.iset_leg_labels(['vL', 'p', 'vR'])
elif len(labels_L) == 0 and len(labels_R) == 0:
B.iset_leg_labels(['vL', 'p', 'vR'])
elif len(labels_L) == 0:
B = B.combine_legs([labels_R], new_axes=[1], qconj=[-1])
B.iset_leg_labels(['p', 'vR'])
else:  # len(labels_R) == 0
B = B.combine_legs([labels_L], new_axes=[0], qconj=[+1])
B.iset_leg_labels(['vL', 'p'])
Bs.append(B)
N = 2**len(labels_R)
Ss.append(np.ones(N) / (N**0.5))
Ts = next_Ts
labels_L = labels_R
return cls([site] * L, Bs, Ss, bc=bc, form=forms)

@property
def L(self):
"""Number of physical sites; for an iMPS the len of the MPS unit cell."""
return len(self.sites)

@property
def dim(self):
"""List of local physical dimensions."""
return [site.dim for site in self.sites]

@property
def finite(self):
"""Distinguish MPS vs iMPS.

True for an MPS (bc='finite', 'segment'), False for an iMPS (bc='infinite').
"""
assert (self.bc in self._valid_bc)
return self.bc != 'infinite'

@property
def chi(self):
"""Dimensions of the (nontrivial) virtual bonds."""
# s.shape[0] == len(s) for 1D numpy array, but works also for a 2D npc Array.
return [min(s.shape) for s in self._S[self.nontrivial_bonds]]

@property
def nontrivial_bonds(self):
"""Slice of the non-trivial bond indices, depending on self.bc."""
if self.bc == 'finite':
return slice(1, self.L)
elif self.bc == 'segment':
return slice(0, self.L + 1)
elif self.bc == 'infinite':
return slice(0, self.L)

def get_B(self, i, form='B', copy=False, cutoff=1.e-16, label_p=None):
"""Return (view of) B at site i in canonical form.

Parameters
----------
i : int
Index choosing the site.
form : 'B' | 'A' | 'C' | 'G' | 'Th' | None | tuple(float, float)
The (canonical) form of the returned B.
For None, return the matrix in whatever form it is.
If any of the tuple entry is None, also don't scale on the corresponding axis.
copy : bool
Whether to return a copy even if form matches the current form.
cutoff : float
During DMRG with a mixer, S may be a matrix for which we need the inverse.
This is calculated as the Penrose pseudo-inverse, which uses a cutoff for the
singular values.
label_p : None | str
Ignored by default (None).
Otherwise replace the physical label 'p' with 'p'+label_p'.
(For derived classes with more than one "physical" leg, replace all the physical leg
labels accordingly.)

Returns
-------
B : :class:~tenpy.linalg.np_conserved.Array
The MPS 'matrix' B at site i with leg labels 'vL', 'p', 'vR'.
May be a view of the matrix (if copy=False),
or a copy (if the form changed or copy=True).

Raises
------
ValueError : if self is not in canoncial form and form is not None.
"""
i = self._to_valid_index(i)
new_form = self._to_valid_form(form)
old_form = self.form[i]
B = self._B[i]
if copy:
B = B.copy()
if new_form is not None and old_form != new_form:
if old_form is None:
raise ValueError("can't convert form of non-canonical state!")
if new_form[0] is not None and new_form[0] - old_form[0] != 0.:
B = self._scale_axis_B(B, self.get_SL(i), new_form[0] - old_form[0], 'vL', cutoff)
if new_form[1] is not None and new_form[1] - old_form[1] != 0.:
B = self._scale_axis_B(B, self.get_SR(i), new_form[1] - old_form[1], 'vR', cutoff)
if label_p is not None:
B = self._replace_p_label(B, label_p)
return B

def set_B(self, i, B, form='B'):
"""Set B at site i.

Parameters
----------
i : int
Index choosing the site.
B : :class:~tenpy.linalg.np_conserved.Array
The 'matrix' at site i. No copy is made!
Should have leg labels 'vL', 'p', 'vR' (not necessarily in that order).
form : 'B' | 'A' | 'C' | 'G' | 'Th' | None | tuple(float, float)
The (canonical) form of the B to set.
None stands for non-canonical form.
"""
i = self._to_valid_index(i)
self.form[i] = self._to_valid_form(form)
self.dtype = np.find_common_type([self.dtype, B.dtype], [])
self._B[i] = B.itranspose(self._B_labels)

def set_svd_theta(self, i, theta, trunc_par=None, update_norm=False):
"""SVD a two-site wave function theta and save it in self.

Parameters
----------
i : int
theta is the wave function on sites i, i+1.
theta : :class:~tenpy.linalg.np_conserved.Array
The two-site wave function with labels combined into "(vL.p0)", "(p1.vR)",
trunc_par : None | dict
Parameters for truncation, see :cfg:config:truncation.
If None, no truncation is done.
update_norm : bool
If True, multiply the norm of theta into :attr:norm.
"""
i0 = self._to_valid_index(i)
i1 = self._to_valid_index(i0 + 1)
self.dtype = np.find_common_type([self.dtype, theta.dtype], [])
qtotal_LR = [self._B[i0].qtotal, None]
if trunc_par is None:
U, S, VH = npc.svd(theta, qtotal_LR=qtotal_LR, inner_labels=['vR', 'vL'])
renorm = np.linalg.norm(S)
S /= renorm
err = None
if update_norm:
self.norm *= renorm
else:
U, S, VH, err, renorm = svd_theta(theta, trunc_par, qtotal_LR)
if update_norm:
self.norm *= renorm
U = U.split_legs().ireplace_label('p0', 'p')
VH = VH.split_legs().ireplace_label('p1', 'p')
self._B[i0] = U.itranspose(self._B_labels)
self.form[i0] = self._valid_forms['A']
self._B[i1] = VH.itranspose(self._B_labels)
self.form[i1] = self._valid_forms['B']
self.set_SR(i, S)
return err

def get_SL(self, i):
"""Return singular values on the left of site i"""
i = self._to_valid_index(i)
return self._S[i]

def get_SR(self, i):
"""Return singular values on the right of site i"""
i = self._to_valid_index(i)
return self._S[i + 1]

def set_SL(self, i, S):
"""Set singular values on the left of site i"""
i = self._to_valid_index(i)
self._S[i] = S
if not self.finite and i == 0:
self._S[self.L] = S

def set_SR(self, i, S):
"""Set singular values on the right of site i"""
i = self._to_valid_index(i)
self._S[i + 1] = S
if not self.finite and i == self.L - 1:
self._S[0] = S

def get_op(self, op_list, i):
"""Given a list of operators, select the one corresponding to site i.

Parameters
----------
op_list : (list of) {str | npc.array}
List of operators from which we choose. We assume that op_list[j] acts on site
j. If the length is shorter than L, we repeat it periodically.
Strings are translated using :meth:~tenpy.networks.site.Site.get_op of site i.
i : int
Index of the site on which the operator acts.

Returns
-------
op : npc.array
One of the entries in op_list, not copied.
"""
if self.finite and i > self.L or i < 0:
raise ValueError("i = {0:d} out of bounds for finite MPS".format(i))
op = op_list[i % len(op_list)]
if (isinstance(op, str)):
op = self.sites[i % self.L].get_op(op)
return op

def get_theta(self, i, n=2, cutoff=1.e-16, formL=1., formR=1.):
"""Calculates the n-site wavefunction on sites[i:i+n].

Parameters
----------
i : int
Site index.
n : int
Number of sites. The result lives on sites[i:i+n].
cutoff : float
During DMRG with a mixer, S may be a matrix for which we need the inverse.
This is calculated as the Penrose pseudo-inverse, which uses a cutoff for the
singular values.
formL : float
Exponent for the singular values to the left.
formR : float
Exponent for the singular values to the right.

Returns
-------
theta : :class:~tenpy.linalg.np_conserved.Array
The n-site wave function with leg labels vL, p0, p1, .... p{n-1}, vR.
In Vidal's notation (with s=lambda, G=Gamma):
theta = s**form_L G_i s G_{i+1} s ... G_{i+n-1} s**form_R.
"""
i = self._to_valid_index(i)
for j in range(i, i + n):
if self.form[j % self.L] is None:
raise ValueError("can't calculate theta for non-canonical form")
if n == 1:
return self.get_B(i, (1., 1.), True, cutoff, '0')
elif n < 1:
raise ValueError("n needs to be larger than 0")
# n >= 2: contract some B's
theta = self.get_B(i, (formL, None), False, cutoff, '0')  # right form as stored
_, old_fR = self.form[i]
for k in range(1, n):  # non-empty range
j = self._to_valid_index(i + k)
new_fR = None if k + 1 < n else formR  # right form as stored, except for last B
B = self.get_B(j, (1. - old_fR, new_fR), False, cutoff, str(k))
_, old_fR = self.form[j]
theta = npc.tensordot(theta, B, axes=['vR', 'vL'])
return theta

def convert_form(self, new_form='B'):
"""Tranform self into different canonical form (by scaling the legs with singular values).

Parameters
----------
new_form : (list of) {'B' | 'A' | 'C' | 'G' | 'Th' | None | tuple(float, float)}
The form the stored 'matrices'. The table in module doc-string.
A single choice holds for all of the entries.

Raises
------
ValueError : if trying to convert from a None form. Use :meth:canonical_form instead!
"""
new_forms = self._parse_form(new_form)
for i, new_form in enumerate(new_forms):
new_B = self.get_B(i, form=new_form, copy=False)  # calculates the desired form.
self.set_B(i, new_B, form=new_form)

def increase_L(self, new_L=None):
"""Modify self inplace to enlarge the MPS unit cell; in place.

.. deprecated:: 0.5.1
This method will be removed in version 1.0.0.
Use the equivalent psi.enlarge_mps_unit_cell(new_L//psi.L) instead of
psi.increase_L(new_L).

Parameters
----------
new_L : int
New number of sites. Needs to be an integer multiple of :attr:L.
Defaults to 2*self.L.
"""
old_L = self.L
if new_L is None:
new_L = 2 * old_L
if new_L % old_L:
raise ValueError("new_L = {0:d} not a multiple of old L={1:d}".format(new_L, old_L))
factor = new_L // old_L
warnings.warn(
"use psi.enlarge_mps_unit_cell(factor=new_L//psi.L) "
"instead of psi.increase_L(new_L).", FutureWarning, 2)
self.enlarge_mps_unit_cell(factor)

def enlarge_mps_unit_cell(self, factor=2):
"""Repeat the unit cell for infinite MPS boundary conditions; in place.

Parameters
----------
factor : int
The new number of sites in the unit cell will be increased from L to factor*L.
"""
if int(factor) != factor:
raise ValueError("factor should be integer!")
if factor <= 1:
raise ValueError("can't shrink!")
if self.bc == 'segment':
raise ValueError("can't enlarge segment MPS")
self.sites = factor * self.sites
self._B = factor * self._B
self._S = factor * self._S[:-1] + [self._S[-1]]
self.form = factor * self.form
self.test_sanity()

def roll_mps_unit_cell(self, shift=1):
"""Shift the section we define as unit cellof an infinite MPS; in place.

Suppose we have a unit cell with tensors [A, B, C, D] (repeated on both sites).
With shift = 1, the new unit cell will be [D, A, B, C],
whereas shift = -1 will give [B, C, D, A].

Parameters
----------
shift : int
By how many sites to move the tensors to the right.
"""
if self.finite:
raise ValueError("makes only sense for infinite boundary conditions")
inds = np.roll(np.arange(self.L), shift)
self.sites = [self.sites[i] for i in inds]
self.form = [self.form[i] for i in inds]
self._B = [self._B[i] for i in inds]
self._S = [self._S[i] for i in inds]
self._S.append(self._S[0])

def group_sites(self, n=2, grouped_sites=None):
"""Modify self inplace to group sites.

Group each n sites together using the :class:~tenpy.networks.site.GroupedSite.
This might allow to do TEBD with a Trotter decomposition,
or help the convergence of DMRG (in case of too long range interactions).

Parameters
----------
n : int
Number of sites to be grouped together.
grouped_sites : None | list of :class:~tenpy.networks.site.GroupedSite
The sites grouped together.

--------
group_split : Reverts the grouping.
"""
self.convert_form('B')
if grouped_sites is None:
grouped_sites = group_sites(self.sites, n, charges='same')
else:
assert grouped_sites[0].n_sites == n
Bs = []
Ss = []
i = 0
B_form = self._valid_forms['B']
for gs in grouped_sites:
n_sites = gs.n_sites
new_B = self.get_theta(i, gs.n_sites, formL=B_form[0], formR=B_form[1])
comb_legs = [[lbl + str(k) for k in range(n_sites)] for lbl in self._p_label]
# comb_legs = [['p0', 'p1', ... ]] for usual MPS
axes = list(range(1, 1 + len(self._p_label)))  # [1]
new_B = new_B.combine_legs(comb_legs, new_axes=axes, qconj=[+1] * len(axes))
new_B.legs[1].test_equal(gs.leg)  # test legcharge compatibility
Bs.append(new_B.iset_leg_labels(self._B_labels))  # ['vL', 'p', 'vR']
Ss.append(self._S[i])
i += n_sites
Ss.append(self._S[-1])  # right-most singular values: need L+1 entries
self._B = Bs
self._S = Ss
self.sites = grouped_sites
self.form = [B_form] * len(grouped_sites)
self.grouped = self.grouped * n

def group_split(self, trunc_par=None):
"""Modify self inplace to split previously grouped sites.

Parameters
----------
trunc_par : dict
Parameters for truncation, see :cfg:config:truncation.
Defaults to {'chi_max': max(self.chi)}.

Returns
-------
trunc_err : :class:~tenpy.algorithms.truncation.TruncationError
The error introduced by the truncation for the splitting.

--------
group_sites : Should have been used before to combine sites.
"""
if trunc_par is None:
trunc_par = {}
self.convert_form('B')
if self.L > 1:
trunc_par.setdefault('chi_max', max(self.chi))
n0 = self.sites[0].n_sites
sites = []
Bs = []
Ss = []
trunc_err = TruncationError()
for i, gs in enumerate(self.sites):
sites.extend(gs.sites)
n = gs.n_sites
Ss_new = []
Bs_new = []
B_gr = self.get_B(i, form='B').transpose(self._B_labels)  # vL, p, vR
B_gr.idrop_labels(self._p_label)  # avoid warning: split label not called '(...)'
n_p_label = len(self._p_label)
split_legs = list(range(1, 1 + n_p_label))
transp = [i for k in range(n) for i in range(1 + k, 1 + n * n_p_label, n)]
transp = ['vL'] + transp + ['vR']
B_gr = B_gr.split_legs(split_legs).itranspose(transp)
theta = self.get_theta(i, n=1)
theta.idrop_labels(self._get_p_label('0'))  # avoid warning
theta = theta.split_legs(split_legs).itranspose(transp)
# for usual MPS, B_gr and theta have legs vL p0 p1 ... p{n-1} vR
# for MPS with legs p, q, they have legs vL p0 q0 p1 q1 ... q{-n-1} vR
combine = [list(range(B_gr.rank - n_p_label - 1)), list(range(-n_p_label - 1, 0, +1))]
# combine = [[0, 1, .... {n-1}], [-2, -1]] for usual MPS
axes_contr = [combine[1], list(range(1, 2 + n_p_label))]
for j in range(n - 1, 0, -1):
# split off the right-most physical leg and vR from theta
# theta: vL p0 ... pj vR
theta = theta.combine_legs(combine, qconj=[+1, -1])
U, S, V, err, _ = svd_theta(theta, trunc_par, inner_labels=['vR', 'vL'])
Ss_new.append(S)
trunc_err += err
theta = U.split_legs(0)  # vL p0 ... pj-1 vR
for _ in range(n_p_label):
combine[0].pop()
B = V.split_legs(1).iset_leg_labels(self._B_labels)  # vL p vR
B_gr = npc.tensordot(B_gr, B.conj(), axes=axes_contr)  # vL p0 ... pj-1 vR
Bs_new.append(B)
Bs_new.append(B_gr.iset_leg_labels(self._B_labels))  # inversion free :)
Ss_new.append(self.get_SL(i))
Bs.extend(Bs_new[::-1])
Ss.extend(Ss_new[::-1])
Ss.append(self._S[-1])
self.sites = sites
self._B = Bs
self._S = Ss
self.grouped = self.grouped // n0
self.form = [self._valid_forms['B']] * len(sites)
self.test_sanity()
return trunc_err

def get_grouped_mps(self, blocklen):
r"""contract blocklen subsequent tensors into a single one and return result as a new MPS.

blocklen = number of subsequent sites to be combined.

Returns
-------
new MPS object with bunched sites.
"""
groupedMPS = self.copy()
groupedMPS.group_sites(n=blocklen)
return groupedMPS

def get_total_charge(self, only_physical_legs=False):
"""Calculate and return the qtotal of the whole MPS (when contracted).

Parameters
----------
only_physical_legs : bool
For 'finite' boundary conditions, the total charge can be gauged away
by changing the LegCharge of the trivial legs on the left and right of the MPS.
This option allows to project out the trivial legs to get the actual "physical"
total charge.

Returns
-------
qtotal : charges
The sum of the qtotal of the individual B tensors.
"""
qtotal = np.sum([B.qtotal for B in self._B], axis=0)
if only_physical_legs:
if self.bc != 'finite':
raise ValueError("only_physical_legs not supported for bc=" + repr(self.bc))
qtotal -= self._B[0].get_leg('vL').get_charge(0)
qtotal -= self._B[-1].get_leg('vR').get_charge(0)  # takes qconj into account
return self.chinfo.make_valid(qtotal)

def gauge_total_charge(self, qtotal=None, vL_leg=None, vR_leg=None):
"""Gauge the legcharges of the virtual bonds such that the MPS has a total qtotal.

Parameters
----------
qtotal : (list of) charges
If a single set of charges is given, it is the desired total charge of the MPS
(which :meth:get_total_charge will return afterwards).
By default (None), use 0 charges, unless vL_leg and vR_leg are specified, in which
case we adjust the total charge to match these legs.
vL_leg : None | LegCharge
Desired new virtual leg on the very left. Needs to have the same block strucuture as
current leg, but can have shifted charge entries.
vR_leg : None | LegCharge
Desired new virtual leg on the very right. Needs to have the same block strucuture as
current leg, but can have shifted charge entries.
Should be vL_leg.conj() for infinite MPS, if qtotal is not given.
"""
if self.chinfo.qnumber == 0:
return
if vL_leg is not None:
vL_chdiff = vL_leg.get_charge(0) - self._B[0].get_leg('vL').get_charge(0)
if vR_leg is not None:
vR_chdiff = vR_leg.get_charge(0) - self._B[-1].get_leg('vR').get_charge(0)
if qtotal is None:
if vL_leg is not None and vR_leg is not None:
qtotal = self.get_total_charge() + vL_chdiff + vR_chdiff
qtotal = self.chinfo.make_valid(qtotal)
if qtotal.ndim == 1:
qtotal_factor = np.array([0] * (self.L - 1) + [1], npc.QTYPE)
qtotal = qtotal_factor[:, np.newaxis] * qtotal[np.newaxis, :]
if qtotal.shape != (self.L, self.chinfo.qnumber):
raise ValueError("wrong shape of qtotal")
if vL_leg is not None:
B = self._B[0]
if np.any(vL_chdiff != 0):
self._B[0] = B.gauge_total_charge('vL', B.qtotal + vL_chdiff, vL_leg.qconj)
B.get_leg('vL').test_equal(vL_leg)
for i in range(self.L):
B = self._B[i]
desired_qtotal = qtotal[i]
chdiff = B.qtotal - desired_qtotal
if np.any(chdiff != 0):
self._B[i] = B.gauge_total_charge('vR', desired_qtotal)
if i + 1 != self.L:  # this 'vR' is contracted with the 'vL' of the next B
# so we need to adjust the next B as well
nextB = self._B[i + 1]
self._B[i + 1] = nextB.gauge_total_charge('vL', nextB.qtotal + chdiff)
self._B[i].get_leg('vR').test_contractible(self._B[i + 1].get_leg('vL'))
# just to check
assert np.all(self.get_total_charge() == self.chinfo.make_valid(np.sum(qtotal, 0)))
if vR_leg is not None:
# check that the charges match
self._B[-1].get_leg('vR').test_equal(vR_leg)
if self.bc == 'infinite':
self._B[0].get_leg('vL').test_contractible(self._B[-1].get_leg('vR'))
# done

def entanglement_entropy(self, n=1, bonds=None, for_matrix_S=False):
r"""Calculate the (half-chain) entanglement entropy for all nontrivial bonds.

Consider a bipartition of the sytem into :math:A = \{ j: j <= i_b \} and
:math:B = \{ j: j > i_b\} and the reduced density matrix :math:\rho_A = tr_B(\rho).
The von-Neumann entanglement entropy is defined as
:math:S(A, n=1) = -tr(\rho_A \log(\rho_A)) = S(B, n=1).
The generalization for n != 1, n>0 are the Renyi entropies:
:math:S(A, n) = \frac{1}{1-n} \log(tr(\rho_A^2)) = S(B, n=1)

This function calculates the entropy for a cut at different bonds i, for which the
the eigenvalues of the reduced density matrix :math:\rho_A and :math:\rho_B is given
by the squared schmidt values S of the bond.

Parameters
----------
n : int/float
Selects which entropy to calculate;
n=1 (default) is the ususal von-Neumann entanglement entropy.
bonds : None | (iterable of) int
Selects the bonds at which the entropy should be calculated.
None defaults to range(0, L+1)[self.nontrivial_bonds].
for_matrix_S : bool
Switch calculate the entanglement entropy even if the _S are matrices.
Since :math:O(\chi^3) is expensive compared to the ususal :math:O(\chi),
we raise an error by default.

Returns
-------
entropies : 1D ndarray
Entanglement entropies for half-cuts.
entropies[j] contains the entropy for a cut at bond bonds[j]
(i.e. left to site bonds[j]).
"""
if bonds is None:
nt = self.nontrivial_bonds
bonds = range(nt.start, nt.stop)
res = []
for ib in bonds:
s = self._S[ib]
if len(s.shape) == 1:
res.append(entropy(s**2, n))
else:
if for_matrix_S:
# explicitly calculate Schmidt values by diagonalizing (s^dagger s)
s = npc.eigvalsh(npc.tensordot(s.conj(), s, axes=[0, 0]))
res.append(entropy(s, n))
else:
raise ValueError("entropy with non-diagonal schmidt values")
return np.array(res)

def entanglement_entropy_segment(self, segment=[0], first_site=None, n=1):
r"""Calculate entanglement entropy for general geometry of the bipartition.

This function is similar as :meth:entanglement_entropy,
but for more general geometry of the region A to be a segment of a *few* sites.

This is acchieved by explicitly calculating the reduced density matrix of A
and thus works only for small segments.

Parameters
----------
segment : list of int
Given a first site i, the region A_i is defined to be [i+j for j in segment].
first_site : None | (iterable of) int
Calculate the entropy for segments starting at these sites.
None defaults to range(L-segment[-1]) for finite
or range(L) for infinite boundary conditions.
n : int | float
Selects which entropy to calculate;
n=1 (default) is the ususal von-Neumann entanglement entropy,
otherwise the n-th Renyi entropy.

Returns
-------
entropies : 1D ndarray
entropies[i] contains the entropy for the the region A_i defined above.
"""
# Side-Remark: there is a trick to calculate the entanglement for large regions A_i
# of consecutive sites (in our notation, segment = range(La))
# To get the entanglement entropy, diagonalize:
#     --theta---
#       | | |
#     --theta*--
#  Diagonalization is O(chi^6), compared to O(d^{3*La})
segment = np.sort(segment)
if first_site is None:
if self.finite:
first_site = range(0, self.L - segment[-1])
else:
first_site = range(self.L)
comb_legs = [
self._get_p_labels(len(segment), False),
self._get_p_labels(len(segment), True)
]
res = []
for i0 in first_site:
rho = self.get_rho_segment(segment + i0)
rho = rho.combine_legs(comb_legs, qconj=[+1, -1])
p = npc.eigvalsh(rho)
res.append(entropy(p, n))
return np.array(res)

def entanglement_spectrum(self, by_charge=False):
r"""return entanglement energy spectrum.

Parameters
----------
by_charge : bool
Wheter we should sort the spectrum on each bond by the possible charges.

Returns
-------
ent_spectrum : list
For each (non-trivial) bond the entanglement spectrum.
If by_charge is False, return (for each bond) a sorted 1D ndarray
with the convention :math:S_i^2 = e^{-\xi_i}, where :math:S_i labels a Schmidt
value and :math:\xi_i labels the entanglement 'energy' in the returned spectrum.
If by_charge is True, return a a list of tuples (charge, sub_spectrum)
for each possible charge on that bond.
"""
if by_charge:
res = []
for i in range(self.L + 1)[self.nontrivial_bonds]:
ss = -2. * np.log(self._S[i])
if i < self.L:
leg = self._B[i].get_leg('vL')
else:  # i == L: segment b.c.
leg = self._B[i - 1].get_leg('vR').conj()
spectrum = [(leg.get_charge(qi), np.sort(ss[leg.get_slice(qi)]))
for qi in range(leg.block_number)]
res.append(spectrum)
return res
else:
return [np.sort(-2. * np.log(ss)) for ss in self._S[self.nontrivial_bonds]]

def get_rho_segment(self, segment):
"""Return reduced density matrix for a segment.

Note that the dimension of rho_A scales exponentially in the length of the segment.

Parameters
----------
segment : iterable of int
Sites for which the reduced density matrix is to be calculated.
Assumed to be sorted.

Returns
-------
rho : :class:~tenpy.linalg.np_conserved.Array
Reduced density matrix of the segment sites.
Labels 'p0', 'p1', ..., 'pk', 'p0*', 'p1*', ..., 'pk*' with k=len(segment).
"""
if len(segment) > 12:
warnings.warn("{0:d} sites in the segment, that's much!".format(len(segment)),
stacklevel=2)
if len(segment) > 20:
raise ValueError("too large segment; this is exponentially expensive!")
segment = np.sort(segment)
if np.all(segment[1:] == segment[:-1] + 1):  # consecutive
theta = self.get_theta(segment[0], segment[-1] - segment[0] + 1)
rho = npc.tensordot(theta, theta.conj(), axes=(['vL', 'vR'], ['vL*', 'vR*']))
return rho
rho = self.get_theta(segment[0], 1)
rho = npc.tensordot(rho, rho.conj(), axes=('vL', 'vL*'))
k = 1
contract_axes = (['vR*'] + self._p_label, ['vL*'] + self._get_p_label('*'))
for i in range(segment[0] + 1, segment[-1]):
B = self.get_B(i)
if i == segment[k]:
B = self._replace_p_label(B, str(k))
k += 1
rho = npc.tensordot(rho, B, axes=('vR', 'vL'))
rho = npc.tensordot(rho, B.conj(), axes=('vR*', 'vL*'))
else:
rho = npc.tensordot(rho, B, axes=('vR', 'vL'))
rho = npc.tensordot(rho, B.conj(), axes=contract_axes)
B = self._replace_p_label(self.get_B(segment[-1]), str(k))
rho = npc.tensordot(rho, B, axes=('vR', 'vL'))
rho = npc.tensordot(rho, B.conj(), axes=(['vR*', 'vR'], ['vL*', 'vR*']))
return rho

def probability_per_charge(self, bond=0):
"""Return probabilites of charge value on the left of a given bond.

For example for particle number conservation, define
:math:N_b = sum_{i<b} n_i for a given bond b.
This function returns the possible values of N_b as rows of charge_values,
and for each row the probabilty that this combination occurs in the given state.

Parameters
----------
bond : int
The bond to be considered. The returned charges are summed on the left of this bond.

Returns
-------
charge_values : 2D array
Columns correspond to the different charges in self.chinfo.
Rows are the different charge fluctuations at this bond
probabilities : 1D array
For each row of charge_values the probablity for these values of charge fluctuations.
"""
if self.bc == 'segment' and bond == self.L:
S = self.get_SR(self.L - 1)**2
leg = self.get_B(self.L - 1, form=None).get_leg('vR').conj()
else:  # usually the case
S = self.get_SL(bond)**2
leg = self.get_B(bond, form=None).get_leg('vL')
assert leg.qconj == +1
if not leg.is_blocked():
raise ValueError("leg not blocked: can have duplicate entries in charge values")
ps = []
for qi in range(leg.block_number):
sl = leg.get_slice(qi)
ps.append(np.sum(S[sl]))
ps = np.array(ps)
if abs(np.sum(ps) - 1.) > 1.e-10:
warnings.warn("Probability_per_charge: Sum of probabilites not 1. Canonical form?",
stacklevel=2)
return leg.charges.copy(), ps

def average_charge(self, bond=0):
r"""Return the average charge for the block on the left of a given bond.

For example for particle number conservation, define
:math:N_b = sum_{i<b} n_i for a given bond b.
Then this function returns :math:<\psi| N_b |\psi>.

Parameters
----------
bond : int
The bond to be considered.
The returned charges are summed over the sites left of bond.

Returns
-------
average_charge : 1D array
For each type of charge in :attr:chinfo
the average value when summing the charge values over sites left of the given bond.
"""
charges, ps = self.probability_per_charge(bond)
return np.sum(ps[:, np.newaxis] * charges, axis=0)

def charge_variance(self, bond=0):
r"""Return the charge variance on the left of a given bond.

For example for particle number conservation, define
:math:N_b = sum_{i<b} n_i for a given bond b.
Then this function returns :math:<\psi| N_b^2 |\psi> - (<\psi| N_b |\psi>)^2.

Parameters
----------
bond : int
The bond to be considered.
The returned charges are summed over the sites left of bond.

Returns
-------
average_charge : 1D array
For each type of charge in :attr:chinfo
the variance of of the charge values left of the given bond.
"""
charges_mean = self.average_charge(bond)
charges, ps = self.probability_per_charge(bond)
return np.sum(ps[:, np.newaxis] * (charges - charges_mean[:, np.newaxis])**2, axis=0)

def mutinf_two_site(self, max_range=None, n=1):
"""Calculate the two-site mutual information :math:I(i:j).

Calculates :math:I(i:j) = S(i) + S(j) - S(i,j),
where :math:S(i) is the single site entropy on site :math:i
and :math:S(i,j) the two-site entropy on sites :math:i,j.

Parameters
----------
max_range : int
Maximal distance |i-j| for which the mutual information should be calculated.
None defaults to L-1.
n : float
Selects the entropy to use, see :func:~tenpy.tools.math.entropy.

Returns
-------
coords : 2D array
Coordinates for the mutinf array.
mutinf : 1D array
mutinf[k] is the mutual information :math:I(i:j) between the
sites i, j = coords[k].
"""
#  Basically the code of get_rho_segment and entanglement_entropy,
#  but optimized to run in O(L*max_range)
if max_range is None:
max_range = self.L
S_i = self.entanglement_entropy_segment(n=n)  # single-site entropy
legs_ij = self._get_p_labels(2, False), self._get_p_labels(2, True)
# = (['p0', 'p1'], ['p0*', 'p1*'])
contr_legs = (
['vR*'] + self._get_p_label('1'),  # ['vL', 'p1']
['vL*'] + self._get_p_label('1*'))  # ['vL*', 'p1*']
mutinf = []
coord = []
for i in range(self.L):
rho = self.get_theta(i, 1)
rho = npc.tensordot(rho, rho.conj(), axes=('vL', 'vL*'))
jmax = i + max_range + 1
if self.finite:
jmax = min(jmax, self.L)
for j in range(i + 1, jmax):
B = self._replace_p_label(self.get_B(j, form='B'), '1')  # 'vL', 'vR', 'p1'
rho = npc.tensordot(rho, B, axes=['vR', 'vL'])
rho_ij = npc.tensordot(rho, B.conj(), axes=(['vR*', 'vR'], ['vL*', 'vR*']))
rho_ij = rho_ij.combine_legs(legs_ij, qconj=[+1, -1])
S_ij = entropy(npc.eigvalsh(rho_ij), n)
mutinf.append(S_i[i] + S_i[j % self.L] - S_ij)
coord.append((i, j))
if j + 1 < jmax:
rho = npc.tensordot(rho, B.conj(), axes=contr_legs)
return np.array(coord), np.array(mutinf)

def overlap(self, other, charge_sector=None, ignore_form=False, **kwargs):
"""Compute overlap <self|other>.

Parameters
----------
other : :class:MPS
An MPS with the same physical sites.
charge_sector : None | charges | 0
Selects the charge sector in which the dominant eigenvector of the TransferMatrix is.
None stands for *all* sectors, 0 stands for the sector of zero charges.
If a sector is given, it *assumes* the dominant eigenvector is in that charge sector.
ignore_form : bool
If False (default), take into account the canonical form :attr:form at each site.
If True, we ignore the canonical form (i.e., whether the MPS is in left, right,
mixed or no canonical form) and just contract all the :attr:_B as they are.
(This can give different results!)
**kwargs :
Further keyword arguments given to :meth:TransferMatrix.eigenvectors;
only used for infinite boundary conditions.

Returns
-------
overlap : dtype.type
The contraction <self|other> * self.norm * other.norm
(i.e., taking into account the :attr:norm of both MPS).
For an infinite MPS, <self|other> is the overlap per unit cell, i.e.,
the largest eigenvalue of the TransferMatrix.
"""
if self.finite:
if ignore_form:
# Use TransferMatrix with option to ignore the form
TM = TransferMatrix(self, other, charge_sector=charge_sector, form=None)
res = TM.matvec(TM.initial_guess(1.))  # apply transfer matrix to identity
return npc.trace(res, 0, 1) * self.norm * other.norm
else:
env = MPSEnvironment(self, other)
return env.full_contraction(0)
else:  # infinite
form = None if ignore_form else 'B'
TM = TransferMatrix(self, other, charge_sector=charge_sector, form=form)
ov, _ = TM.eigenvectors(**kwargs)
return ov[0] * self.norm * other.norm

def expectation_value(self, ops, sites=None, axes=None):
"""Expectation value <psi|ops|psi>/<psi|psi> of (n-site) operator(s).

Given the MPS in canonical form, it calculates n-site expectation values.
For example the contraction for a two-site (n = 2) operator on site i would look like::

|          .--S--B[i]--B[i+1]--.
|          |     |     |       |
|          |     |-----|       |
|          |     | op  |       |
|          |     |-----|       |
|          |     |     |       |
|          .--S--B*[i]-B*[i+1]-.

Parameters
----------
ops : (list of) { :class:~tenpy.linalg.np_conserved.Array | str }
The operators, for wich the expectation value should be taken,
All operators should all have the same number of legs (namely 2 n).
If less than self.L operators are given, we repeat them periodically.
Strings (like 'Id', 'Sz') are translated into single-site operators defined by
:attr:sites.
sites : None | list of int
List of site indices. Expectation values are evaluated there.
If None (default), the entire chain is taken (clipping for finite b.c.)
axes : None | (list of str, list of str)
Two lists of each n leg labels giving the physical legs of the operator used for
contraction. The first n legs are contracted with conjugated B,
the second n legs with the non-conjugated B.
None defaults to (['p'], ['p*']) for single site operators (n = 1), or
(['p0', 'p1', ... 'p{n-1}'], ['p0*', 'p1*', .... 'p{n-1}*']) for n > 1.

Returns
-------
exp_vals : 1D ndarray
Expectation values, exp_vals[i] = <psi|ops[i]|psi>, where ops[i] acts on
site(s) j, j+1, ..., j+{n-1} with j=sites[i].

Examples
--------
One site examples (n=1):

>>> psi.expectation_value('Sz')
[Sz0, Sz1, ..., Sz{L-1}]
>>> psi.expectation_value(['Sz', 'Sx'])
[Sz0, Sx1, Sz2, Sx3, ... ]
>>> psi.expectation_value('Sz', sites=[0, 3, 4])
[Sz0, Sz3, Sz4]

Two site example (n=2), assuming homogeneous sites:

>>> SzSx = npc.outer(psi.sites[0].Sz.replace_labels(['p', 'p*'], ['p0', 'p0*']),
psi.sites[1].Sx.replace_labels(['p', 'p*'], ['p1', 'p1*']))
>>> psi.expectation_value(SzSx)
[Sz0Sx1, Sz1Sx2, Sz2Sx3, ... ]   # with len L-1 for finite bc, or L for infinite

Example measuring <psi|SzSx|psi2> on each second site, for inhomogeneous sites:

>>> SzSx_list = [npc.outer(psi.sites[i].Sz.replace_labels(['p', 'p*'], ['p0', 'p0*']),
psi.sites[i+1].Sx.replace_labels(['p', 'p*'], ['p1', 'p1*']))
for i in range(0, psi.L-1, 2)]
>>> psi.expectation_value(SzSx_list, range(0, psi.L-1, 2))
[Sz0Sx1, Sz2Sx3, Sz4Sx5, ...]
"""
ops, sites, n, (op_ax_p, op_ax_pstar) = self._expectation_value_args(ops, sites, axes)
ax_p = ['p' + str(k) for k in range(n)]
ax_pstar = ['p' + str(k) + '*' for k in range(n)]
E = []
for i in sites:
op = self.get_op(ops, i)
op = op.replace_labels(op_ax_p + op_ax_pstar, ax_p + ax_pstar)
theta = self.get_theta(i, n)
C = npc.tensordot(op, theta, axes=[ax_pstar, ax_p])  # C has same labels as theta
E.append(npc.inner(theta, C, axes='labels', do_conj=True))
return np.real_if_close(np.array(E))

def expectation_value_term(self, term, autoJW=True):
r"""Expectation value  <psi|op_{i0}op_{i1}...op_{iN}|psi>/<psi|psi>.

Calculates the expectation value of a tensor product of single-site operators
acting on different sites i0, i1, ... (not necessarily next to each other).
In other words, evaluate the expectation value of a term op0_i0 op1_i1 op2_i2 ....

For example the contraction of three one-site operators on sites i0,
i1=i0+1, i2=i0+3 would look like::

|          .--S--B[i0]---B[i0+1]--B[i0+2]--B[i0+3]--.
|          |     |       |        |        |        |
|          |    op1     op2       |       op3       |
|          |     |       |        |        |        |
|          .--S--B*[i0]--B*[i0+1]-B*[i0+2]-B*[i0+3]-.

Parameters
----------
term : list of (str, int)
List of tuples op, i where i is the MPS index of the site the operator
named op acts on.
The order inside term determines the order in which they act
(in the mathematical convention: the last operator in term is right-most,
so it acts first on a Ket).
autoJW : bool
If True (default), automatically insert Jordan Wigner strings for Fermions as needed.

Returns
-------
exp_val : float/complex
The expectation value of the tensorproduct of the given onsite operators,
<psi|op_i0 op_i1 ... op_iN |psi>/<psi|psi>,
where |psi> is the represented MPS.

--------
correlation_function : efficient way to evaluate many correlation functions.

Examples
--------
>>> a = psi.expectation_value_term([('Sx', 2), ('Sz', 4)])
>>> b = psi.expectation_value_term([('Sz', 4), ('Sx', 2)])
>>> c = psi.expectation_value_multi_sites(['Sz', 'Id', 'Sz'], i0=2)
>>> assert a == b == c
"""
# strategy: translate term into a list "ops" to be used for expectation_value_multi_sites
term = list(term)
i_min = min([t[1] for t in term])
i_max = max([t[1] for t in term])
ops = [None] * (i_max - i_min + 1)
count_JW = 0
for op, i in term:
j = i - i_min  # index in ops
if ops[j] is not None:
ops[j] = ops[j] + " " + op
else:
ops[j] = op
if autoJW and self.sites[self._to_valid_index(i)].op_needs_JW(op):
count_JW += 1
for k in range(j):
if ops[k] is not None:
ops[k] = ops[k] + ' JW'
if ops[k].endswith(' JW JW'):
ops[k] = ops[k][:-len(' JW JW')]
else:
ops[k] = 'JW'
for i in range(len(ops)):
if ops[i] is None:
ops[i] = 'Id'
if count_JW % 2 == 1:
raise ValueError("Odd number of operators which need a Jordan Wigner string")
return self.expectation_value_multi_sites(ops, i_min)

def expectation_value_multi_sites(self, operators, i0):
r"""Expectation value  <psi|op0_{i0}op1_{i0+1}...opN_{i0+N}|psi>/<psi|psi>.

Calculates the expectation value of a tensor product of single-site operators
acting on different sites next to each other.
In other words, evaluate the expectation value of a term
op0_i0 op1_{i0+1} op2_{i0+2} ....

.. warning ::
This function does *not* automatically add Jordan-Wigner strings!
For correct handling of fermions, use :meth:expectation_value_term instead.

Parameters
----------
operators : List of { :class:~tenpy.linalg.np_conserved.Array | str }
List of one-site operators. This method calculates the
expectation value of the n-sites operator given by their tensor
product.
i0 : int
The left most index on which an operator acts, i.e.,
operators[i] acts on site i + i0.

Returns
-------
exp_val : float/complex
The expectation value of the tensorproduct of the given onsite operators,
<psi|operators[0]_{i0} operators[1]_{i0+1} ... |psi>/<psi|psi>,
where |psi> is the represented MPS.
"""
op = operators[0]
if (isinstance(op, str)):
op = self.sites[self._to_valid_index(i0)].get_op(op)
theta = self.get_B(i0, 'Th')
C = npc.tensordot(op, theta, axes=['p*', 'p'])
axes = [['vL*'] + self._get_p_label('*'), ['vL'] + self._p_label]
C = npc.tensordot(theta.conj(), C, axes=axes)
axes[1][0] = 'vR*'
for j in range(1, len(operators)):
op = operators[j]  # the operator
i = i0 + j  # the site it acts on
B = self.get_B(i, form='B')
C = npc.tensordot(C, B, axes=['vR', 'vL'])
if op != 'Id':
if (isinstance(op, str)):
op = self.sites[self._to_valid_index(i)].get_op(op)
C = npc.tensordot(op, C, axes=['p*', 'p'])
C = npc.tensordot(B.conj(), C, axes=axes)
exp_val = npc.trace(C, 'vR*', 'vR')
return np.real_if_close(exp_val)

def expectation_value_terms_sum(self, term_list, prefactors=None):
"""Calculate expectation values for a bunch of terms and sum them up.

This is equivalent to the following expression::

sum([self.expectation_value_term(term)*strength for term, strength in term_list])

However, for effiency, the term_list is converted to an MPO and the expectation value
of the MPO is evaluated.

.. note ::
Due to the way MPO expectation values are evaluated for infinite systems,
it works only if all terms in the term_list start within the MPS unit cell.

.. deprecated:: 0.4.0
prefactor will be removed in version 1.0.0.
Instead, directly give just TermList(term_list, prefactors) as argument.

Parameters
----------
term_list : :class:~tenpy.networks.terms.TermList
The terms and prefactors (strength) to be summed up.
prefactors :
Instead of specifying a :class:~tenpy.networks.terms.TermList,
one can also specify the term_list and strength separately.
This is deprecated.

Returns
-------
terms_sum : list of (complex) float
Equivalent to the expression
sum([self.expectation_value_term(term)*strength for term, strength in term_list]).
_mpo :
Intermediate results: the generated MPO.
For a finite MPS, terms_sum = _mpo.expectation_value(self), for an infinite MPS
terms_sum = _mpo.expectation_value(self) * self.L

--------
expectation_value_term : evaluates a single term.
tenpy.networks.mpo.MPO.expectation_value : expectation value density of an MPO.
"""
from . import mpo, terms
if prefactors is not None:
warnings.warn(
"Deprecated argument prefactors: replace arguments with "
"TermList(term_list, prefactors).", FutureWarning, 2)
term_list = terms.TermList(term_list, prefactors)
L = self.L
if not self.finite:
for term in term_list.terms:
if not 0 <= min([i for _, i in term]) < L:
raise ValueError("term doesn't start in MPS unit cell: " + repr(term))
# conversion
ot, ct = term_list.to_OnsiteTerms_CouplingTerms(self.sites)
bc = 'finite' if self.finite else 'infinite'
mpo_graph = mpo.MPOGraph.from_terms(ot, ct, self.sites, bc)
mpo_ = mpo_graph.build_MPO()
terms_sum = mpo_.expectation_value(self, max_range=ct.max_range())
if not self.finite:
terms_sum = terms_sum * self.L
return terms_sum, mpo_

def correlation_function(self,
ops1,
ops2,
sites1=None,
sites2=None,
opstr=None,
str_on_first=True,
hermitian=False,
autoJW=True):
r"""Correlation function  <psi|op1_i op2_j|psi>/<psi|psi> of single site operators.

Given the MPS in canonical form, it calculates 2-site correlation functions.
For examples the contraction for a two-site operator on site i would look like::

|          .--S--B[i]--B[i+1]--...--B[j]---.
|          |     |     |            |      |
|          |     |     |            op2    |
|          |     op1   |            |      |
|          |     |     |            |      |
|          .--S--B*[i]-B*[i+1]-...--B*[j]--.

Onsite terms are taken in the order <psi | op1 op2 | psi>.

If opstr is given and str_on_first=True, it calculates::

|           for i < j                               for i > j
|
|          .--S--B[i]---B[i+1]--...- B[j]---.     .--S--B[j]---B[j+1]--...- B[i]---.
|          |     |      |            |      |     |     |      |            |      |
|          |     opstr  opstr        op2    |     |     op2    |            |      |
|          |     |      |            |      |     |     |      |            |      |
|          |     op1    |            |      |     |     opstr  opstr        op1    |
|          |     |      |            |      |     |     |      |            |      |
|          .--S--B*[i]--B*[i+1]-...- B*[j]--.     .--S--B*[j]--B*[j+1]-...- B*[i]--.

For i==j, no opstr is included.
For str_on_first=False, the opstr on site min(i, j) is always left out.

Strings (like 'Id', 'Sz') in the arguments are translated into single-site
operators defined by the :class:~tenpy.networks.site.Site on which they act.
Each operator should have the two legs 'p', 'p*'.

Parameters
----------
ops1 : (list of) { :class:~tenpy.linalg.np_conserved.Array | str }
First operator of the correlation function (acting after ops2).
If a list is given, ops1[i] acts on site i of the MPS.
ops2 : (list of) { :class:~tenpy.linalg.np_conserved.Array | str }
Second operator of the correlation function (acting before ops1).
If a list is given, ops2[j] acts on site j of the MPS.
sites1 : None | int | list of int
List of site indices i; a single int is translated to range(0, sites1).
None defaults to all sites range(0, L).
Is sorted before use, i.e. the order is ignored.
sites2 : None | int | list of int
List of site indices; a single int is translated to range(0, sites2).
None defaults to all sites range(0, L).
Is sorted before use, i.e. the order is ignored.
opstr : None | (list of) { :class:~tenpy.linalg.np_conserved.Array | str }
Ignored by default (None).
Operator(s) to be inserted between ops1 and ops2.
If less than :attr:L operators are given, we repeat them periodically.
If given as a list, opstr[r] is inserted at site r (independent of sites1 and
sites2).
str_on_first : bool
Whether the opstr is included on the site min(i, j).
Note the order, which is chosen that way to handle fermionic Jordan-Wigner strings
correctly. (In other words: choose str_on_first=True for fermions!)
hermitian : bool
Optimization flag: if sites1 == sites2 and Ops1[i]^\dagger == Ops2[i]
(which is not checked explicitly!), the resulting C[x, y] will be hermitian.
We can use that to avoid calculations, so hermitian=True will run faster.
autoJW : bool
*Ignored* if opstr is given.
If True, auto-determine if a Jordan-Wigner string is needed.
Works only if exclusively strings were used for op1 and op2.

Returns
-------
C : 2D ndarray
The correlation function C[x, y] = <psi|ops1[i] ops2[j]|psi>,
where ops1[i] acts on site i=sites1[x] and ops2[j] on site j=sites2[y].
If opstr is given, it gives (for str_on_first=True):

- For i < j: C[x, y] = <psi|ops1[i] prod_{i <= r < j} opstr[r] ops2[j]|psi>.
- For i > j: C[x, y] = <psi|prod_{j <= r < i} opstr[r] ops1[i] ops2[j]|psi>.
- For i = j: C[x, y] = <psi|ops1[i] ops2[j]|psi>.

The condition <= r is replaced by a strict < r, if str_on_first=False.

Examples
--------
For a spin chain:

>>> psi.correlation_function("A", "B")
[[A0B0,     A0B1, ..., A0B{L-1}],
[A1B0,     A1B1, ..., A1B{L-1]],
...,
[A{L-1}B0, ALB1, ..., A{L-1}B{L-1}],
]

To evaluate the correlation function for a single i, you can use sites1=[i]:

>>> psi.correlation_function("A", "B", [3])
[[A3B0,     A3B1, ..., A3B{L-1}]]

For fermions, it auto-determines that/whether a Jordan Wigner string is needed:

>>> CdC = psi.correlation_function("Cd", "C")  # optionally: use hermitian=True
>>> psi.correlation_function("C", "Cd")[1, 2] == -CdC[1, 2]
True
>>> np.all(np.diag(CdC) == psi.expectation_value("Cd C"))  # "Cd C" is equivalent to "N"
True

--------
expectation_value_term : best for a single combination of i and j.
"""
if opstr is not None:
autoJW = False
ops1, ops2, sites1, sites2, opstr = self._correlation_function_args(
ops1, ops2, sites1, sites2, opstr)
if autoJW and not all([isinstance(op1, str) for op1 in ops1]):
warnings.warn("Non-string operator: can't auto-determine Jordan-Wigner!", stacklevel=2)
autoJW = False
if autoJW:
need_JW = []
for i in sites1:
need_JW.append(self.sites[i % self.L].op_needs_JW(ops1[i % len(ops1)]))
for j in sites2:
need_JW.append(self.sites[j % self.L].op_needs_JW(ops1[j % len(ops1)]))
if any(need_JW):
if not all(need_JW):
raise ValueError("Some, but not any operators need 'JW' string!")
if not str_on_first:
raise ValueError("Need Jordan Wigner string, but str_on_first=False")
opstr = 'JW'
if hermitian and np.any(sites1 != sites2):
warnings.warn("MPS correlation function can't use the hermitian flag", stacklevel=2)
hermitian = False
C = np.empty((len(sites1), len(sites2)), dtype=np.complex)
for x, i in enumerate(sites1):
# j > i
j_gtr = sites2[sites2 > i]
if len(j_gtr) > 0:
C_gtr = self._corr_up_diag(ops1, ops2, i, j_gtr, opstr, str_on_first, True)
C[x, (sites2 > i)] = C_gtr
if hermitian:
C[x + 1:, x] = np.conj(C_gtr)
# j == i
j_eq = sites2[sites2 == i]
if len(j_eq) > 0:
# on-site correlation function
op12 = npc.tensordot(self.get_op(ops1, i), self.get_op(ops2, i), axes=['p*', 'p'])
C[x, (sites2 == i)] = self.expectation_value(op12, i, [['p'], ['p*']])
if not hermitian:
#  j < i
for y, j in enumerate(sites2):
i_gtr = sites1[sites1 > j]
if len(i_gtr) > 0:
C[(sites1 > j), y] = self._corr_up_diag(ops2, ops1, j, i_gtr, opstr,
str_on_first, False)
# exchange ops1 and ops2 : they commute on different sites,
# but we apply opstr after op1 (using the last argument = False)
return np.real_if_close(C)

def norm_test(self):
"""Check that self is in canonical form.

Returns
-------
norm_error: array, shape (L, 2)
For each site the norm error to the left and right.
The error norm_error[i, 0] is defined as the norm-difference between
the following networks::

|   --theta[i]---.       --s[i]--.
|       |        |    vs         |
|   --theta*[i]--.       --s[i]--.

Similarly, norm_errror[i, 1] is the norm-difference of::

|   .--theta[i]---         .--s[i+1]--
|   |    |          vs     |
|   .--theta*[i]--         .--s[i+1]--
"""
err = np.empty((self.L, 2), dtype=np.float)
lbl_R = (self._get_p_label('0') + ['vR'], self._get_p_label('0*') + ['vR*'])
lbl_L = (['vL'] + self._get_p_label('0'), ['vL*'] + self._get_p_label('0*'))
for i in range(self.L):
th = self.get_theta(i, 1)
rho_L = npc.tensordot(th, th.conj(), axes=lbl_R)
S = self.get_SL(i)
if isinstance(S, npc.Array):  # during DMRG with mixer, S may be a 2D npc.Array
if S.rank != 2:
raise ValueError("Expect 2D npc.Array or 1D numpy ndarray")
rho_L2 = npc.tensordot(S, S.conj(), axes=['vR', 'vR*'])
else:
rho_L2 = npc.diag(S**2, rho_L.get_leg('vL'), dtype=rho_L.dtype)
err[i, 0] = npc.norm(rho_L - rho_L2)
rho_R = npc.tensordot(th, th.conj(), axes=lbl_L)
S = self.get_SR(i)
if isinstance(S, npc.Array):
if S.rank != 2:
raise ValueError("Expect 2D npc.Array or 1D numpy ndarray")
rho_R2 = npc.tensordot(S, S.conj(), axes=['vL', 'vL*'])
else:
rho_R2 = npc.diag(S**2, rho_R.get_leg('vR'), dtype=rho_L.dtype)
err[i, 1] = npc.norm(rho_R - rho_R2)
return err

def canonical_form(self, renormalize=True):
"""Bring self into canonical 'B' form, (re-)calculate singular values.

Simply calls :meth:canonical_form_finite or :meth:canonical_form_infinite.
"""
if self.finite:
return self.canonical_form_finite(renormalize)
else:
self.canonical_form_infinite(renormalize)

def canonical_form_finite(self, renormalize=True, cutoff=0.):
"""Bring a finite (or segment) MPS into canonical form (in place).

If any site is in :attr:form None, it does *not* use any of the singular values S
(for 'finite' boundary conditions, or only the very left S for 'segment' b.c.).
If all sites have a form, it respects the form to ensure
that one S is included per bond.
The final state is always in right-canonical 'B' form.

Performs one sweep left to right doing QR decompositions, and one sweep right to left
doing SVDs calculating the singular values.

Parameters
----------
renormalize: bool
Whether a change in the norm should be discarded or used to update :attr:norm.
cutoff : float | None
Cutoff of singular values used in the SVDs.

Returns
-------
U_L, V_R : :class:~tenpy.linalg.np_conserved.Array
Only returned for 'segment' boundary conditions.
The unitaries defining the new left and right Schmidt states in terms of the old ones,
with legs 'vL', 'vR'.
"""
assert (self.finite)
L = self.L
assert (L > 2)  # otherwise implement yourself...
# normalize very left singular values
S = self.get_SL(0)
if self.bc == 'segment':
if S is None:
raise ValueError("Need S[0] and S[L] for segment boundary conditions.")
self.set_SL(0, S / np.linalg.norm(S))
S = self.get_SR(L - 1)
self.set_SR(L - 1, S / np.linalg.norm(S))
else:  # bc == 'finite':
self.set_SL(0, np.array([1.]))  # trivial singular value on very left/right
self.set_SR(L - 1, np.array([1.]))
# sweep from left to right to bring it into left canonical form.
if any([(f is None) for f in self.form]):
# ignore any 'S' and canonical form
M = self.get_B(0, form=None)
form = None
else:
# we actually had a canonical form before, so we should *not* ignore the 'S'
M = self.get_B(0, form='Th')
form = 'B'  # for other 'M'
if self.bc == 'segment':
M.iscale_axis(self.get_SL(0), axis='vL')
Q, R = npc.qr(M.combine_legs(['vL'] + self._p_label), inner_labels=['vR', 'vL'])
# Q = unitary, R has to be multiplied to the right
self.set_B(0, Q.split_legs(0), form='A')
for i in range(1, L - 1):
M = self.get_B(i, form)
M = npc.tensordot(R, M, axes=['vR', 'vL'])
Q, R = npc.qr(M.combine_legs(['vL'] + self._p_label), inner_labels=['vR', 'vL'])
# Q is unitary, i.e. left canonical, R has to be multiplied to the right
self.set_B(i, Q.split_legs(0), form='A')
M = self.get_B(L - 1, None)
M = npc.tensordot(R, M, axes=['vR', 'vL'])
if self.bc == 'segment':
# also neet to calculate new singular values on the very right
U, S, VR_segment = npc.svd(M.combine_legs(['vL'] + self._p_label),
cutoff=cutoff,
inner_labels=['vR', 'vL'])
S /= np.linalg.norm(S)
self.set_SR(L - 1, S)
M = U.scale_axis(S, 1).split_legs(0)
# sweep from right to left, calculating all the singular values
U, S, V = npc.svd(M.combine_legs(['vR'] + self._p_label, qconj=-1),
cutoff=cutoff,
inner_labels=['vR', 'vL'])
if not renormalize:
self.norm = self.norm * np.linalg.norm(S)
S = S / np.linalg.norm(S)  # normalize
self.set_SL(L - 1, S)
self.set_B(L - 1, V.split_legs(1), form='B')
for i in range(L - 2, -1, -1):
M = self.get_B(i, 'A')
M = npc.tensordot(M, U.scale_axis(S, 'vR'), axes=['vR', 'vL'])
U, S, V = npc.svd(M.combine_legs(['vR'] + self._p_label, qconj=-1),
cutoff=cutoff,
inner_labels=['vR', 'vL'])
S = S / np.linalg.norm(S)  # normalize
self.set_SL(i, S)
self.set_B(i, V.split_legs(1), form='B')
if self.bc == 'finite':
assert len(S) == 1
self._B[0] *= U[0, 0]  # just a trivial phase factor, but better keep it
# done. Discard the U for segment bc, although it might be a non-trivial unitary.
# and just re-shuffling of the states left for 'segment' bc)
if self.bc == 'segment':
return U, VR_segment

def canonical_form_infinite(self, renormalize=True, tol_xi=1.e6):
"""Bring an infinite MPS into canonical form (in place).

If any site is in :attr:form None, it does *not* use any of the singular values S.
If all sites have a form, it respects the form to ensure
that one S is included per bond.
The final state is always in right-canonical 'B' form.

Proceeds in three steps, namely 1) diagonalize right and left transfermatrix on a given
bond to bring that bond into canonical form, and then
2) sweep right to left, and 3) left to right to bringing other bonds into canonical form.

.. warning :
You might *loose* precision when calling this function.
When we diagonalize the transfermatrix, we get the singular values squared as
eigenvalues, with numerical noise on the order of machine precision (usually ~1.e-15).
Taking the square root, the new singular values are only precise to *half* the machine
precision (usually ~1.e-7).

Parameters
----------
renormalize: bool
Whether a change in the norm should be discarded or used to update :attr:norm.
tol_xi : float
Raise an error if the correlation length is larger than that
(which indicates a degenerate "cat" state, e.g., for spontaneous symmetry breaking).
"""
assert not self.finite
i1 = np.argmin(self.chi)  # start at this bond
if any([(f is None) for f in self.form]):
# ignore any 'S' and canonical form, just state that we are in 'B' form
self.form = self._parse_form('B')
self._S[i1] = np.ones(self.chi[i1], dtype=np.float)  # (is later used for guess of Gl)
else:
# was in canonical form before; bring back into canonical form
# -> make sure we don't use multiple S on one bond in our definition of the MPS
self.convert_form('B')
L = self.L
Wr_list = [None] * L  # right eigenvectors of TM on each bond after ..._correct_right

# phase 1: bring bond (i1-1, i1) in canonical form
# find dominant right eigenvector
norm, Gr = self._canonical_form_dominant_gram_matrix(i1, False, tol_xi)
self._B[i1] /= np.sqrt(norm)  # correct norm
if not renormalize:
self.norm *= np.sqrt(norm)
# make Gr diagonal to Wr
Wr, Gl = self._canonical_form_correct_right(i1, Gr, return_Gl_guess=True)
# find dominant left eigenvector
norm, Gl = self._canonical_form_dominant_gram_matrix(i1, True, tol_xi, Gl)
if abs(1. - norm) > 1.e-13:
warnings.warn("Although we renormalized the TransferMatrix, "
"the largest eigenvalue is not 1")  # (this shouldn't happen)
self._B[i1] /= np.sqrt(norm)  # correct norm again
if not renormalize:
self.norm *= np.sqrt(norm)
# bring bond to canonical form
Gl, Wr = self._canonical_form_correct_left(i1, Gl, Wr)
# now the bond (i1-1,i1) is in canonical form

# phase 2: sweep from right to left; find other right eigenvectors and make them diagonal
Wr_list[i1] = Wr  # diag(Wr) is right eigenvector on bond (i1-1, i1)
for j1 in range(i1 - 1, i1 - L, -1):
B1 = self.get_B(j1, 'B')
axes = [self._p_label + ['vR'], self._get_p_label('*') + ['vR*']]
Gr = npc.tensordot(B1.scale_axis(Wr, 'vR'), B1.conj(), axes=axes)
Wr_list[j1 % L] = Wr = self._canonical_form_correct_right(j1, Gr)

# phase 3: sweep from left to right; find other left eigenvectors,
# bring each bond into canonical form
for j1 in range(i1 - L + 1, i1, +1):
# find Gl on bond j1-1, j1
B1 = self.get_B(j1 - 1, 'B')
Gl = npc.tensordot(
B1.conj(),  # old B1; now on site j1-1
npc.tensordot(Gl, B1, axes=['vR', 'vL']),
axes=[self._get_p_label('*') + ['vL*'], self._p_label + ['vR*']])
# axes=[['p*', 'vL*'], ['p', 'vR*']])
Gl, Wr = self._canonical_form_correct_left(j1, Gl, Wr_list[j1 % L])

def correlation_length(self, target=1, tol_ev0=1.e-8, charge_sector=0):
r"""Calculate the correlation length by diagonalizing the transfer matrix.

Assumes that self is in canonical form.

Works only for infinite MPS, where the transfer matrix is a useful concept.
Assuming a single-site unit cell, any correlation function splits into
:math:C(A_i, B_j) = A'_i T^{j-i-1} B'_j
with some parts left and right and the :math:j-i-1-th power of the transfer matrix in
between. The largest eigenvalue is 1 (if self is properly normalized)
and gives the dominant contribution of
:math:A'_i E_1 * 1^{j-i-1} * E_1^T B'_j = <A> <B>,
and the second largest one gives a contribution :math:\propto \lambda_2^{j-i-1}.
Thus :math:\lambda_2 = \exp(-\frac{1}{\xi}).

More general for a L-site unit cell we get :math:\lambda_2 = \exp(-\frac{L}{\xi}),
where the xi is given in units of 1 lattice spacing in the MPS.

.. warning ::
For a higher-dimensional lattice (which the MPS class doesn't know about),
the correct unit is the lattice spacing in x-direction, and the correct formula is
:math:\lambda_2 = \exp(-\frac{L_x}{\xi}),
where L_x is the number of lattice spacings in the infinite direction within the
MPS unit cell, e.g. the number of "rings" of a cylinder in the MPS unit cell.
To get to these units, divide the returned xi by the number of sites within a "ring",
for a lattice given in :attr:~tenpy.networks.lattice.N_sites_per_ring.

Parameters
----------
target : int
We look for the target + 1 largest eigenvalues.
tol_ev0 : float
Print warning if largest eigenvalue deviates from 1 by more than tol_ev0.
charge_sector : None | charges | 0
Selects the charge sector in which the dominant eigenvector of the TransferMatrix is.
None stands for *all* sectors, 0 stands for the zero-charge sector.
Defaults to 0, i.e., *assumes* the dominant eigenvector is in charge sector 0.

Returns
-------
xi : float | 1D array
If target=1, return just the correlation length,
otherwise an array of the target largest correlation lengths.
It is measured in units of a single lattice spacing in the MPS language,
see the warning above.
"""
assert (not self.finite)
T = TransferMatrix(self, self, charge_sector=charge_sector, form='B')
num = max(target + 1, self._transfermatrix_keep)
E, _ = T.eigenvectors(num, which='LM')
E = E[np.argsort(-np.abs(E))]  # sort descending by magnitude
if charge_sector is not None and charge_sector != 0:
# need also dominant eigenvector: include 0 charge sector to results
del T
T = TransferMatrix(self, self, charge_sector=0, form='B')
E0, _ = T.eigenvectors(num, which='LM')
assert abs(E0[0]) > abs(E[0]), "dominant eigenvector in zero charge sector?"
E = np.array([E0[0]] + list(E))
if abs(E[0] - 1.) > tol_ev0:
warnings.warn(
"Correlation length: largest eigenvalue not one. "
"Not in canonical form/normalized?",
stacklevel=2)
if len(E) < 2:
return 0.  # only a single eigenvector: zero correlation length
if target == 1:
return -1. / np.log(abs(E[1] / E[0])) * self.L
return -1. / np.log(np.abs(E[1:target + 1] / E[0])) * self.L

def add(self, other, alpha, beta, cutoff=1.e-15):
"""Return an MPS which represents alpha|self> + beta |others>.

Works only for 'finite', 'segment' boundary conditions.
For 'segment' boundary conditions, the virtual legs on the very left/right are
assumed to correspond to each other (i.e. self and other have the same state outside of
the considered segment).
Takes into account :attr:norm.

Parameters
----------
other : :class:MPS
Another MPS of the same length to be added with self.
alpha, beta : complex float
Prefactors for self and other. We calculate
alpha * |self> + beta * |other>
cutoff : float | None
Cutoff of singular values used in the SVDs.

Returns
-------
sum : :class:MPS
An MPS representing alpha|self> + beta |other>.
Has same total charge as self.
U_L, V_R : :class:~tenpy.linalg.np_conserved.Array
Only returned for 'segment' boundary conditions.
The unitaries defining the new left and right Schmidt states in terms of the old ones,
with legs 'vL', 'vR'.
"""
L = self.L
assert (other.L == L and L >= 2)  # (if you need this, generalize this function...)
assert self.finite
self._gauge_compatible_vL_vR(other)
legs = ['vL', 'vR'] + self._p_label
# alpha and beta appear only on the first site
alpha = alpha * self.norm
beta = beta * other.norm
Bs = [
npc.grid_concat(
[[alpha * self.get_B(0).transpose(legs), beta * other.get_B(0).transpose(legs)]],
axes=[0, 1])
]
for i in range(1, L - 1):
B1 = self.get_B(i).transpose(legs)
B2 = other.get_B(i).transpose(legs)
grid = [[B1, npc.zeros([B1.get_leg('vL'), B2.get_leg('vR')] + B1.legs[2:])],
[npc.zeros([B2.get_leg('vL'), B1.get_leg('vR')] + B1.legs[2:]), B2]]
Bs.append(npc.grid_concat(grid, [0, 1]))
Bs.append(
npc.grid_concat(
[[self.get_B(L - 1).transpose(legs)], [other.get_B(L - 1).transpose(legs)]],
axes=[0, 1]))

Ss = [np.ones(1)] + [np.ones(B.shape[1]) for B in Bs]
psi = self.__class__(self.sites, Bs, Ss, 'finite', form=None)  # new class instance
# bring to canonical form, calculate Ss
if self.bc == 'segment':
U_L, V_R = psi.canonical_form_finite(renormalize=False, cutoff=cutoff)
return psi, U_L, V_R
else:
psi.canonical_form_finite(renormalize=False, cutoff=cutoff)
return psi

def apply_local_op(self, i, op, unitary=None, renormalize=False, cutoff=1.e-13):
"""Apply a local (one or multi-site) operator to self.

Note that this destroys the canonical form if the local operator is non-unitary.
Therefore, this function calls :meth:canonical_form if necessary.

Parameters
----------
i : int
(Left-most) index of the site(s) on which the operator should act.
op : str | npc.Array
A physical operator acting on site i, with legs 'p', 'p*' for a single-site
operator or with legs ['p0', 'p1', ...], ['p0*', 'p1*', ...] for an operator
acting on n>=2 sites.
Strings (like 'Id', 'Sz') are translated into single-site operators defined by
:attr:sites.
unitary : None | bool
Whether op is unitary, i.e., whether the canonical form is preserved (True)
or whether we should call :meth:canonical_form (False).
None checks whether norm(op dagger(op) - identity) is smaller than cutoff.
renormalize : bool
Whether the final state should keep track of the norm (False, default) or be
renormalized to have norm 1 (True).
cutoff : float
Cutoff for singular values if op acts on more than one site (see :meth:from_full).
(And used as cutoff for a unspecified unitary.)
"""
i = self._to_valid_index(i)
if isinstance(op, str):
op = self.sites[i].get_op(op)
n = op.rank // 2  # same as int(rank/2)
if n == 1:
pstar, p = 'p*', 'p'
else:
p = self._get_p_labels(n, False)
pstar = self._get_p_labels(n, True)
if unitary is None:
op_op_dagger = npc.tensordot(op, op.conj(), axes=[pstar, p])
if n > 1:
op_op_dagger = op_op_dagger.combine_legs([p, pstar], qconj=[+1, -1])
unitary = npc.norm(op_op_dagger - npc.eye_like(op_op_dagger)) < cutoff
if n == 1:
opB = npc.tensordot(op, self._B[i], axes=['p*', 'p'])
self.set_B(i, opB, self.form[i])
else:
th = self.get_theta(i, n)
th = npc.tensordot(op, th, axes=[pstar, p])
# use MPS.from_full to split the sites
split_th = self.from_full(self.sites[i:i + n], th, None, cutoff, False, 'segment',
(self.get_SL(i), self.get_SR(i + n - 1)))
for j in range(n):
self.set_B(i + j, split_th._B[j], split_th.form[j])
for j in range(n - 1):
self.set_SR(i + j, split_th._S[j + 1])
if not unitary:
self.canonical_form(renormalize)

def swap_sites(self, i, swap_op='auto', trunc_par={}):
"""Swap the two neighboring sites i and i+1 (inplace).

Exchange two neighboring sites: form theta, 'swap' the physical legs and split
with an svd. While the 'swap' is just a transposition/relabeling for bosons, one needs to
be careful about the sign for fermions.

Parameters
----------
i : int
Swap the two sites at positions i and i+1.
swap_op : None | 'auto' | :class:~tenpy.linalg.np_conserved.Array
The operator used to swap the phyiscal legs of the two-site wave function theta.
For None, just transpose/relabel the legs, for 'auto' also take care of
fermionic signs. Alternative give an npc :class:~tenpy.linalg.np_conserved.Array
which represents the full operator used for the swap.
Should have legs ['p0', 'p1', 'p0*', 'p1*'] whith 'p0', 'p1*' contractible.
trunc_par : dict
Parameters for truncation, see :cfg:config:truncation.
Defaults to {'chi_max': max(self.chi)}.

Returns
-------
trunc_err : :class:~tenpy.algorithms.truncation.TruncationError
The error of the represented state introduced by the truncation after the swap.
"""
trunc_par.setdefault('chi_max', max(self.chi))
siteL, siteR = self.sites[self._to_valid_index(i)], self.sites[self._to_valid_index(i + 1)]
if swap_op == 'auto':
# get sign for Fermions.
# If we write the wave function as
# psi = sum_{ [n_i]} psi_[n_i] prod_i (c^dagger_i)^{n_i}  |vac>
# we see that switching i <-> i+1 the phase to be introduced is by commuting
# (c^dagger_i)^{n_i} with (c^dagger_{i+1})^{n_{i+1}}
# This gives a sign (-1)^{n_i * n_{i+1}}.
# site.JW_exponent is the n_i in the above equations, for each physical index.
sign = siteL.JW_exponent[:, np.newaxis] * siteR.JW_exponent[np.newaxis, :]
if np.any(sign):
dL, dR = siteL.dim, siteR.dim
sign = np.real_if_close(np.exp(1.j * np.pi * sign.reshape([dL * dR])))
swap_op = np.diag(sign).reshape([dL, dR, dL, dR])
legs = [siteL.leg, siteR.leg, siteL.leg.conj(), siteR.leg.conj()]
swap_op = npc.Array.from_ndarray(swap_op, legs, labels=['p1', 'p0', 'p0*', 'p1*'])
else:  # no sign necessary
swap_op = None  # continue with transposition as for Bosons
theta = self.get_theta(i, n=2)
if swap_op is None:
# just replace the labels, effectively this is a transposition.
theta.ireplace_labels(['p0', 'p1'], ['p1', 'p0'])
C.ireplace_labels(['p0', 'p1'], ['p1', 'p0'])
elif isinstance(swap_op, npc.Array):
theta = npc.tensordot(swap_op, theta, axes=[['p0*', 'p1*'], ['p0', 'p1']])
C = npc.tensordot(swap_op, C, axes=(['p0*', 'p1*'], ['p0', 'p1']))
else:
raise ValueError("Invalid swap_op: got " + repr(swap_op))
theta = theta.combine_legs([('vL', 'p0'), ('vR', 'p1')], qconj=[+1, -1])
U, S, V, err, renormalize = svd_theta(theta, trunc_par, inner_labels=['vR', 'vL'])
B_R = V.split_legs(1).ireplace_label('p1', 'p')
B_L = npc.tensordot(C.combine_legs(('vR', 'p1'), pipes=theta.legs[1]),
V.conj(),
axes=['(vR.p1)', '(vR*.p1*)'])
B_L.ireplace_labels(['vL*', 'p0'], ['vR', 'p'])
B_L /= renormalize  # re-normalize to <psi|psi> = 1
self.set_SR(i, S)
self.set_B(i, B_L, 'B')
self.set_B(i + 1, B_R, 'B')
self.sites[self._to_valid_index(i)] = siteR  # swap 'sites' as well
self.sites[self._to_valid_index(i + 1)] = siteL
return err

def permute_sites(self, perm, swap_op='auto', trunc_par={}, verbose=0):
"""Applies the permutation perm to the state (inplace).

Parameters
----------
perm : ndarray[ndim=1, int]
The applied permutation, such that psi.permute_sites(perm)[i] = psi[perm[i]]
(where [i] indicates the i-th site).
swap_op : None | 'auto' | :class:~tenpy.linalg.np_conserved.Array
The operator used to swap the phyiscal legs of a two-site wave function theta,
see :meth:swap_sites.
trunc_par : dict
Parameters for truncation, see :cfg:config:truncation.
Defaults to {'chi_max': max(self.chi)}.
verbose : float
Level of verbosity, print status messages if verbose > 0.

Returns
-------
trunc_err : :class:~tenpy.algorithms.truncation.TruncationError
The error of the represented state introduced by the truncation after the swaps.
"""
perm = list(perm)  # gets modified, so we should copy
# In order to keep sites close together, we always scan from the left,
# keeping everything up to i in strictly ascending order.
# => more or less an 'insertion' sort algorithm.
# Works nicely for permutations like [1,2,3,0,6,7,8,5] (swapping the 0 and 5 around).
# For [ 2 3 4 5 6 7 0 1], it splits 0 and 1 apart (first swapping the 0 down, then the 1)
trunc_par.setdefault('chi_max', max(self.chi))
trunc_err = TruncationError()
num_swaps = 0
i = 0
while i < self.L - 1:
if perm[i] > perm[i + 1]:
if verbose > 1:
print(i, ": chi = ", self._S[i + 1].shape[0], end='')
trunc = self.swap_sites(i, swap_op, trunc_par)
if verbose > 1:
print("->", self._S[i + 1].shape[0], ". eps = ", trunc.eps)
num_swaps += 1
x, y = perm[i], perm[i + 1]
perm[i + 1], perm[i] = x, y
# restart from very left; but we know it's already sorted up to i-1
if i > 0:
i -= 1
trunc_err += trunc
else:
i += 1
if verbose > 0:
print("Total swaps in permute_sites:", num_swaps, repr(trunc_err))
return trunc_err

def compute_K(self, perm, swap_op='auto', trunc_par=None, canonicalize=1.e-6, verbose=0):
r"""Compute the momentum quantum numbers of the entanglement spectrum for 2D states.

Works for an infinite MPS living on a cylinder, infinitely long in x direction and with
periodic boundary conditions in y directions.
If the state is invariant under 'rotations' around the cylinder axis, one can find the
momentum quantum numbers of it. (The rotation is nothing more than a translation in y.)
This function permutes some sites (on a copy of self) to enact the rotation, and then
finds the dominant eigenvector of the mixed transfer matrix to get the quantum numbers,
along the lines of [PollmannTurner2012]_, see also (the appendix and Fig. 11 in the arXiv
version of) [CincioVidal2013]_.

Parameters
----------
perm : 1D ndarray | :class:~tenpy.models.lattice.Lattice
Permuation to be applied to the physical indices, see :meth:permute_sites.
If a lattice is given, we use it to read out the lattice structure and shift
each site by one lattice-vector in y-direction (assuming periodic boundary conditions).
(If you have a :class:~tenpy.models.model.CouplingModel,
give its lat attribute for this argument)
swap_op : None | 'auto' | :class:~tenpy.linalg.np_conserved.Array
The operator used to swap the phyiscal legs of a two-site wave function theta,
see :meth:swap_sites.
trunc_par : dict
Parameters for truncation, see :cfg:config:truncation.
If not set, chi_max defaults to max(self.chi).
canonicalize : float
Check that self is in canonical form; call :meth:canonical_form
if :meth:norm_test yields np.linalg.norm(self.norm_test()) > canonicalize.
verbose : float
Level of verbosity, print status messages if verbose > 0.

Returns
-------
U : :class:~tenpy.linalg.np_conserved.Array
Unitary representation of the applied permutation on left Schmidt states.
W : ndarray
1D array of the form S**2 exp(i K), where S are the Schmidt values
on the left bond. You can use :func:np.abs and :func:np.angle to extract the
Schmidt values S and momenta K from W.
q : :class:~tenpy.linalg.charges.LegCharge
LegCharge corresponding to W.
ov : complex
The eigenvalue of the mixed transfer matrix <psi|T|psi> per :attr:L sites.
An absolute value different smaller than 1 indicates that the state is not invariant
under the permutation or that the truncation error trunc_err was too large!
trunc_err : :class:~tenpy.algorithms.truncation.TruncationError
The error of the represented state introduced by the truncation after swaps when
performing the truncation.
"""
from ..models.lattice import Lattice  # dynamical import to avoid import loops
if self.finite:
raise ValueError("Works only for infinite b.c.")
if trunc_par is None:
trunc_par = {}
trunc_par.setdefault('chi_max', max(self.chi))
trunc_par.setdefault('verbose', verbose)

if isinstance(perm, Lattice):
lat = perm
assert lat.dim >= 2  # ensure that the lattice is at least 2D
assert lat.N_sites == self.L
shifted_lat_order = lat.order.copy()
shifted_lat_order[:, 1] = np.mod(shifted_lat_order[:, 1] + 1, lat.Ls[1])
perm = lat.lat2mps_idx(shifted_lat_order)
if verbose > 1:
print("permutation: ", perm)
# preliminary: check canonical form
self.convert_form('B')
norm_err = np.linalg.norm(self.norm_test())
if norm_err > canonicalize:
warnings.warn("self.norm_test() = {0!s} ==> canonicalize".format(self.norm_test()))
self.canonical_form()
# get copy of self
psi_t = self.copy()
# apply permutation
perm = np.asarray(perm)
trunc_err = psi_t.permute_sites(perm, swap_op, trunc_par, verbose / 10.)
# re-check canonical form
norm_err = np.linalg.norm(psi_t.norm_test())
if norm_err > canonicalize:
warnings.warn("psi_t.norm_test() = {0!s} ==> canonicalize".format(psi_t.norm_test()))
psi_t.convert_form('B')
TM = TransferMatrix(self, psi_t, transpose=True, charge_sector=0)
# Find left dominant eigenvector of this mixed transfer matrix.
# Because we are in B form and get the left eigenvector,
# the resulting vector should be sUs up to a scaling.
ov, sUs = TM.eigenvectors(num_ev=self._transfermatrix_keep)
if verbose > 0:
print("compute_K: overlap ", ov[0], ", |o| = 1. -", 1. - np.abs(ov[0]))
# (should be 1 if state is invariant under translations)
print("compute_K: truncation error ", trunc_err.eps)
sUs = sUs[0].split_legs(0)
_, sUs_blocked = sUs.as_completely_blocked()
W = npc.eigvals(sUs_blocked, sort='m>')
# W = s^2 exp(i K ) up to overall scaling
# Strip S's from U
inv_S = 1. / self.get_SL(0)
U = sUs.scale_axis(inv_S, 0).iscale_axis(inv_S, 1)
# U should be unitary - scale it
U *= (np.sqrt(U.shape[0]) / npc.norm(U))
return U, W / np.sum(np.abs(W)), sUs_blocked.legs[0], ov[0], trunc_err

def __str__(self):
"""Some status information about the MPS."""
res = ["MPS, L={L:d}, bc={bc!r}.".format(L=self.L, bc=self.bc)]
res.append("chi: " + str(self.chi))
if self.L > 10:
res.append("first two sites: " + repr(self.sites[0]) + " " + repr(self.sites[1]))
res.append("first two forms:" + " ".join([repr(f) for f in self.form[:2]]))
else:
res.append("sites: " + " ".join([repr(s) for s in self.sites]))
res.append("forms: " + " ".join([repr(f) for f in self.form]))
return "\n".join(res)

def compress(self, options):
"""Compresss an MPS.

Options
-------
.. cfg:config :: MPS_compress
:include: VariationalCompression

compression_method : 'SVD' | 'variational'
Mandatory.
Selects the method to be used for compression.
For the SVD compression, trunc_params is the only other option used.
trunc_params : dict
Truncation parameters as described in :cfg:config:truncation.
"""
options = asConfig(options, "MPS_compress")
method = options['compression_method']
trunc_params = options.subconfig('trunc_params')
if method == 'SVD':
return self.compress_svd(trunc_params)
elif method == 'variational':
from ..algorithms.mps_common import VariationalCompression
return VariationalCompression(self, options).run()
raise ValueError("Unknown compression method: " + repr(method))

def compress_svd(self, trunc_par):
"""Compress self with a single sweep of SVDs; in place.

Perform a single right-sweep of QR/SVD without truncation, followed by a left-sweep with
truncation, very much like :meth:canonical_form_finite.

.. warning ::
In case of a strong compression, this does not find the optimal, global solution.

Parameters
----------
trunc_par : dict
Parameters for truncation, see :cfg:config:truncation.
"""
trunc_err = TruncationError()
if self.bc == 'finite':
# Do QR starting from the left
B = self.get_B(0, form='Th')
for i in range(self.L - 1):
B = B.combine_legs(['vL', 'p'])
q, r = npc.qr(B, inner_labels=['vR', 'vL'])
B = q.split_legs()
self.set_B(i, B, form=None)
B = self.get_B(i + 1, form='B')
B = npc.tensordot(r, B, axes=('vR', 'vL'))
# Do SVD from right to left & truncate
for i in range(self.L - 1, 0, -1):
B = B.combine_legs(['p', 'vR'])
U, S, VH, err, norm_new = svd_theta(B, trunc_par)
trunc_err += err
self.norm *= norm_new
VH = VH.split_legs()
self.set_B(i, VH, form='B')
B = self.get_B(i - 1, form=None)
B = npc.tensordot(B, U, axes=('vR', 'vL'))
B.iscale_axis(S, 'vR')
self.set_SL(i, S)
self.set_B(0, B, form='Th')
elif self.bc == 'infinite':
for i in range(self.L):
theta = self.get_theta(i, n=2)
theta = theta.combine_legs([['vL', 'p0'], ['p1', 'vR']], qconj=[+1, -1])
self.set_svd_theta(i, theta, update_norm=False)
for i in range(self.L - 1, -1, -1):
theta = self.get_theta(i, n=2)
theta = theta.combine_legs([['vL', 'p0'], ['p1', 'vR']], qconj=[+1, -1])
trunc_err += self.set_svd_theta(i, theta, trunc_par, update_norm=False)
else:
raise NotImplementedError("unsupported boundary conditions " + repr(self.bc))
return trunc_err

def _to_valid_index(self, i):
"""Make sure i is a valid index (depending on self.bc)."""
if not self.finite:
return i % self.L
if i < 0:
i += self.L
if i >= self.L or i < 0:
raise ValueError("i = {0:d} out of bounds for finite MPS".format(i))
return i

def _parse_form(self, form):
"""Parse form = (list of) {tuple | key of _valid_forms} to list of tuples"""
if isinstance(form, tuple):
return [form] * self.L
form = to_iterable(form)
if len(form) == 1:
form = [form[0]] * self.L
if len(form) != self.L:
raise ValueError("Wrong len of form: " + repr(form))
return [self._to_valid_form(f) for f in form]

def _to_valid_form(self, form):
"""Parse form = {tuple | key of _valid_forms} to a tuple"""
if isinstance(form, tuple):
return form
return self._valid_forms[form]

def _scale_axis_B(self, B, S, form_diff, axis_B, cutoff):
"""Scale an axis of B with S to bring it in desired form.

If S is just 1D (as usual, e.g. during TEBD), this function just performs
B.scale_axis(S**form_diff, axis_B).

However, during the DMRG with mixer, S might acutally be a 2D matrix.
For form_diff = -1, we need to calculate the inverse of S, more precisely the
(Moore-Penrose) pseudo inverse, see :func:~tenpy.linalg.np_conserved.pinv.
The cutoff is only used in that case.

Returns scaled B.
"""
if form_diff == 0:
return B  # nothing to do
if not isinstance(S, npc.Array):
# the usual case: S is a 1D array with singular values
if form_diff != 1.:
S = S**form_diff
return B.scale_axis(S, axis_B)
else:
# e.g. during DMRG with a DensityMatrixMixer
if S.rank != 2:
raise ValueError("Expect 2D npc.Array or 1D numpy ndarray")
if form_diff == -1:
S = npc.pinv(S, cutoff)
elif form_diff != 1.:
raise ValueError("Can't scale/tensordot a 2D S for non-integer form_diff")

# Hack: mpo.MPOEnvironment.full_contraction uses axis_B == 'vL*'
if axis_B == 'vL' or axis_B == 'vL*':
B = npc.tensordot(S, B, axes=[1, axis_B]).replace_label(0, axis_B)
elif axis_B == 'vR' or axis_B == 'vR*':
B = npc.tensordot(B, S, axes=[axis_B, 0]).replace_label(-1, axis_B)
else:
raise ValueError("This should never happen: unexpected leg for scaling with S")
return B

def _replace_p_label(self, A, s):
"""Return npc Array A with replaced label, 'p' -> 'p'+s.

This is done for each of the 'physical labels' in :attr:_p_label. With a clever use of
this function, the re-implementation of various functions (like get_theta) in derived
classes with multiple legs per site can be avoided.
"""
return A.replace_label('p', 'p' + s)
#  return A.replace_labels(self._p_label, self._get_p_label(s))

def _get_p_label(self, s):
"""return  self._p_label with additional string s."""
return ['p' + s]
#  return [lbl + s for lbl in self._p_label]

def _get_p_labels(self, ks, star=False):
"""Join self._get_p_label(str(k)) for k in range(ks) to a single list."""
if star:
return ['p' + str(k) + '*' for k in range(ks)]
#  return [lbl + str(k) + '*' for k in range(ks) for lbl in self._p_label]
else:
return ['p' + str(k) for k in range(ks)]
#  return [lbl + str(k) for k in range(ks) for lbl in self._p_label]

def _expectation_value_args(self, ops, sites, axes):
"""parse the arguments of self.expectation_value()"""
ops = npc.to_iterable_arrays(ops)
if any(isinstance(op, str) for op in ops):
n = 1
else:
s = 0 if sites is None else to_iterable(sites)[0]
n = ops[s % len(ops)].rank // 2  # same as int(rank/2)
L = self.L
if sites is None:
if self.finite:
sites = range(L - (n - 1))
else:
sites = range(L)
sites = to_iterable(sites)
if axes is None:
if n == 1:
axes = (['p'], ['p*'])
else:
axes = (self._get_p_labels(n), self._get_p_labels(n, True))
# check number of axes
ax_p, ax_pstar = axes
if len(ax_p) != n or len(ax_pstar) != n:
raise ValueError("Len of axes does not match to n-site operator with n=" + str(n))
return ops, sites, n, axes

def _correlation_function_args(self, ops1, ops2, sites1, sites2, opstr):
"""get default arguments of self.correlation_function()"""
if sites1 is None:
sites1 = range(0, self.L)
elif isinstance(sites1, int):
sites1 = range(0, sites1)
if sites2 is None:
sites2 = range(0, self.L)
elif isinstance(sites2, int):
sites2 = range(0, sites2)
ops1 = npc.to_iterable_arrays(ops1)
ops2 = npc.to_iterable_arrays(ops2)
opstr = npc.to_iterable_arrays(opstr)
sites1 = np.sort(sites1)
sites2 = np.sort(sites2)
return ops1, ops2, sites1, sites2, opstr

def _corr_up_diag(self, ops1, ops2, i, j_gtr, opstr, str_on_first, apply_opstr_first):
"""correlation function above the diagonal: for fixed i and all j in j_gtr, j > i."""
op1 = self.get_op(ops1, i)
opstr1 = self.get_op(opstr, i)
if opstr1 is not None and str_on_first:
axes = ['p*', 'p'] if apply_opstr_first else ['p', 'p*']
op1 = npc.tensordot(op1, opstr1, axes=axes)
theta = self.get_theta(i, n=1)
C = npc.tensordot(op1, theta, axes=['p*', 'p0'])
C = npc.tensordot(theta.conj(), C, axes=[['p0*', 'vL*'], ['p', 'vL']])
# C has legs 'vR*', 'vR'
js = list(j_gtr[::-1])  # stack of j, sorted *descending*
res = []
for r in range(i + 1, js[0] + 1):  # js[0] is the maximum
B = self.get_B(r, form='B')
C = npc.tensordot(C, B, axes=['vR', 'vL'])
if r == js[-1]:
Cij = npc.tensordot(self.get_op(ops2, r), C, axes=['p*', 'p'])
Cij = npc.inner(B.conj(), Cij, axes=[['vL*', 'p*', 'vR*'], ['vR*', 'p', 'vR']])
res.append(Cij)
js.pop()
if len(js) > 0:
op = self.get_op(opstr, r)
if op is not None:
C = npc.tensordot(op, C, axes=['p*', 'p'])
C = npc.tensordot(B.conj(), C, axes=[['vL*', 'p*'], ['vR*', 'p']])
return res

def _canonical_form_dominant_gram_matrix(self, bond0, transpose, tol_xi, guess=None):
"""Find dominant eigenvector of the transfer matrix starting between sites (bond0-1,bond0).

Find right (transpose=False) or left (transpose=True) eigenvector of the transfermatrix.
"""
TM = TransferMatrix(self, self, bond0, bond0, transpose=transpose, charge_sector=0)
if guess is None:
diag = self.get_SL(bond0)**2 if transpose else 1.
guess = TM.initial_guess(diag)
guess = guess.combine_legs([0, 1], pipes=TM.pipe)
eta, V = TM.eigenvectors(self._transfermatrix_keep, v0=guess, which='LM')
self._transfermatrix_keep = len(eta)
if len(eta) > 1:
if np.abs(eta[0]) > np.abs(eta[1]):
xi = -self.L / np.log(np.abs((eta[1] / eta[0])))
else:
xi = np.inf
if xi > tol_xi:
raise ValueError("Degenerate spectrum of TransferMatrix "
"(corr length xi={xi:.3e})".format(xi=xi))
eta, G = eta[0], V[0]
G = G.split_legs()
# note: the dominant eigenvector should be hermitian and positive
# removes phase (arbitrary for eigenvectors!) and normalize
# right eigenvectors should have trace chi, left ones trace 1
# (since we expect something close to eye(chi) for right and diag(S**2) for left G)
norm = 1. if transpose else G.shape[0]
G *= norm / npc.trace(G, 0, 1)
if self.dtype.kind != 'c':  # psi is real -> G should be real
eta = np.abs(eta)
G.iunary_blockwise(np.real)
return eta, G  # G has legs vL, vL* or vR, vR*

def _canonical_form_correct_right(self,
i1,
Gr,
return_Gl_guess=False,
eps=2. * np.finfo(np.double).eps):
"""Given the right gram matrix Gr, updated the bond (i0, i1), where i0 = i1 - 1.

Diagonalise Gr = X^H Wr X and update
B[i0] -> B[i0] X^H / norm,
B[i1] -> X B[i1] * norm, where norm = sqrt(chi/sum(Wr)) == sqrt(chi/tr(Gr))
If Wr has (almost) zero entries, reduce the bond dimension at the given bond.
Then Gr -> Wr.
Return Wr normalized to sum(Wr) = chi.
If return_Gl_guess, return also Gl_guess = X S[i1]**2 X^H.
"""
Gr.itranspose(['vL', 'vL*'])
W, XH = npc.eigh(Gr)  # -> XH has legs vL vL* = vL vR
if np.sign(W[np.argmax(np.abs(W))]) == -1:  # fix sign
W = -W  # should actually never happen:  we initially normalize tr(Gr) = chi > 0
# discard small values on order of machine precision
proj = (W > eps)
if np.count_nonzero(proj) < len(W):
# project into non-degenerate subspace, reducing the bond dimensions!
warnings.warn("canonical_form_infinite: project to smaller bond dimension",
stacklevel=3)
XH.iproject(proj, axes=1)
W = W[proj]
norm = len(W) / np.sum(W)
W *= norm
norm = np.sqrt(norm)  # (norm doesn't change eigenvalue of TM)
Kl = XH.iset_leg_labels(['vL', 'vR']) * (1. / norm)
Kr = XH.transpose().iconj().iset_leg_labels(['vL', 'vR']) * norm
i0 = i1 - 1
self.set_B(i0, npc.tensordot(self.get_B(i0), Kl, axes=['vR', 'vL']))
self.set_B(i1, npc.tensordot(Kr, self.get_B(i1), axes=['vR', 'vL']))
if return_Gl_guess:
guess = npc.tensordot(Kr.scale_axis(self.get_SL(i1)**2, 1), Kl, axes=['vR', 'vL'])
return W, guess.iset_leg_labels(['vR*', 'vR'])
return W

def _canonical_form_correct_left(self, i1, Gl, Wr, eps=2. * np.finfo(np.double).eps):
"""Bring into canonical form on bond (i0, i1) where i0= i1 - 1.

Given the left Gram matrix Gl (with legs 'vR*', 'vR')
and right diag(Wr), compute and diagonalize the density matrix
rho = sqrt(Wr) Gl sqrt(Wr) -> Y^H S^2 Y (Y acting on ket, Y^H on bra).
Then we can update
B[i0] -> B[i0] sqrt(Wr) Y^H
B[i1] -> Y 1/sqrt(Wr) B[i1]
Thus the new dominant left eigenvector is
Gl -> Y sqrt(Wr) Gl sqrt(Wr) Y^H = S^2
and dominant right eigenvector is
diag(Wr) -> Y 1/sqrt(Wr) diag(Wr) 1/sqrt(W) Y^H = diag(1)
i.e., we brought the bond to canonical form and S is the Schmidt spectrum.
"""
sqrt_Wr = np.sqrt(Wr)
Gl.itranspose(['vR*', 'vR'])
rhor = Gl.scale_axis(sqrt_Wr, 0).iscale_axis(sqrt_Wr, 1)
S2, YH = npc.eigh(rhor, sort='>')  # YH has legs 'vR*', 'vR'
S2 /= np.sum(S2)  # equivalent to normalizing tr(rhor)=1
s_norm = 1.
# discard small values on order of machine precision
proj = (S2 > eps)
if np.count_nonzero(proj) < len(S2):
# project into non-degenerate subspace, reducing the bond dimensions!
warnings.warn("canonical_form_infinite: project to smaller bond dimension",
stacklevel=2)
YH.iproject(proj, axes=1)
S2 = S2[proj]
s_norm = np.sqrt(np.sum(S2))
S = np.sqrt(S2) / s_norm
self.set_SL(i1, S)
Yl = YH.scale_axis(sqrt_Wr / s_norm, 0).iset_leg_labels(['vL', 'vR'])
Yr = YH.transpose().iconj().scale_axis(1. / sqrt_Wr, 1).iset_leg_labels(['vL', 'vR'])
i0 = i1 - 1
self.set_B(i0, npc.tensordot(self.get_B(i0), Yl, axes=['vR', 'vL']))
self.set_B(i1, npc.tensordot(Yr, self.get_B(i1), axes=['vR', 'vL']))
Gl = npc.tensordot(Gl, Yl, axes=['vR', 'vL'])
Gl = npc.tensordot(Yl.conj(), Gl, axes=['vL*', 'vR*'])  # labels 'vR*', 'vR'
Gl /= npc.trace(Gl)
# Gl is diag(S**2) up to numerical errors...
return Gl, np.ones(Yr.legs[0].ind_len, np.float)

def _gauge_compatible_vL_vR(self, other):
"""If necessary, gauge total charge of other to match the vL, vR legs of self."""
if self.chinfo.qnumber == 0:
return
from tenpy.tools import optimization
need_gauge = False
with optimization.temporary_level(optimization.OptimizationFlag.default):
try:
other._B[0].get_leg('vL').test_equal(self._B[0].get_leg('vL'))
other._B[-1].get_leg('vR').test_equal(self._B[-1].get_leg('vR'))
except ValueError:
need_gauge = True
if need_gauge:
other.gauge_total_charge(None, self._B[0].get_leg('vL'), self._B[-1].get_leg('vR'))
if any(self.get_total_charge() != other.get_total_charge()):
raise ValueError("self and other have different total charges!")

class MPSEnvironment:
"""Stores partial contractions of :math:<bra|Op|ket> for local operators Op.

The network for a contraction :math:<bra|Op|ket> of a local operator Op, say exemplary
at sites i, i+1 looks like::

|     .-----M[0]--- ... --M[1]---M[2]--- ... ->--.
|     |     |             |      |               |
|     |     |             |------|               |
|     LP[0] |             |  Op  |               RP[-1]
|     |     |             |------|               |
|     |     |             |      |               |
|     .-----N[0]*-- ... --N[1]*--N[2]*-- ... -<--.

Of course, we can also calculate the overlap <bra|ket> by using the special case Op = Id.

We use the following label convention (where arrows indicate qconj)::

|    .-->- vR           vL ->-.
|    |                        |
|    LP                       RP
|    |                        |
|    .--<- vR*         vL* -<-.

To avoid recalculations of the whole network e.g. in the DMRG sweeps,
we store the contractions up to some site index in this class.
For bc='finite','segment', the very left and right part LP[0] and
RP[-1] are trivial and don't change,
but for bc='infinite' they are might be updated
(by inserting another unit cell to the left/right).

The MPS bra and ket have to be in canonical form.
All the environments are constructed without the singular values on the open bond.
In other words, we contract left-canonical A to the left parts LP
and right-canonical B to the right parts RP.
Thus, the special case ket=bra should yield identity matrices for LP and RP.

Parameters
----------
bra : :class:~tenpy.networks.mps.MPS
The MPS to project on. Should be given in usual 'ket' form;
we call conj() on the matrices directly.
Stored in place, without making copies.
If necessary to match charges, we call :meth:~tenpy.networks.mps.MPS.gauge_total_charge.
ket : :class:~tenpy.networks.mpo.MPO | None
The MPS on which the local operator acts.
Stored in place, without making copies.
If None, use bra.
init_LP : None | :class:~tenpy.linalg.np_conserved.Array
Initial very left part LP. If None, build trivial one with :meth:init_LP.
init_RP : None | :class:~tenpy.linalg.np_conserved.Array
Initial very right part RP. If None, build trivial one with :meth:init_RP.
age_LP : int
The number of physical sites involved into the contraction yielding firstLP.
age_RP : int
The number of physical sites involved into the contraction yielding lastRP.

Attributes
----------
L : int
Number of physical sites involved into the Environment, i.e. the least common multiple
of bra.L and ket.L.
bra, ket : :class:~tenpy.networks.mps.MPS
The two MPS for the contraction.
dtype : type
The data type.
_finite : bool
Whether the boundary conditions of the MPS are finite.
_LP : list of {None | :class:~tenpy.linalg.np_conserved.Array}
Left parts of the environment, len L.
LP[i] contains the contraction strictly left of site i
(or None, if we don't have it calculated).
_RP : list of {None | :class:~tenpy.linalg.np_conserved.Array}
Right parts of the environment, len L.
RP[i] contains the contraction strictly right of site i
(or None, if we don't have it calculated).
_LP_age : list of int | None
Used for book-keeping, how large the DMRG system grew:
_LP_age[i] stores the number of physical sites invovled into the contraction
network which yields self._LP[i].
_RP_age : list of int | None
Used for book-keeping, how large the DMRG system grew:
_RP_age[i] stores the number of physical sites invovled into the contraction
network which yields self._RP[i].
"""
def __init__(self, bra, ket, init_LP=None, init_RP=None, age_LP=0, age_RP=0):
if ket is None:
ket = bra
if ket is not bra:
ket._gauge_compatible_vL_vR(bra)  # ensure matching charges
self.bra = bra
self.ket = ket
self.dtype = np.find_common_type([bra.dtype, ket.dtype], [])
self.L = L = lcm(bra.L, ket.L)
self._finite = bra.finite
self._LP = [None] * L
self._RP = [None] * L
self._LP_age = [None] * L
self._RP_age = [None] * L
self._finite = self.ket.finite  # just for _to_valid_index
if init_LP is None:
init_LP = self.init_LP(0)
self.set_LP(0, init_LP, age=age_LP)
if init_RP is None:
init_RP = self.init_RP(L - 1)
self.set_RP(L - 1, init_RP, age=age_RP)
self.test_sanity()

def test_sanity(self):
"""Sanity check, raises ValueErrors, if something is wrong."""
assert (self.bra.finite == self.ket.finite == self._finite)
# check that the network is contractible
for i in range(self.L):
b_s = self.bra.sites[i % self.bra.L]
k_s = self.ket.sites[i % self.ket.L]
b_s.leg.test_equal(k_s.leg)
assert any([LP is not None for LP in self._LP])
assert any([RP is not None for RP in self._RP])

def init_LP(self, i):
"""Build initial left part LP.

Parameters
----------
i : int
Build LP left of site i.

Returns
-------
init_LP : :class:~tenpy.linalg.np_conserved.Array
Identity contractible with the vL leg of ket.get_B(i), labels 'vR*', 'vR'.
"""
leg_ket = self.ket.get_B(i, None).get_leg('vL')
leg_bra = self.bra.get_B(i, None).get_leg('vL')
leg_ket.test_equal(leg_bra)
init_LP = npc.diag(1., leg_ket, dtype=self.dtype, labels=['vR*', 'vR'])
return init_LP

def init_RP(self, i):
"""Build initial right part RP for an MPS/MPOEnvironment.

Parameters
----------
i : int
Build RP right of site i.

Returns
-------
init_RP : :class:~tenpy.linalg.np_conserved.Array
Identity contractible with the vR leg of ket.get_B(i), labels 'vL*', 'vL'.
"""
leg_ket = self.ket.get_B(i, None).get_leg('vR')
leg_bra = self.bra.get_B(i, None).get_leg('vR')
leg_ket.test_equal(leg_bra)
init_RP = npc.diag(1., leg_ket, dtype=self.dtype, labels=['vL*', 'vL'])
return init_RP

def get_LP(self, i, store=True):
"""Calculate LP at given site from nearest available one (including i).

The returned LP_i corresponds to the following contraction,
where the M's and the N's are in the 'A' form::

|     .-------M[0]--- ... --M[i-1]--->-   'vR'
|     |       |             |
|     LP[0]   |             |
|     |       |             |
|     .-------N[0]*-- ... --N[i-1]*--<-   'vR*'

Parameters
----------
i : int
The returned LP will contain the contraction *strictly* left of site i.
store : bool
Wheter to store the calculated LP in self (True) or discard them (False).

Returns
-------
LP_i : :class:~tenpy.linalg.np_conserved.Array
Contraction of everything left of site i,
with labels 'vR*', 'vR' for bra, ket.
"""
# find nearest available LP to the left.
for i0 in range(i, i - self.L, -1):
LP = self._LP[self._to_valid_index(i0)]
if LP is not None:
break
# (for finite, LP[0] should always be set, so we should abort at latest with i0=0)
else:  # no break called
raise ValueError("No left part in the system???")
age_i0 = self.get_LP_age(i0)
for j in range(i0, i):
LP = self._contract_LP(j, LP)
if store:
self.set_LP(j + 1, LP, age=age_i0 + j - i0 + 1)
return LP

def get_RP(self, i, store=True):
"""Calculate RP at given site from nearest available one (including i).

The returned RP_i corresponds to the following contraction,
where the M's and the N's are in the 'B' form::

|     'vL'  ->---M[i+1]-- ... --M[L-1]----.
|                |              |         |
|                |              |         RP[-1]
|                |              |         |
|     'vL*' -<---N[i+1]*- ... --N[L-1]*---.

Parameters
----------
i : int
The returned RP will contain the contraction *strictly* right of site i.
store : bool
Wheter to store the calculated RP in self (True) or discard them (False).

Returns
-------
RP_i : :class:~tenpy.linalg.np_conserved.Array
Contraction of everything left of site i,
with labels 'vL*', 'vL' for bra, ket.
"""
# find nearest available RP to the right.
for i0 in range(i, i + self.L):
RP = self._RP[self._to_valid_index(i0)]
if RP is not None:
break
# (for finite, RP[-1] should always be set, so we should abort at latest with i0=L-1)
age_i0 = self.get_RP_age(i0)
for j in range(i0, i, -1):
RP = self._contract_RP(j, RP)
if store:
self.set_RP(j - 1, RP, age=age_i0 + i0 - j + 1)
return RP

def get_LP_age(self, i):
"""Return number of physical sites in the contractions of get_LP(i).

Might be None.
"""
return self._LP_age[self._to_valid_index(i)]

def get_RP_age(self, i):
"""Return number of physical sites in the contractions of get_RP(i).

Might be None.
"""
return self._RP_age[self._to_valid_index(i)]

def set_LP(self, i, LP, age):
"""Store part to the left of site i."""
i = self._to_valid_index(i)
self._LP[i] = LP
self._LP_age[i] = age

def set_RP(self, i, RP, age):
"""Store part to the right of site i."""
i = self._to_valid_index(i)
self._RP[i] = RP
self._RP_age[i] = age

def del_LP(self, i):
"""Delete stored part strictly to the left of site i."""
i = self._to_valid_index(i)
self._LP[i] = None
self._LP_age[i] = None

def del_RP(self, i):
"""Delete storde part scrictly to the right of site i."""
i = self._to_valid_index(i)
self._RP[i] = None
self._RP_age[i] = None

def get_initialization_data(self):
"""Return data for (re-)initialization.

The returned parameters are collected in a dictionary with the following names.

Returns
-------
init_LP, init_RP : :class:~tenpy.linalg.np_conserved.Array
LP on the left of site 0 and RP on the right of site L-1, which can be
used as init_LP and init_RP for the initialization of a new environment.
age_LP, age_RP : int
The number of physical sites involved into the contraction yielding init_LP and
init_RP, respectively.
"""
L = self.ket.L
data = {'init_LP': self.get_LP(0, True), 'init_RP': self.get_RP(L - 1, True)}
data['age_LP'] = self.get_LP_age(0)
data['age_RP'] = self.get_RP_age(L - 1)
return data

def full_contraction(self, i0):
"""Calculate the overlap by a full contraction of the network.

The full contraction of the environments gives the overlap <bra|ket>,
taking into account :attr:MPS.norm of both bra and ket.
For this purpose, this function contracts
get_LP(i0+1, store=False) and get_RP(i0, store=False) with appropriate singular
values in between.

Parameters
----------
i0 : int
Site index.
"""
if self.ket.finite and i0 + 1 == self.L:
# special case to handle _to_valid_index correctly:
# get_LP(L) is not valid for finite b.c, so we use need to calculate it explicitly.
LP = self.get_LP(i0, store=False)
LP = self._contract_LP(i0, LP)
else:
LP = self.get_LP(i0 + 1, store=False)
# multiply with S: a bit of a hack: use 'private' MPS._scale_axis_B
S_bra = self.bra.get_SR(i0).conj()
LP = self.bra._scale_axis_B(LP, S_bra, form_diff=1., axis_B='vR*', cutoff=0.)
# cutoff is not used for form_diff = 1
S_ket = self.ket.get_SR(i0)
LP = self.bra._scale_axis_B(LP, S_ket, form_diff=1., axis_B='vR', cutoff=0.)
RP = self.get_RP(i0, store=False)
contr = npc.inner(LP, RP, axes=[['vR*', 'vR'], ['vL*', 'vL']], do_conj=False)
return contr * self.bra.norm * self.ket.norm

def expectation_value(self, ops, sites=None, axes=None):
"""Expectation value <bra|ops|ket> of (n-site) operator(s).

Calculates n-site expectation values of operators sandwiched between bra and ket.
For examples the contraction for a two-site operator on site i would look like::

|          .--S--B[i]--B[i+1]--.
|          |     |     |       |
|          |     |-----|       |
|          LP[i] | op  |       RP[i+1]
|          |     |-----|       |
|          |     |     |       |
|          .--S--B*[i]-B*[i+1]-.

Here, the B are taken from ket, the B* from bra.
The call structure is the same as for :meth:MPS.expectation_value.

.. warning :
In contrast to :meth:MPS.expectation_value, this funciton does not normalize,
thus it also takes into account :attr:MPS.norm of both bra and ket.

Parameters
----------
ops : (list of) { :class:~tenpy.linalg.np_conserved.Array | str }
The operators, for wich the expectation value should be taken,
All operators should all have the same number of legs (namely 2 n).
If less than len(sites) operators are given, we repeat them periodically.
Strings (like 'Id', 'Sz') are translated into single-site operators defined by
:attr:sites.
sites : list
List of site indices. Expectation values are evaluated there.
If None (default), the entire chain is taken (clipping for finite b.c.)
axes : None | (list of str, list of str)
Two lists of each n leg labels giving the physical legs of the operator used for
contraction. The first n legs are contracted with conjugated B,
the second n legs with the non-conjugated B.
None defaults to (['p'], ['p*']) for single site (n=1), or
(['p0', 'p1', ... 'p{n-1}'], ['p0*', 'p1*', .... 'p{n-1}*']) for n > 1.

Returns
-------
exp_vals : 1D ndarray
Expectation values, exp_vals[i] = <bra|ops[i]|ket>, where ops[i] acts on
site(s) j, j+1, ..., j+{n-1} with j=sites[i].

Examples
--------
One site examples (n=1):

>>> env.expectation_value('Sz')
[Sz0, Sz1, ..., Sz{L-1}]
>>> env.expectation_value(['Sz', 'Sx'])
[Sz0, Sx1, Sz2, Sx3, ... ]
>>> env.expectation_value('Sz', sites=[0, 3, 4])
[Sz0, Sz3, Sz4]

Two site example (n=2), assuming homogeneous sites:

>>> SzSx = npc.outer(psi.sites[0].Sz.replace_labels(['p', 'p*'], ['p0', 'p0*']),
psi.sites[1].Sx.replace_labels(['p', 'p*'], ['p1', 'p1*']))
>>> env.expectation_value(SzSx)
[Sz0Sx1, Sz1Sx2, Sz2Sx3, ... ]   # with len L-1 for finite bc, or L for infinite

Example measuring <bra|SzSx|ket> on each second site, for inhomogeneous sites:

>>> SzSx_list = [npc.outer(psi.sites[i].Sz.replace_labels(['p', 'p*'], ['p0', 'p0*']),
...                        psi.sites[i+1].Sx.replace_labels(['p', 'p*'], ['p1', 'p1*']))
...              for i in range(0, psi.L-1, 2)]
>>> env.expectation_value(SzSx_list, range(0, psi.L-1, 2))
[Sz0Sx1, Sz2Sx3, Sz4Sx5, ...]
"""
ops, sites, n, (op_ax_p, op_ax_pstar) = self.ket._expectation_value_args(ops, sites, axes)
ax_p = ['p' + str(k) for k in range(n)]
ax_pstar = ['p' + str(k) + '*' for k in range(n)]
E = []
for i in sites:
LP = self.get_LP(i, store=True)
RP = self.get_RP(i + n - 1, store=True)
op = self.ket.get_op(ops, i)
op = op.replace_labels(op_ax_p + op_ax_pstar, ax_p + ax_pstar)
C = self.ket.get_theta(i, n)
C = npc.tensordot(op, C, axes=[ax_pstar, ax_p])  # same labels
C = npc.tensordot(LP, C, axes=['vR', 'vL'])  # axes_p + (vR*, vR)
C = npc.tensordot(C, RP, axes=['vR', 'vL'])  # axes_p + (vR*, vL*)
C.ireplace_labels(['vR*', 'vL*'], ['vL', 'vR'])  # back to original theta labels
theta_bra = self.bra.get_theta(i, n)
E.append(npc.inner(theta_bra, C, axes='labels', do_conj=True))
return np.real_if_close(np.array(E)) * self.bra.norm * self.ket.norm

def _contract_LP(self, i, LP):
"""Contract LP with the tensors on site i to form self._LP[i+1]"""
LP = npc.tensordot(LP, self.ket.get_B(i, form='A'), axes=('vR', 'vL'))
axes = (self.ket._get_p_label('*') + ['vL*'], self.ket._p_label + ['vR*'])
# for a ususal MPS, axes = (['p*', 'vL*'], ['p', 'vR*'])
LP = npc.tensordot(self.bra.get_B(i, form='A').conj(), LP, axes=axes)
return LP  # labels 'vR*', 'vR'

def _contract_RP(self, i, RP):
"""Contract RP with the tensors on site i to form self._RP[i-1]"""
RP = npc.tensordot(self.ket.get_B(i, form='B'), RP, axes=('vR', 'vL'))
axes = (self.ket._p_label + ['vL*'], self.ket._get_p_label('*') + ['vR*'])
# for a ususal MPS, axes = (['p', 'vL*'], ['p*', 'vR*'])
RP = npc.tensordot(RP, self.bra.get_B(i, form='B').conj(), axes=axes)
return RP  # labels 'vL', 'vL*'

def _to_valid_index(self, i):
"""Make sure i is a valid index (depending on finite)."""
if not self._finite:
return i % self.L
if i < 0:
i += self.L
if i >= self.L or i < 0:
raise ValueError("i = {0:d} out of bounds for MPSEnvironment".format(i))
return i

class TransferMatrix(sparse.NpcLinearOperator):
r"""Transfer matrix of two MPS (bra & ket).

For an iMPS in the thermodynamic limit, we often need to find the 'dominant RP' (and LP).
This mean nothing else than to take the transfer matrix of the unit cell and find the
(right/left) eigenvector with the largest (magnitude) eigenvalue, since it will dominate
:math:(TM)^n RP (or :math:LP (TM)^n) in the limit :math:n \rightarrow \infty - whatever
the initial RP is. This class provides exactly that functionality with :meth:eigenvectors.

Given two MPS, we define the transfer matrix as::

|    ---M[i]---M[i+1]- ... --M[i+L]---
|       |      |             |
|    ---N[j]*--N[j+1]* ... --N[j+L]*--

Here the M denotes the matrices of the bra and N the ones of the ket, respectively.
To view it as a matrix, we combine the left and right indices to pipes::

|  (vL.vL*) ->-TM->- (vR.vR*)   acting on  (vL.vL*) ->-RP

Note that we keep all M and N as copies.

.. deprecated :: 0.6.0
The default for shift_ket was the value of shift_bra, this will be changed to 0.

Parameters
----------
bra : MPS
The MPS which is to be (complex) conjugated.
ket : MPS
The MPS which is not (complex) conjugated.
shift_bra : int
We start the N of the bra at site shift_bra (i.e. the j in the above network).
shift_ket : int | None
We start the M of the ket at site shift_ket (i.e. the i in the above network).
None is deprecated, default will be changed to 0 in the future.
transpose : bool
Wheter self.matvec acts on RP (False) or LP (True).
charge_sector : None | charges | 0
Selects the charge sector of the vector onto which the Linear operator acts.
None stands for *all* sectors, 0 stands for the zero-charge sector.
Defaults to 0, i.e., **assumes** the dominant eigenvector is in charge sector 0.
form : 'B' | 'A' | 'C' | 'G' | 'Th' | None | tuple(float, float)
In which canonical form we take the M and N matrices.

Attributes
----------
L : int
Number of physical sites involved in the transfer matrix, i.e. the least common multiple
of bra.L and ket.L.
shift_bra : int
We start the N of the bra at site shift_bra.
shift_ket : int | None
We start the M of the ket at site shift_ket. None defaults to shift_bra.
transpose : bool
Wheter self.matvec acts on RP (True) or LP (False).
qtotal : charges
Total charge of the transfer matrix (which is gauged away in matvec).
form : tuple(float, float) | None
In which canonical form (all of) the M and N matrices are.
flat_linop : :class:~tenpy.linalg.sparse.FlatLinearOperator
Class lifting :meth:matvec to ndarrays in order to use :func:~tenpy.tools.math.speigs.
pipe : :class:~tenpy.linalg.charges.LegPipe
Pipe corresponding to '(vL.vL*)' for transpose=False
or to '(vR.vR*)' for transpose=True.
label_split :
['vL', 'vL*'] if tranpose=False or ['vR', 'vR*'] if transpose=True.
_bra_N : list of npc.Array
Complex conjugated matrices of the bra, transposed for fast matvec.
_ket_M : list of npc.Array
The matrices of the ket, transposed for fast matvec.
_contract_legs : int
Number of physical legs per site + 1.
"""
def __init__(self,
bra,
ket,
shift_bra=0,
shift_ket=None,
transpose=False,
charge_sector=0,
form='B'):
L = self.L = lcm(bra.L, ket.L)
if shift_ket is None:
if shift_bra != 0:
warnings.warn("default for shift_ket will change to 0. Specify both explicitly!",
FutureWarning, 2)
shift_ket = shift_bra
self.shift_bra = shift_bra
self.shift_ket = shift_ket
self.transpose = transpose
if ket.chinfo != bra.chinfo:
raise ValueError("incompatible charges")
form = ket._to_valid_form(form)
p = ket._p_label  # for ususal MPS just ['p']
assert p == bra._p_label
pstar = ket._get_p_label('*')  # ['p*']
if not transpose:  # right to left
label = '(vL.vL*)'  # what we act on
label_split = ['vL', 'vL*']
M = self._ket_M = [
ket.get_B(i, form=form).itranspose(['vL'] + p + ['vR'])
for i in reversed(range(shift_ket, shift_ket + L))
]
N = self._bra_N = [
bra.get_B(i, form=form).conj().itranspose(pstar + ['vR*', 'vL*'])
for i in reversed(range(shift_bra, shift_bra + L))
]
pipe = npc.LegPipe([M[0].get_leg('vR'), N[0].get_leg('vR*')], qconj=-1).conj()
else:  # left to right
label = '(vR*.vR)'  # mathematically more natural
label_split = ['vR*', 'vR']
M = self._ket_M = [
ket.get_B(i, form=form).itranspose(['vL'] + p + ['vR'])
for i in range(shift_ket, shift_ket + L)
]
N = self._bra_N = [
bra.get_B(i, form=form).conj().itranspose(['vR*', 'vL*'] + pstar)
for i in range(shift_bra, shift_bra + L)
]
pipe = npc.LegPipe([N[0].get_leg('vL*'), M[0].get_leg('vL')], qconj=+1).conj()
dtype = np.promote_types(bra.dtype, ket.dtype)
self.pipe = pipe
self.label_split = label_split
self.flat_linop = sparse.FlatLinearOperator(self.matvec, pipe, dtype, charge_sector, label)
self.qtotal = bra.chinfo.make_valid(np.sum([B.qtotal for B in M + N], axis=0))
self._contract_legs = len(ket._p_label) + 1  # for a ususal MPS: 2

def matvec(self, vec):
"""Given vec as an npc.Array, apply the transfer matrix.

Parameters
----------
vec : :class:~tenpy.linalg.np_conserved.Array
Vector to act on with the transfermatrix.
If not transposed, vec is the right part RP of an environment,
with legs '(vL.vL*)' in a pipe or splitted.
If transposed, the left part LP of an environment with legs '(vR*.vR)'.

Returns
-------
mat_vec : :class:~tenpy.linalg.np_conserved.Array
The tranfer matrix acted on vec, in the same form as given.
"""
pipe = None
if vec.rank == 1:
vec = vec.split_legs(0)
pipe = self.pipe
vec.itranspose(self.label_split)  # ['vL', 'vL*'] or ['vR*', 'vR']
qtotal = vec.qtotal
legs = vec.legs
contract = self._contract_legs  # number of physical legs per site + 1
# the actual work
if not self.transpose:  # right to left
for N, M in zip(self._bra_N, self._ket_M):
vec = npc.tensordot(M, vec, axes=1)  # axes=['vR', 'vL']
vec = npc.tensordot(vec, N, axes=contract)  # [['p', 'vL*'], ['p*', 'vR*']]
else:  # left to right
for N, M in zip(self._bra_N, self._ket_M):
vec = npc.tensordot(vec, M, axes=1)  # axes=['vR', 'vL']
vec = npc.tensordot(N, vec, axes=contract)  # [['vL*', 'p*'], ['vR*', 'p']])
if np.any(self.qtotal != 0):
# Hack: replace leg charges and qtotal -> effectively gauge self.qtotal away.
vec.qtotal = qtotal
vec.legs = legs
vec.test_sanity()  # Should be fine, but who knows...
if pipe is not None:
vec = vec.combine_legs([0, 1], pipes=pipe)
return vec

def initial_guess(self, diag=1.):
"""Return a diagonal matrix as initial guess for the eigenvector.

Parameters
----------
diag : float | 1D ndarray
Should be 1. for the identity or some singular values squared.

Returns
-------
mat : :class:~tenpy.linalg.np_conserved.Array
A 2D array with diag on the diagonal such that :meth:matvec can act on it.
"""
return npc.diag(diag, self.pipe.legs[0], labels=self.label_split)

def eigenvectors(self,
num_ev=1,
max_num_ev=None,
max_tol=1.e-12,
which='LM',
v0=None,
**kwargs):
"""Find (dominant) eigenvector(s) of self using :mod:scipy.sparse.

If no charge_sector was selected, we look in *all* charge sectors.

Parameters
----------
num_ev : int
Number of eigenvalues/vectors to look for.
max_num_ev : int
:func:scipy.sparse.linalg.speigs somtimes raises a NoConvergenceError for small
num_ev, which might be avoided by increasing num_ev. As a work-around,
we try it again in the case of an error, just with larger num_ev up to max_num_ev.
None defaults to num_ev + 2.
max_tol : float
After the first NoConvergenceError we increase the tol argument to that value.
which : str
Which eigenvalues to look for, see scipy.sparse.linalg.speigs.
**kwargs :
Further keyword arguments given to :func:~tenpy.tools.math.speigs.

Returns
-------
eta : 1D ndarray
The eigenvalues, sorted according to which.
w : list of :class:~tenpy.linalg.np_conserved.Array
The eigenvectors corresponding to eta, as npc.Array with LegPipe.
"""
if max_num_ev is None:
max_num_ev = num_ev + 2
flat_linop = self.flat_linop
if flat_linop.charge_sector is None:
# Try for all charge sectors
eta = []
A = []
for chsect in flat_linop.possible_charge_sectors:
flat_linop.charge_sector = chsect
eta_cs, A_cs = self.eigenvectors(num_ev, max_num_ev, max_tol, which, **kwargs)
eta.extend(eta_cs)
A.extend(A_cs)
flat_linop.charge_sector = None
else:
if v0 is not None:
kwargs['v0'] = self.flat_linop.npc_to_flat(v0)
# for given charge sector
for k in range(num_ev, max_num_ev + 1):
if k > num_ev:
warnings.warn("TransferMatrix: increased num_ev to " + str(k + 1))
try:
eta, A = speigs(flat_linop, k=k, which='LM', **kwargs)
A = np.real_if_close(A)
A = [flat_linop.flat_to_npc(A[:, j]) for j in range(A.shape[1])]
break
except scipy.sparse.linalg.eigen.arpack.ArpackNoConvergence:
if k == max_num_ev:
raise
# just retry with larger k and 'tol'
kwargs['tol'] = max(max_tol, kwargs.get('tol', 0))
# sort
perm = argsort(eta, which)
return np.array(eta)[perm], [A[j] for j in perm]

def build_initial_state(size, states, filling, mode='random', seed=None):
"""Build an "initial state" list.

Uses two iterables ('states' and 'filling') to determine how to fill the
state. The two lists should have the same length as every element in 'filling' gives the filling
fraction for the corresponding state in 'states'.

Example:
size = 6, states = [0, 1, 2], filling = [1./3, 2./3, 0.]
n_states = size * filling = [2, 4, 0]
==> Two sites will get state 0, 4 sites will get state 1, 0 sites will
get state 2.

.. todo ::
Make more general: it should be possible to specify states as strings.

Parameters
----------
size : int
length of state
states : iterable
Containing the possible local states
filling : iterable
Fraction of the total number of sites to get a certain state. If
infinite fractions (e.g. 1/3) are needed, one should supply a
fraction (1./3.)
mode : str | None
State filling pattern. Only 'random' is implemented
seed : int | None
Seed for random number generators

Returns
-------
initial_state (list) : the initial state

Raises
------
ValueError
If fractonal fillings are incommensurate with system size.
AssertionError
If the total filling is not equal to 1, or the length of filling
does not equal the length of states.
"""

random.seed(seed)

# Do some safety checks
assert sum(filling) == 1
assert len(states) == len(filling)

# Get number of sites for each local state
n_states = np.array(filling) * size
for num in n_states:
if ((num - round(num)) < 1e-12):
num = int(round(num))
if type(num) != int and not num.is_integer():
raise ValueError("Cannot create model of length {} with filling {}".format(
size, filling))

# Randomly assign local states
initial_state = [0] * size
all_sites = list(range(size))  # To avoid having two types on same site.
for state, fill in zip(states, filling):
sites = random.sample(set(all_sites),
int(fill * size))  # pick fill*size sites to put state
for site in sites:
initial_state[site] = state
all_sites.remove(site)

return initial_state
`