Qbit added with quantum gate applicable [FIX: a little bug on how the…

… gate is performed (seems to invert in some way) ]

quantum/gate.py

@@ -24,4 +24,11 @@
 H = 1/np.sqrt(2) * np.array([
     [1., 1.],
     [1., -1.]
-])
+])
+
+CX = np.array([
+    [1., 0., 0., 0.],
+    [0., 1., 0., 0.],
+    [0., 0., 0., 1.],
+    [0., 0., 1., 0.],
+]).reshape((2,2,2,2))

quantum/qbit.py

@@ -1,20 +1,59 @@
 
 
+from matplotlib.pyplot import isinteractive
+from syngular.tensor.matrix_product_operator import MatrixProductOperator
 from syngular.tensor.matrix_product_state import MatrixProductState
 from syngular.quantum import gate
 
+from typing import Tuple, List, Union
+import numpy as np
+
 class Qbit:
 
-    def __init__(self, size):
+    def __init__(self, size, init=True):
         self.size = size
         self.dim = 2**self.size
 
-        self.state = MatrixProductState.zeros((2,)*size, (2,)*(size-1)).decompose()
-        self.state.real_parameters_number = self.dim
+        self.state = None
 
-        for idx in range(self.state.sites_number-1):
-            self.state.sites[idx][0] = gate.I
-        self.state.sites[-1][0][0] = 1.
+        if init:
+            self.state = MatrixProductState.zeros((2,)*size, (2,)*(size-1)).decompose()
+            self.state.real_parameters_number = self.dim
+
+            for idx in range(self.state.sites_number-1):
+                self.state.sites[idx][0] = gate.I
+            self.state.sites[-1][0][0] = 1.
+
+    def __matmul__(self, operator: Union[Tuple, Union[List, MatrixProductOperator]]):
+        if isinstance(operator, tuple):
+            return Qbit.from_mps(self.state.apply(*operator))
+        else:
+            operator = MatrixProductOperator(operator, bond_shape=()).decompose()
+            return Qbit.from_mps(operator @ self.state)
+
+    def apply(self, gate):
+        gate = MatrixProductOperator(gate, bond_shape=()).decompose()
+        return Qbit.from_mps(gate @ self.state)
 
     def to_tensor(self):
-        return self.state.to_tensor().reshape(self.dim)
+        return self.state.to_tensor().reshape(self.dim)
+
+    def to_binary(self):
+        tensor = self.to_tensor().astype(int)
+        # print(tensor)
+        # pos = 0
+        # value = self.state[(0,) * self.size]
+        # while value != 1 and pos < 2**self.size-1:
+        #     pos += 1
+        #     print(pos, 2**self.size-1)
+        #     print((0,) * pos + (1,) + (0,) * (self.size-pos-1))
+        #     value = self.state[(0,) * pos + (1,) + (0,) * (self.size-pos-1)]
+            # print((0,) * pos + (1,) + (0,) * (self.size-pos-1))
+            # print(value)
+        return bin(np.where(tensor == 1)[0][0])[2:].zfill(self.size)
+
+    def from_mps(mps):
+        # print(mps.sites_number)
+        qbit = Qbit(mps.sites_number, init=False)
+        qbit.state = mps
+        return qbit

tensor/matrix_product_operator.py

@@ -4,6 +4,7 @@
 import gc
 import itertools
 from typing import Tuple, List, Union
+from numpy.core.numeric import outer
 
 from opt_einsum import contract
 
@@ -29,20 +30,21 @@ def __init__(self, tensor=None, bond_shape: Union[Union[Tuple, List], np.ndarray
             self.sites_number = len(self.bond_shape)+1
             self.sites = [None] * self.sites_number
 
-            if self.bond_shape == (): self.bond_shape = (1,)
-
             self.input_shape    = tuple(self.tensor_shape[:self.sites_number])
             self.output_shape   = tuple(self.tensor_shape[self.sites_number:])
-
+    
             if len(self.input_shape) != len(self.output_shape):
                 raise Exception("input_shape and output_shape of the tensor must have the same length")
 
             if len(self.input_shape) != len(self.bond_shape)+1:
                 raise Exception("dimensions of bond indices do not match input dimension - 1")
 
-            self.shape = [(1, self.tensor_shape[0], self.tensor_shape[self.sites_number], self.bond_shape[0])]
-            self.shape += [(self.bond_shape[i-1], self.tensor_shape[i], self.tensor_shape[self.sites_number+i], self.bond_shape[i]) for i in range(1, self.sites_number-1)]
-            self.shape += [(self.bond_shape[self.sites_number-2], self.tensor_shape[self.sites_number-1], self.tensor_shape[2*self.sites_number-1], 1)]
+            if self.bond_shape != ():
+                self.shape = [(1, self.tensor_shape[0], self.tensor_shape[self.sites_number], self.bond_shape[0])]
+                self.shape += [(self.bond_shape[i-1], self.tensor_shape[i], self.tensor_shape[self.sites_number+i], self.bond_shape[i]) for i in range(1, self.sites_number-1)]
+                self.shape += [(self.bond_shape[self.sites_number-2], self.tensor_shape[self.sites_number-1], self.tensor_shape[2*self.sites_number-1], 1)]
+            else:
+                self.shape = [(1, self.tensor_shape[0], self.tensor_shape[1], 1)]
 
             self.shape_indices = [
                 (2*self.sites_number+i, i, self.sites_number+i, 2*self.sites_number+i+1) for i in range(self.sites_number)
@@ -296,7 +298,14 @@ def to_tensor(self):
         return tensor
 
     def decompose(self, mode="left"):
-
+        if self.bond_shape == ():
+            self.sites = [self.tensor.reshape(1, self.tensor_shape[0], self.tensor_shape[1], 1)]
+
+            self.decomposed = True
+
+            del self.tensor
+            gc.collect()
+
         if not self.decomposed:
             if mode == "left":
 

tensor/matrix_product_state.py

@@ -43,12 +43,14 @@ def __init__(self, tensor=None, bond_shape: Union[Union[Tuple, List], np.ndarray
 
             if self.order != self.sites_number:
                 raise Exception("dimensions of bond indices do not match order - 1")
-
-
-            self.shape = [(1, self.tensor_shape[0],self.bond_shape[0])]
-            self.shape += [(self.bond_shape[i-1], self.tensor_shape[i], self.bond_shape[i]) for i in range(1, self.sites_number-1)]
-            self.shape += [(self.bond_shape[self.sites_number-2], self.tensor_shape[self.sites_number-1], 1)]
-
+
+            if self.bond_shape != ():
+                self.shape = [(1, self.tensor_shape[0],self.bond_shape[0])]
+                self.shape += [(self.bond_shape[i-1], self.tensor_shape[i], self.bond_shape[i]) for i in range(1, self.sites_number-1)]
+                self.shape += [(self.bond_shape[self.sites_number-2], self.tensor_shape[self.sites_number-1], 1)]
+            else:
+                self.shape = [(1, self.tensor_shape[0], 1)]
+
         self.verbose = verbose
         self.decomposed = False
 
@@ -185,6 +187,7 @@ def from_sites(sites: List[np.ndarray], orthogonality=None, real_parameters_numb
         mp = MatrixProductState()
         mp.sites = sites
         mp.sites_number = len(sites)
+        # mp.order = len()
 
         mp.decomposed = True
         mp.orthonormalized = orthogonality
@@ -213,11 +216,13 @@ def empty():
     def zeros(input_shape, bond_shape):
         n = len(input_shape)
 
-        shape = [(1, input_shape[0],bond_shape[0])]
-        shape += [(bond_shape[i-1], input_shape[i], bond_shape[i]) for i in range(1, n-1)]
-        shape += [(bond_shape[n-2], input_shape[n-1], 1)]
-
-        print(shape)
+        if bond_shape != ():
+            shape = [(1, input_shape[0],bond_shape[0])]
+            shape += [(bond_shape[i-1], input_shape[i], bond_shape[i]) for i in range(1, n-1)]
+            shape += [(bond_shape[n-2], input_shape[n-1], 1)]
+        else:
+            shape = [(1, input_shape[0], 1)]
+        # print(shape)
         sites = []
         for idx in range(n):
             sites.append(np.zeros(shape=shape[idx]))
@@ -240,12 +245,12 @@ def norm(self):
     @returns: the corresponding numpy tensor
     """
     def to_tensor(self):
-        tensor = np.zeros(shape=self.input_shape)
+        tensor = np.zeros(shape=self.input_shape, dtype=np.complex)
 
         range_inp = [range(inp) for inp in self.input_shape]
 
         for inp in itertools.product(*range_inp):
-                tensor[inp] = self[inp]
+            tensor[inp] = self[inp]
 
         return tensor
 
@@ -396,6 +401,50 @@ def compress(self, dim: int, mode="left", strict=False):
 
             return that
 
+    def apply(self, operator, index, strict=True):
+        if strict: that = self.copy()
+        # print("> Apply")
+        m = len(operator.shape) // 2
+        n = that.sites_number
+
+        # index = n-1-index
+
+        struct = []
+        struct += [operator, [*list(range(2*n+index, 2*n+m+index)), *list(range(3*n+index, 3*n+m+index))]]
+        for idx in range(index, m+index):
+            struct += [that.sites[idx], [idx, 2*n+idx, idx+1]]
+        struct += [[index, *list(range(3*n+index, 3*n+m+index)), m+index]]
+        # print(struct)
+
+        T = contract(*struct)
+        # print('T', T.shape)
+
+        # print(index, m, m+index-1, list(range(index, m+index-1)))
+        for k in range(index, m+index-1):
+            print("K", k)
+            lrank = T.shape[0]
+            input_shape = T.shape[1]
+
+            L = T.reshape(input_shape*lrank, -1)
+
+            Q, R = np.linalg.qr(L, mode="complete")
+            # print(Q.shape, R.shape)
+            rank = that.shape[k][2]
+            Q = Q[:,:rank]
+            R = R[:rank, :]
+
+            that.sites[k] = that.tensoricization(Q, k)
+            T = R
+
+        # print(that.sites[m+index-1])
+        # print('---')
+        # print(T)
+        # print(that.sites[m+index-1].shape)
+        that.sites[m+index-1] = that.tensoricization(T, m+index-1)
+        # print(that.sites[m+index-1].shape)
+
+        return that
+
     def retrieve(self, indices):
         einsum_structure = []
         for idx in range(self.sites_number):

test/test_mps.py

@@ -6,6 +6,18 @@
 
 np.set_printoptions(suppress=True)
 
+def test_apply():
+    x = np.arange(8).reshape((2,2,2))
+    X = MatrixProductState(x, bond_shape=(2,2)).decompose()
+
+    from syngular.quantum import gate
+
+    g = gate.CX.reshape((2,2,2,2))
+    # print(g)
+    X.apply(g, 1)
+
+    print(X)
+
 
 def test_augment():
     x = np.arange(8).reshape((2,2,2))
@@ -90,4 +102,5 @@ def test_random():
 # test_random()
 # test_zeros()
 
-test_augment()
+# test_augment()
+test_apply()

test/test_quantum.py

@@ -1,17 +1,102 @@
 from syngular.quantum import Circuit, Qbit
 import syngular.quantum.gate as gate
 
+import time
+
 circ = Circuit(initializer="ground", size=2, structure=[
     (gate.X, 0)
 ])
 
-size = 30
+size = 1
+
+# ground = Qbit(size)
+# print("> Ground state")
+# print(ground.to_tensor())
+
+# terminal = ground @ (1.j * gate.X) @ (-1.j * gate.X)
+# print("> Terminal state")
+# print(terminal.to_tensor())
+
+
+size = 4
 
 ground = Qbit(size)
 
+print("> Ground state")
+print(ground.to_tensor())
+
+# ground = ground @ gate.X
+
+print("> Intermediate state : X-gate:0 & X-gate:1")
+intermediate = ground @ (gate.X, 1) @ (gate.X, 0) @ (gate.X, 2) @ (gate.X, 3)
+# print(intermediate.to_tensor())
+print(intermediate.to_binary())
+print("> Intermediate state : X-gate:1")
+intermediate = intermediate @ (gate.X, 1)
+# print(intermediate.to_tensor())
+print(intermediate.to_binary())
+print("> Intermediate state : X-gate:1")
+intermediate = intermediate @ (gate.X, 1)
+# print(intermediate.to_tensor())
+print(intermediate.to_binary())
+print("> Intermediate state : X-gate:0")
+intermediate = intermediate @ (gate.X, 0)
+# print(intermediate.to_tensor())
+print(intermediate.to_binary())
+print("> Intermediate state : X-gate:1")
+intermediate = intermediate @ (gate.X, 1)
+# print(intermediate.to_tensor())
+print(intermediate.to_binary())
+# ground.state.apply(gate.X, 1)
 # print(ground.to_tensor())
-print(ground.state[(1,)*size])
-print(ground.state[(0,)*size])
+print("> Intermediate state : X-gate:3")
+intermediate = intermediate @ (gate.X, 3)
+# print(intermediate.to_tensor())
+print(intermediate.to_binary())
+
+print("> Terminal state : CX-gate")
+terminal = intermediate @ (gate.CX, 2)
+print(terminal.to_binary())
+# ground @ (gate.X, 3)
+
+# start = time.time()
+# print(ground.state[(0,)*size])
+# print(ground.to_binary())
+# end = time.time()
+# print(f"> Execution time : {end-start:.8f}sec")
+
+
+def timing():
+    import pickle
+    import numpy as np
+    import matplotlib.pyplot as plt
+
+    times = []
+    curve = []
+    mx = 100
+    for s in [100*i for i in range(1,mx)]:
+        print(s)
+        ground = Qbit(s)
+        start = time.time()
+        print(ground.state[(0,)*s])
+        end = time.time()
+        print(f"> Execution time : {end-start:.8f}sec")
+        times.append(end-start)
+        # curve.append(np.exp(s*np.log(5.13302064)/10000)-(1-0.005))
+
+    b = 1-times[0]
+    a = np.log((times[-1]+b))/(mx*100)
+    for s in [100*i for i in range(1,mx)]:
+        curve.append(np.exp(a*s)-b)
+
+    print(times[-1], times[0])
+
+    with open("times.txt", "wb") as fp:
+        pickle.dump(times, fp)
+
+    plt.plot([100*i for i in range(1,mx)], curve)
+    plt.plot([100*i for i in range(1,mx)], times)
+    plt.show()
 # print(ground.state.real_parameters_number, ground.state.parameters_number)
 # import numpy as np
 # from syngular.tensor import MatrixProductState