include nkeep in mps.__init__

src/qailo/dispatch.py

@@ -3,21 +3,21 @@
 from . import state_vector as sv
 
 
-def apply(v, p, pos=None, maxdim=None):
+def apply(v, p, pos=None):
     if sv.is_state_vector(v):
         v = sv.apply(v, p, pos)
     elif mps.is_mps(v):
-        v = mps.apply(v, p, pos, maxdim)
+        v = mps.apply(v, p, pos)
     else:
         assert False
     return v
 
 
-def apply_seq(v, seq, maxdim=None):
+def apply_seq(v, seq):
     if sv.is_state_vector(v):
         v = sv.apply_seq(v, seq)
     elif mps.is_mps(v):
-        v = mps.apply_seq(v, seq, maxdim)
+        v = mps.apply_seq(v, seq)
     else:
         assert False
     return v

src/qailo/mps/apply.py

@@ -5,29 +5,29 @@
 from . import type as mps
 
 
-def _swap_tensors(m, s, maxdim=None):
+def _swap_tensors(m, s):
     """
     swap neighboring two tensors at s and s+1
     """
     assert s in range(0, mps.num_qubits(m) - 1)
-    m._apply_two(swap(), s, maxdim=maxdim)
+    m._apply_two(swap(), s)
     p0, p1 = m.t2q[s], m.t2q[s + 1]
     m.q2t[p0], m.q2t[p1] = s + 1, s
     m.t2q[s], m.t2q[s + 1] = p1, p0
 
 
-def _move_qubit(m, p, s, maxdim=None):
+def _move_qubit(m, p, s):
     if m.q2t[p] != s:
         # print(f"moving qubit {p} at {m.q2t[p]} to {s}")
         for u in range(m.q2t[p], s):
             # print(f"swap tensors {u} and {u+1}")
-            _swap_tensors(m, u, maxdim=maxdim)
+            _swap_tensors(m, u)
         for u in range(m.q2t[p], s, -1):
             # print(f"swap tensors {u-1} and {u}")
-            _swap_tensors(m, u - 1, maxdim=maxdim)
+            _swap_tensors(m, u - 1)
 
 
-def _apply(m, p, pos=None, maxdim=None):
+def _apply(m, p, pos=None):
     assert op.is_operator(p)
     n = mps.num_qubits(m)
     if pos is None:
@@ -41,24 +41,24 @@ def _apply(m, p, pos=None, maxdim=None):
         assert ss[0] != ss[1]
         if ss[0] < ss[1]:
             _move_qubit(m, pos[1], ss[0] + 1)
-            m._apply_two(p, ss[0], maxdim=maxdim)
+            m._apply_two(p, ss[0])
         else:
             _move_qubit(m, pos[0], ss[1] + 1)
-            m._apply_two(p, ss[1], maxdim=maxdim, reverse=True)
+            m._apply_two(p, ss[1], reverse=True)
     else:
         raise ValueError
     return m
 
 
-def apply(m, p, pos=None, maxdim=None):
-    return _apply(deepcopy(m), p, pos, maxdim)
+def apply(m, p, pos=None):
+    return _apply(deepcopy(m), p, pos)
 
 
-def _apply_seq(m, seq, maxdim=None):
+def _apply_seq(m, seq):
     for p, qubit in seq:
-        _apply(m, p, qubit, maxdim)
+        _apply(m, p, qubit)
     return m
 
 
-def apply_seq(m, seq, maxdim=None):
-    return _apply_seq(deepcopy(m), seq, maxdim)
+def apply_seq(m, seq):
+    return _apply_seq(deepcopy(m), seq)

src/qailo/mps/mps_c.py

@@ -22,9 +22,10 @@ class MPS_C(MPS):
         tensors [cp(1)+1...n-1]: right canonical
     """
 
-    def __init__(self, tensors):
+    def __init__(self, tensors, nkeep=None):
         assert isinstance(tensors, list)
         self.tensors = deepcopy(tensors)
+        self.nkeep = nkeep
         n = len(self.tensors)
         self.q2t = list(range(n))
         self.t2q = list(range(n))
@@ -108,7 +109,7 @@ def _apply_one(self, p, s):
         self.cp[0] = min(self.cp[0], s)
         self.cp[1] = max(self.cp[1], s)
 
-    def _apply_two(self, p, s, maxdim=None, reverse=False):
+    def _apply_two(self, p, s, reverse=False):
         """
         apply 2-qubit operator on neighboring tensors at s and s+1
         """
@@ -121,6 +122,6 @@ def _apply_two(self, p, s, maxdim=None, reverse=False):
             t = np.einsum(t, [0, 4, 5, 3], p, [1, 2, 4, 5])
         else:
             t = np.einsum(t, [0, 4, 5, 3], p, [2, 1, 5, 4])
-        L, R = tensor_svd(t, [[0, 1], [2, 3]], nkeep=maxdim)
+        L, R = tensor_svd(t, [[0, 1], [2, 3]], nkeep=self.nkeep)
         self.tensors[s] = L
         self.tensors[s + 1] = R

src/qailo/mps/product_state.py

@@ -23,16 +23,16 @@ def tensor_decomposition(v, nkeep=None, tol=1e-12):
         return tensors
 
 
-def product_state(states, mps=MPS_C):
+def product_state(states, nkeep=None, mps=MPS_C):
     tensors = []
     for s in states:
-        tensors = tensors + tensor_decomposition(s)
-    return mps(tensors)
+        tensors = tensors + tensor_decomposition(s, nkeep)
+    return mps(tensors, nkeep)
 
 
-def zero(n=1, mps=MPS_C):
-    return product_state([sv_zero()] * n, mps)
+def zero(n=1, nkeep=None, mps=MPS_C):
+    return product_state([sv_zero()] * n, nkeep, mps)
 
 
-def one(n=1, mps=MPS_C):
-    return product_state([sv_one()] * n, mps)
+def one(n=1, nkeep=None, mps=MPS_C):
+    return product_state([sv_one()] * n, nkeep, mps)

src/qailo/mps/svd.py

@@ -20,7 +20,7 @@ def tensor_svd(T, partition, canonical="center", nkeep=None, tol=1e-12):
     m = np.einsum(T, partition[0] + partition[1]).reshape(
         np.prod(dims0), np.prod(dims1)
     )
-    S, U, V = compact_svd(m, nkeep, tol)
+    S, U, V = compact_svd(m, nkeep=nkeep, tol=tol)
     L = U
     R = V.conj().T
     if canonical == "center":

src/qailo/mps_p/mps_p.py

@@ -23,9 +23,10 @@ class MPS_P(MPS):
         tensors [cp(1)+1...n-1]: bottom canonical
     """
 
-    def __init__(self, tensors):
+    def __init__(self, tensors, nkeep=None):
         assert isinstance(tensors, list)
         self.tensors = deepcopy(tensors)
+        self.nkeep = nkeep
         n = len(self.tensors)
         self.q2t = list(range(n))
         self.t2q = list(range(n))
@@ -110,7 +111,7 @@ def _apply_one(self, p, s):
         self.cp[0] = min(self.cp[0], s)
         self.cp[1] = max(self.cp[1], s)
 
-    def _apply_two(self, p, s, maxdim=None, reverse=False):
+    def _apply_two(self, p, s, reverse=False):
         """
         apply 2-qubit operator on neighboring tensors, s and s+1
         """
@@ -126,6 +127,6 @@ def _apply_two(self, p, s, maxdim=None, reverse=False):
             t1 = np.einsum(t1, [0, 4, 3], p0, [2, 4, 1])
         tt0 = np.einsum(self.env[s], [0, 4], t0, [4, 1, 2, 3])
         tt1 = np.einsum(t1, [0, 1, 2, 4], self.env[s + 2], [4, 3])
-        _, WLh, WR = projector(tt0, [0, 1, 4, 5], tt1, [5, 4, 2, 3], maxdim)
+        _, WLh, WR = projector(tt0, [0, 1, 4, 5], tt1, [5, 4, 2, 3], nkeep=self.nkeep)
         self.tensors[s] = np.einsum(t0, [0, 1, 3, 4], WR, [3, 4, 2])
         self.tensors[s + 1] = np.einsum(WLh, [3, 4, 0], t1, [4, 3, 1, 2])

test/mps/test_apply.py

@@ -4,11 +4,11 @@
 from pytest import approx
 
 
-def apply(op, m0, m1, m2, m3, v, seq, pos, maxdim=None):
+def apply(op, m0, m1, m2, m3, v, seq, pos):
     m0 = q.mps.apply(m0, op, pos)
     m1 = q.mps.apply(m1, op, pos)
     m2 = q.mps.apply(m2, op, pos)
-    m3 = q.mps.apply(m3, op, pos, maxdim)
+    m3 = q.mps.apply(m3, op, pos)
     v = q.sv.apply(v, op, pos)
     seq.append([op, pos])
     return m0, m1, m2, m3, v, seq
@@ -17,19 +17,19 @@ def apply(op, m0, m1, m2, m3, v, seq, pos, maxdim=None):
 def test_apply():
     n = 8
     p = 64
-    maxdim = 2
+    nkeep = 2
 
-    m0 = q.mps.zero(n, q.mps.MPS_C)
-    m1 = q.mps.zero(n, q.mps_p.MPS_P)
-    m2 = q.mps.zero(n, q.mps.MPS_C)
-    m3 = q.mps.zero(n, q.mps_p.MPS_P)
+    m0 = q.mps.zero(n, mps=q.mps.MPS_C)
+    m1 = q.mps.zero(n, mps=q.mps_p.MPS_P)
+    m2 = q.mps.zero(n, nkeep=nkeep, mps=q.mps.MPS_C)
+    m3 = q.mps.zero(n, nkeep=nkeep, mps=q.mps_p.MPS_P)
     v = q.sv.zero(n)
     seq = []
 
     i = 4
     j = 0
     print("apply cz on {} and {}".format(i, j))
-    m0, m1, m2, m3, v, seq = apply(q.op.cz(), m0, m1, m2, m3, v, seq, [i, j], maxdim)
+    m0, m1, m2, m3, v, seq = apply(q.op.cz(), m0, m1, m2, m3, v, seq, [i, j])
 
     for _ in range(p):
         i = random.randrange(n)
@@ -38,30 +38,24 @@ def test_apply():
             t = random.randrange(3)
             if t == 0:
                 print("apply h on {}".format(i))
-                m0, m1, m2, m3, v, seq = apply(
-                    q.op.h(), m0, m1, m2, m3, v, seq, [i], maxdim
-                )
+                m0, m1, m2, m3, v, seq = apply(q.op.h(), m0, m1, m2, m3, v, seq, [i])
             elif t == 1:
                 print("apply x on {}".format(i))
-                m0, m1, m2, m3, v, seq = apply(
-                    q.op.s(), m0, m1, m2, m3, v, seq, [i], maxdim
-                )
+                m0, m1, m2, m3, v, seq = apply(q.op.s(), m0, m1, m2, m3, v, seq, [i])
             elif t == 2:
                 print("apply s on {}".format(i))
-                m0, m1, m2, m3, v, seq = apply(
-                    q.op.t(), m0, m1, m2, m3, v, seq, [i], maxdim
-                )
+                m0, m1, m2, m3, v, seq = apply(q.op.t(), m0, m1, m2, m3, v, seq, [i])
         else:
             t = random.randrange(2)
             if t == 0:
                 print("apply cx on {} and {}".format(i, j))
                 m0, m1, m2, m3, v, seq = apply(
-                    q.op.cx(), m0, m1, m2, m3, v, seq, [i, j], maxdim
+                    q.op.cx(), m0, m1, m2, m3, v, seq, [i, j]
                 )
             elif t == 1:
                 print("apply cz on {} and {}".format(i, j))
                 m0, m1, m2, m3, v, seq = apply(
-                    q.op.cz(), m0, m1, m2, m3, v, seq, [i, j], maxdim
+                    q.op.cz(), m0, m1, m2, m3, v, seq, [i, j]
                 )
         m0._is_canonical()
         m1._is_canonical()
@@ -74,7 +68,6 @@ def test_apply():
         print("fidelity = {} {} {} {}".format(f0, f1, f2, f3))
         assert f0 == approx(1)
         assert f1 == approx(1)
-        assert f2 == approx(1)
 
     f0 = q.sv.fidelity(q.mps.state_vector(m0), v)
     f1 = q.sv.fidelity(q.mps.state_vector(m1), v)
@@ -83,10 +76,8 @@ def test_apply():
     print("final fidelity = {} {} {} {}".format(f0, f1, f2, f3))
     assert f0 == approx(1)
     assert f1 == approx(1)
-    assert f2 == approx(1)
-    # assert f2 == approx(f3)
 
-    m4 = q.mps.zero(n, q.mps_p.MPS_P)
+    m4 = q.mps.zero(n, mps=q.mps_p.MPS_P)
     f4 = q.sv.fidelity(q.mps.state_vector(q.mps.apply_seq(m4, seq)), v)
     assert f4 == approx(1)