add zero one for mps

example/grover.py

@@ -27,9 +27,9 @@ def diffusion(v):
 
 def grover(n, target, iter, use_mps):
     if use_mps:
-        v = q.mps.MPS_P(q.mps.product_state(n))
+        v = q.mps.zero(n)
     else:
-        v = q.sv.state_vector(n)
+        v = q.sv.zero(n)
     for i in range(n):
         v = q.apply(v, q.op.h(), [i])
     for k in range(iter):

example/qpe.py

@@ -1,13 +1,12 @@
-from copy import deepcopy
-
 import numpy as np
 import qailo as q
 
 
 def qpe(n, u, ev):
     m = q.num_qubits(u)
-    assert q.num_qubits(ev) == n + m
-    v = deepcopy(ev)
+    assert q.num_qubits(ev) == m
+    v = q.sv.product_state([q.sv.zero(n), ev])
+    assert q.num_qubits(v) == n + m
 
     for p in range(n):
         v = q.apply(v, q.op.h(), [p])
@@ -28,8 +27,8 @@ def qpe(n, u, ev):
     n = 3
     phi = 2 * np.pi * (1 / 3)
     u = q.op.p(phi)
-    ev = q.sv.state_vector(n + 1)
-    ev = q.apply(ev, q.op.x(), [n])
+    ev = q.sv.zero()
+    ev = q.apply(ev, q.op.x())
     v = qpe(n, u, ev)
     prob = np.diag(q.op.matrix(q.op.trace(q.sv.pure_state(v), [n])).real)
     for i in range(len(prob)):

example/simple.py

@@ -3,9 +3,9 @@
 
 def main(use_mps=False):
     if use_mps:
-        v = q.mps.MPS_C(q.mps.tensor_decomposition(q.sv.state_vector(3)))
+        v = q.mps.zero(3)
     else:
-        v = q.sv.state_vector(3)
+        v = q.sv.zero(3)
     print("input:")
     print("state vector:", q.vector(v))
     print("probabitily:", q.probability(v))

example/test_grover.py

@@ -3,7 +3,7 @@
 from pytest import approx
 
 
-def test_simple():
+def test_grover():
     n = 4
     target = q.util.str2binary("0000")
     iter = 2 ** (n // 2)

example/test_qpe.py

@@ -1,16 +1,16 @@
 import numpy as np
 import qailo as q
-import qpe
 from pytest import approx
+from qpe import qpe
 
 
 def test_qpe():
     n = 3
     phi = 2 * np.pi * (1 / 8)
     u = q.op.p(phi)
-    ev = q.sv.state_vector(n + 1)
-    ev = q.apply(ev, q.op.x(), [n])
-    v = qpe.qpe(n, u, ev)
+    ev = q.sv.zero()
+    ev = q.apply(ev, q.op.x())
+    v = qpe(n, u, ev)
     prob = np.diag(q.op.matrix(q.op.trace(q.sv.pure_state(v), [n])).real)
     assert prob[4] == approx(1)
     assert prob[0] == approx(0)

src/qailo/mps/__init__.py

@@ -2,7 +2,7 @@
 from .mps_c import MPS_C
 from .mps_p import MPS_P
 from .norm import norm
-from .product_state import product_state, tensor_decomposition
+from .product_state import one, product_state, tensor_decomposition, zero
 from .projector import projector
 from .state_vector import state_vector
 from .svd import compact_svd, tensor_svd
@@ -14,8 +14,10 @@
     MPS_C,
     MPS_P,
     norm,
+    one,
     product_state,
     tensor_decomposition,
+    zero,
     projector,
     state_vector,
     compact_svd,

src/qailo/mps/product_state.py

@@ -1,28 +1,34 @@
-import numpy as np
-
+from ..state_vector.state_vector import one as sv_one
+from ..state_vector.state_vector import zero as sv_zero
 from ..state_vector.type import num_qubits
 from ..state_vector.vector import vector
+from .mps_c import MPS_C
 from .svd import tensor_svd
 
 
-def product_state(n, c=0):
-    assert n > 0
-    tensors = []
-    for t in range(n):
-        tensor = np.zeros((1, 2, 1))
-        tensor[0, (c >> (n - t - 1)) & 1, 0] = 1
-        tensors.append(tensor)
-    return tensors
-
-
 def tensor_decomposition(v, nkeep=None, tol=1e-12):
     n = num_qubits(v)
-    vv = vector(v).reshape((1, 2**n))
+    w = vector(v).reshape((1, 2**n))
     tensors = []
     for t in range(n - 1):
-        dims = vv.shape
-        vv = vv.reshape(dims[0], 2, dims[1] // 2)
-        t, vv = tensor_svd(vv, [[0, 1], [2]], "left", nkeep, tol)
+        dims = w.shape
+        w = w.reshape(dims[0], 2, dims[1] // 2)
+        t, w = tensor_svd(w, [[0, 1], [2]], "left", nkeep, tol)
         tensors.append(t)
-    tensors.append(vv.reshape(vv.shape + (1,)))
+    tensors.append(w.reshape(w.shape + (1,)))
     return tensors
+
+
+def product_state(states, mps=MPS_C):
+    tensors = []
+    for s in states:
+        tensors = tensors + tensor_decomposition(s)
+    return mps(tensors)
+
+
+def zero(n, mps=MPS_C):
+    return product_state([sv_zero()] * n, mps)
+
+
+def one(n, mps=MPS_C):
+    return product_state([sv_one()] * n, mps)

src/qailo/state_vector/__init__.py

@@ -3,7 +3,7 @@
 from .is_close import is_close
 from .probability import probability
 from .pure_state import pure_state
-from .state_vector import one, product_state, state_vector, zero
+from .state_vector import one, product_state, zero
 from .type import is_state_vector, num_qubits
 from .vector import vector
 
@@ -16,7 +16,6 @@
     pure_state,
     one,
     product_state,
-    state_vector,
     zero,
     is_state_vector,
     num_qubits,

src/qailo/state_vector/state_vector.py

@@ -1,26 +1,36 @@
 import numpy as np
 
-
-def zero():
-    return np.array((1, 0)).reshape((2, 1))
-
-
-def one():
-    return np.array((0, 1)).reshape((2, 1))
+from .type import is_state_vector, num_qubits
 
 
 def product_state(states):
-    n = len(states)
-    assert n > 0
+    m = len(states)
+    assert m > 0
     v = states[0]
-    print(0, v.shape)
-    for i in range(1, n):
-        v = np.einsum(v.reshape((2,) * (i)), list(range(i)), states[i], [i + 1, i + 2])
-        print(i, v.shape)
+    assert is_state_vector(v)
+    for i in range(1, m):
+        n0 = num_qubits(v)
+        n1 = num_qubits(states[i])
+        v = np.einsum(
+            v.reshape((2,) * n0),
+            list(range(n0)),
+            states[i],
+            list(range(n0, n0 + n1 + 1)),
+        )
     return v
 
 
-def state_vector(n, c=0):
-    v = np.zeros(2**n)
-    v[c] = 1
-    return v.reshape((2,) * n + (1,))
+def zero(n=1):
+    assert n > 0
+    if n == 1:
+        return np.array((1, 0)).reshape((2, 1))
+    else:
+        return product_state([zero(1)] * n)
+
+
+def one(n=1):
+    assert n > 0
+    if n == 1:
+        return np.array((0, 1)).reshape((2, 1))
+    else:
+        return product_state([one(1)] * n)

test/alg/test_qft.py

@@ -8,7 +8,7 @@ def test_qft():
     target = 5
 
     # generate result of qft of target
-    v = q.sv.state_vector(n)
+    v = q.sv.zero(n)
     for p in range(n):
         v = q.apply(v, q.op.h(), [p])
     v = q.apply(v, q.op.p(target * np.pi / 4), [0])

test/mps/test_apply.py

@@ -2,6 +2,8 @@
 
 import qailo as q
 from pytest import approx
+from qailo.mps.mps_c import MPS_C
+from qailo.mps.mps_p import MPS_P
 
 
 def apply(op, m0, m1, m2, m3, v, seq, pos, maxdim=None):
@@ -19,11 +21,11 @@ def test_apply():
     p = 64
     maxdim = 2
 
-    m0 = q.mps.MPS_C(q.mps.product_state(n))
-    m1 = q.mps.MPS_P(q.mps.product_state(n))
-    m2 = q.mps.MPS_C(q.mps.product_state(n))
-    m3 = q.mps.MPS_P(q.mps.product_state(n))
-    v = q.sv.state_vector(n)
+    m0 = q.mps.zero(n, MPS_C)
+    m1 = q.mps.zero(n, MPS_P)
+    m2 = q.mps.zero(n, MPS_C)
+    m3 = q.mps.zero(n, MPS_P)
+    v = q.sv.zero(n)
     seq = []
 
     i = 4
@@ -84,7 +86,7 @@ def test_apply():
     assert f1 == approx(1)
     # assert f2 == approx(f3)
 
-    m4 = q.mps.MPS_P(q.mps.product_state(n))
+    m4 = q.mps.zero(n, MPS_P)
     f4 = q.sv.fidelity(q.mps.state_vector(q.mps.apply_seq(m4, seq)), v)
     assert f4 == approx(1)
 

test/mps/test_mps.py

@@ -1,16 +1,18 @@
 import numpy as np
 import qailo as q
+from qailo.mps.mps_c import MPS_C
+from qailo.mps.mps_p import MPS_P
 
 
 def test_mps():
-    n = 4
-    c = q.util.str2binary("1100")
-    m0 = q.mps.MPS_C(q.mps.product_state(n, c))
-    m1 = q.mps.MPS_P(q.mps.product_state(n, c))
+    states = [q.sv.one(), q.sv.one(), q.sv.zero(), q.sv.zero()]
+    v = q.sv.product_state(states)
+    m0 = q.mps.product_state(states, MPS_C)
+    m1 = q.mps.product_state(states, MPS_P)
     v0 = q.mps.state_vector(m0)
     v1 = q.mps.state_vector(m1)
-    assert np.allclose(v0, q.sv.state_vector(n, c))
-    assert np.allclose(v1, q.sv.state_vector(n, c))
+    assert np.allclose(v0, v)
+    assert np.allclose(v1, v)
 
 
 if __name__ == "__main__":

test/operator/test_swap.py

@@ -12,9 +12,9 @@ def test_swap():
     p = q.op.multiply(p, q.op.cx(), [0, 1])
     assert q.op.is_close(p, q.op.swap())
 
-    v = q.sv.state_vector(2)
+    v = q.sv.zero(2)
     v = q.sv.apply(v, q.op.h(), [0])
     v = q.sv.apply(v, q.op.swap(), [0, 1])
     v = q.sv.apply(v, q.op.h(), [1])
     v = q.sv.apply(v, q.op.swap(), [0, 1])
-    assert q.op.is_close(v, q.sv.state_vector(2))
+    assert q.op.is_close(v, q.sv.zero(2))

test/state_vector/test_apply_seq.py

@@ -4,9 +4,9 @@
 
 
 def test_apply_seq():
-    for n in range(8):
-        v0 = q.sv.state_vector(n)
-        v1 = q.sv.state_vector(n)
+    for n in range(1, 8):
+        v0 = q.sv.zero(n)
+        v1 = q.sv.zero(n)
         seq = []
         for i in range(n):
             v1 = q.apply(v1, q.op.h(), [i])

test/state_vector/test_fidelity.py

@@ -4,8 +4,8 @@
 
 
 def test_fidelity():
-    for n in range(8):
-        sv0 = q.sv.state_vector(n)
+    for n in range(1, 8):
+        sv0 = q.sv.zero(n)
         sv1 = sv0
         for i in range(n):
             sv1 = q.sv.apply(sv1, q.op.h(), [i])

test/state_vector/test_pure_state.py

@@ -3,7 +3,7 @@
 
 def test_pure_state():
     n = 2
-    sv = q.sv.state_vector(n)
+    sv = q.sv.product_state([q.sv.zero()] * n)
     dm = q.sv.pure_state(sv)
     assert q.op.is_density_matrix(dm)
 

test/state_vector/test_state_vector.py

@@ -15,6 +15,8 @@ def test_product_state():
     n = 8
     v0 = q.sv.product_state([q.sv.zero()] * n)
     v1 = q.sv.product_state([q.sv.one()] * n)
+    assert np.allclose(v0, q.sv.zero(n))
+    assert np.allclose(v1, q.sv.one(n))
     for p in range(n):
         v1 = q.sv.apply(v1, q.op.x(), [p])
     assert np.allclose(v0, v1)