Merge branch 'master' of github.com:dseuss/mpnum

mpnum/_tools.py

@@ -13,7 +13,7 @@
 from six.moves import range, zip
 
 
-def global_to_local(array, sites):
+def global_to_local(array, sites, left_skip=0, right_skip=0):
     """Converts a general `sites`-local array with fixed number p of physical
     legs per site from the global form
 
@@ -27,6 +27,8 @@ def global_to_local(array, sites):
 
     :param np.ndarray array: Array with ndim, such that ndim % sites = 0
     :param int sites: Number of distinct sites
+    :param int left_skip: Ignore that many axes on the left
+    :param int right_skip: Ignore that many axes on the right
     :returns: Array with same ndim as array, but reshaped
 
     >>> global_to_local(np.zeros((1, 2, 3, 4, 5, 6)), 3).shape
@@ -35,10 +37,14 @@ def global_to_local(array, sites):
     (1, 3, 5, 2, 4, 6)
 
     """
-    assert array.ndim % sites == 0, \
+    skip = left_skip + right_skip
+    ndim = array.ndim - skip
+    assert ndim % sites == 0, \
         "ndim={} is not a multiple of {}".format(array.ndim, sites)
-    plegs = array.ndim // sites
-    order = [i // plegs + sites * (i % plegs) for i in range(plegs * sites)]
+    plegs = ndim // sites
+    order = (left_skip + i + sites * j for i in range(sites) for j in range(plegs))
+    order = tuple(it.chain(
+        range(left_skip), order, range(array.ndim - right_skip, array.ndim)))
     return np.transpose(array, order)
 
 
@@ -73,7 +79,7 @@ def local_to_global(array, sites, left_skip=0, right_skip=0):
     plegs = ndim // sites
     order = (left_skip + plegs*i + j for j in range(plegs) for i in range(sites))
     order = tuple(it.chain(
-        range(left_skip), order, range(array.ndim-right_skip, array.ndim)))
+        range(left_skip), order, range(array.ndim - right_skip, array.ndim)))
     return np.transpose(array, order)
 
 

mpnum/linalg.py

@@ -186,7 +186,7 @@ def _mineig_local_op(leftvec, mpo_ltens, rightvec):
 
 
 def _mineig_minimize_locally(leftvec, mpo_ltens, rightvec, eigvec_ltens,
-                             eigs_opts=None):
+                             user_eigs_opts=None):
     """Perform the local eigenvalue minimization on one site on one site.
 
     Return a new (expectedly smaller) eigenvalue and a new local
@@ -213,11 +213,12 @@ def _mineig_minimize_locally(leftvec, mpo_ltens, rightvec, eigvec_ltens,
     Middle row: MPO matrices with row (column) indices to bottom (top)
 
     """
-    if eigs_opts is None:
-        eigs_opts = {'k': 1, 'which': 'SR', 'tol': 1e-6}
-    else:
-        eigs_opts = eigs_opts.copy()
-        eigs_opts['k'] = 1
+    eigs_opts = {'k': 1, 'which': 'SR', 'tol': 1e-6}
+    if user_eigs_opts is not None:
+        eigs_opts.update(user_eigs_opts)
+    if eigs_opts['k'] != 1:
+        raise ValueError('Supplying k != 1 in requires changes in the code, '
+                         'k={} was requested'.format(user_eigs_opts['k']))
     op = _mineig_local_op(leftvec, mpo_ltens, rightvec)
     eigvec_bonddim = max(lten.shape[0] for lten in eigvec_ltens)
     eigvec_lten = eigvec_ltens[0]
@@ -257,8 +258,9 @@ def mineig(mpo,
         no start vector is given. (default: Use the bond dimension of `mpo`)
     :param randstate: numpy.random.RandomState instance or None
     :param max_num_sweeps: Maximum number of sweeps to do (default 5)
-    :param eigs_opts: kwargs for scipy.sparse.linalg.eigs(), k is always set
-        equal to 1!
+    :param eigs_opts: kwargs for `scipy.sparse.linalg.eigs()`. If you
+        supple `which`, you will probably not obtain the minimal
+        eigenvalue. `k` different from one is not supported at the moment.
     :param int minimize_sites: Number of connected sites minimization should
         be performed on (default 1)
     :returns: mineigval, mineigval_eigvec_mpa
@@ -306,10 +308,20 @@ def mineig(mpo,
         pdims = max(dim[0] for dim in mpo.pdims)
         if startvec_bonddim is None:
             startvec_bonddim = max(mpo.bdims)
+        if startvec_bonddim == 1:
+            raise ValueError('startvec_bonddim must be at least 2')
 
         startvec = random_mpa(nr_sites, pdims, startvec_bonddim,
                               randstate=randstate)
         startvec /= mp.norm(startvec)
+    # Can we avoid this overly complex check by improving
+    # _mineig_minimize_locally()? eigs() will fail under the excluded
+    # conditions because of too small matrices.
+    assert not any(bdim12 == (1, 1) for bdim12 in
+                   zip((1,) + startvec.bdims, startvec.bdims + (1,))), \
+        'startvec must not contain two consecutive bonds of dimension 1, ' \
+        'bdims including dummy bonds = (1,) + {!r} + (1,)' \
+            .format(startvec.bdims)
     # For
     #
     #   pos in range(nr_sites - minimize_sites),

tests/conftest.py

@@ -5,3 +5,9 @@
 @pytest.fixture(scope="module")
 def rgen():
     return numpy.random.RandomState(seed=52973992)
+
+
+@pytest.fixture(scope='function', autouse=True)
+def bug_workaround():
+    # Workaround for https://github.com/pytest-dev/pytest/issues/1832
+    pass

tests/linalg_test.py

@@ -7,6 +7,7 @@
 import pytest as pt
 from numpy.testing import assert_almost_equal
 
+import mpnum as mp
 import mpnum.linalg
 import mpnum.factory as factory
 
@@ -67,3 +68,30 @@ def test_mineig_minimize_sites(nr_sites, local_dim, bond_dim, rgen):
     overlap = np.inner(mineig_eigvec.conj(), mineig_eigvec_mp)
     assert_almost_equal(mineig, mineig_mp)
     assert_almost_equal(1, abs(overlap))
+
+
+@pt.mark.parametrize('nr_sites, local_dim, bond_dim', MP_TEST_PARAMETERS)
+def test_mineig_eigs_opts(nr_sites, local_dim, bond_dim, rgen):
+    """Verify correct operation if eigs_opts() is specified
+
+    This test mainly verifies correct operation if the user specifies
+    eigs_opts() but does not include which='SR' herself. It also tests
+    minimization on another example MPO (a rank-1 MPO in this case).
+
+    """
+    # Need at least two sites
+    if nr_sites < 2:
+        return
+
+    mps = factory.random_mps(nr_sites, local_dim, bond_dim, rgen)
+    mpo = mp.mps_to_mpo(mps)
+    # mineig does not support startvec_bonddim = 1
+    bond_dim = 2 if bond_dim == 1 else bond_dim
+    eigval, eigvec = mpnum.linalg.mineig(
+        mpo, startvec_bonddim=bond_dim, randstate=rgen, max_num_sweeps=10,
+        eigs_opts=dict(tol=1e-5), minimize_sites=1)
+    # Check correct eigenvalue
+    assert_almost_equal(eigval, 0)
+    # Check for orthogonal eigenvectors and `eigvec` being in the
+    # kernel of `mpo`
+    assert_almost_equal(abs(mp.inner(eigvec, mps)), 0)