Little change and bug in model prediction after update of weights to …

…be fixed

core/model.py

@@ -1,4 +1,5 @@
-from syngular.tensor import MatrixProductOperator, MatrixProductState
+from numpy.matrixlib.defmatrix import matrix
+from syngular.tensor import MatrixProductOperator, MatrixProductState, matrix_product_operator
 
 import numpy as np
 from opt_einsum import contract
@@ -22,8 +23,16 @@ def build(self):
                 layer.build(None)
                 layer.is_built = True
 
-    def train(self):
-        pass
+    def train(self, x, y, batchsize=32, epochs=1, verbose=1):
+        if verbose: print("[START] Training ")
+        for e in range(epochs):
+            if verbose: print(f"Epoch {str(e+1)} : ")
+            for layer in self.layers:
+                for weight in layer.trainable_tensor_weights:
+                    # print(weight)
+                    weight["weight"] += MatrixProductOperator.random((2,2), (2,2), (4,))
+                    print(weight["weight"])
+        if verbose: print("[END] Training ")
 
     def draw(self):
         repr = ''
@@ -62,17 +71,19 @@ def draw(self):
         return repr
 
 
-    def add_weight(self, input_shape, output_shape, bond, name=None, initializer="normal"):
+    def add_weight(self, input_shape, output_shape, bond_shape, name=None, initializer="normal"):
         if name == None:
             name = f'weight_{np.random.randint(0,999999)}'
 
-        if initializer == "normal":
-            weight = np.random.normal(size=(*self._input_shape, *self._output_shape))
-        else:
-            weight = np.zeros(shape=(*self._input_shape, *self._output_shape))
+        # if initializer == "normal":
+        #     weight = np.random.normal(size=(*self._input_shape, *self._output_shape))
+        # else:
+        #     weight = np.zeros(shape=(*self._input_shape, *self._output_shape))
+
+        # matrix_product_weight = MatrixProductOperator(weight, bond_shape=bond)
+        # matrix_product_weight.decompose()
 
-        matrix_product_weight = MatrixProductOperator(weight, bond_shape=bond)
-        matrix_product_weight.decompose()
+        matrix_product_weight = MatrixProductOperator.random(input_shape, output_shape, bond_shape)
 
         self.trainable_tensor_weights.append({'name': name, 'weight': matrix_product_weight})
 
@@ -122,18 +133,23 @@ def __init__(self,
     def build(self, input_shape):
         # self.add_bias(self._output_shape, name="bias", initializer="normal")
         # print(self._input_shape, self._output_shape)
-        self.add_weight(self._input_shape, self._output_shape, bond=self._bond_shape, name="weight", initializer="normal")
+        if self.weights_initializer == "normal":
+            self.add_weight(self._input_shape, self._output_shape, bond_shape=self._bond_shape, name="weight", initializer="normal")
+        else:
+            self.trainable_tensor_weights.append({'name': '', 'weight': self.weights_initializer})
 
     def __call__(self, inputs):
         super(Linear, self).__call__(inputs)
-
+        print("okokokok")
         weight = self.trainable_tensor_weights[0]["weight"]
 
         print("Weight", weight)
         print("input", inputs)
-        print("contract", MatrixProductOperator.contract(weight, inputs))
-
-        return MatrixProductOperator.contract(weight, inputs)
+        print([site.shape for site in weight.sites], weight.bond_shape)
+        print([site.shape for site in inputs.sites], inputs.bond_shape)
+        print("contract", weight @ inputs)
+        # print(inputs, weight, weight @ inputs)
+        return weight @ inputs
 
 class Output(Layer):
 

quantum/gate.py

@@ -31,6 +31,11 @@
     [0., np.exp(1.j*np.pi/4)]
 ])
 
+S = np.array([
+    [1., 0.],
+    [0., np.exp(1.j*np.pi/2)]
+])
+
 CX = np.array([
     [1., 0., 0., 0.],
     [0., 1., 0., 0.],

quantum/qbit.py

@@ -34,13 +34,29 @@ def __matmul__(self, operator: Union[Tuple, Union[List, MatrixProductOperator]])
             if len(operator) == 2:
                 return Qbit.from_mps(self.state.apply(*operator))
             else:
-                qbit = self.swap_in(operator[1], operator[2]-2)
-                print('->', qbit.to_binary())
-                op = (operator[0], operator[2]-1)
-                print(op[1])
+                # print('->', self.to_binary())
+
+                min_index = min(operator[1], operator[2])
+                max_index = max(operator[1], operator[2])
+
+                qbit = self
+                if abs(min_index - max_index) != 1:
+                    qbit = self.swap_in(min_index,max_index-1)
+                    # print("swap in OKAY")
+                    # print('->', qbit.to_binary())
+                    # print("idx", op[1])
+                    print(min_index, operator[1], operator[2])
+                    if min_index != operator[1]:
+                        qbit @= (gate.SWAP, max_index-1)
+
+                op = (operator[0], max_index-1)
                 qbit = Qbit.from_mps(qbit.state.apply(*op))
-                print('->', qbit.to_binary())
-                qbit = qbit.swap_out(operator[1], operator[2]-2)
+                # print("gate applied OKAY")
+                # print('->', qbit.to_binary())
+                if abs(min_index - max_index) != 1:
+                    qbit = qbit.swap_out(min_index, max_index-1)
+                # print("swap out OKAY")
+                # print('->', qbit.to_binary())
                 return qbit
         else:
             operator = MatrixProductOperator(operator, bond_shape=()).decompose()
@@ -76,10 +92,12 @@ def swap_out(self, idx1, idx2):
         _idx1 = min(idx1, idx2)
         _idx2 = max(idx1, idx2)
 
+        # print(list(range(_idx2, _idx1-1, -1)))
+
         # self @= (gate.SWAP, _idx2)
         # print("ok", self.to_binary())
         # print(list(range(_idx2, _idx1-1, -1)))
-        for i in range(_idx2-2, _idx1-1, -1):
+        for i in range(_idx2, _idx1-1, -1):
             # print(i)
             # print(self.to_binary())
             self @= (gate.SWAP, i)
@@ -116,4 +134,16 @@ def from_mps(mps):
         # print(mps.sites_number)
         qbit = Qbit(mps.sites_number, init=False)
         qbit.state = mps
+        return qbit
+
+    @staticmethod
+    def from_binary(bin):
+        size = len(bin)
+
+        qbit = Qbit(size)
+
+        for i in range(size): 
+            if bin[i] == '1': 
+                qbit @= (gate.X, i)
+
         return qbit

tensor/matrix_product_operator.py

@@ -15,6 +15,9 @@
 
 class MatrixProductOperator:
 
+    COMPRESS = True
+    DECOMPOSE = False
+
     def __init__(self, tensor=None, bond_shape: Union[Union[Tuple, List], np.ndarray]=(), verbose=0) -> None:
 
         self.parameters_number = 0
@@ -69,7 +72,11 @@ def __init__(self, tensor=None, bond_shape: Union[Union[Tuple, List], np.ndarray
     @returns: a float representing the sum of the two MatrixProductOperator
     """
     def __add__(self, mpo):
-        print(mpo)
+        # print(mpo)
+        min_bond = min(min(self.bond_shape), min(mpo.bond_shape))
+        print("min", min_bond, self.bond_shape, mpo.bond_shape)
+        print(min_bond)
+
         if self.decomposed and mpo.decomposed:
             sites = []
             for idx in range(self.sites_number):
@@ -96,7 +103,7 @@ def __add__(self, mpo):
 
                 sites.append(site)
 
-            return MatrixProductOperator.from_sites(sites)
+            return MatrixProductOperator.from_sites(sites) >> min_bond
         else:
             raise Exception("Both Matrix Product Operator must be in canonical form (use .decompose()")
 
@@ -113,7 +120,9 @@ def __add__(self, mpo):
     @returns: a float representing the sum of the two MatrixProductOperator
     """
     def __mul__(self, mp: Union[MatrixProductOperator, MatrixProductState]):
-        print(mp)
+        # print(mp)
+        min_bond = min(min(self.bond_shape), min(mp.bond_shape))
+        print("min", min_bond, self.bond_shape, mp.bond_shape)
         if not self.decomposed and mp.decomposed:
             raise Exception("Both Matrix Product Operator must be in canonical form (use .decompose()")
 
@@ -140,7 +149,7 @@ def __mul__(self, mp: Union[MatrixProductOperator, MatrixProductState]):
 
                 sites.append(site)
 
-            return MatrixProductOperator.from_sites(sites)
+            return MatrixProductOperator.from_sites(sites) >> min_bond
 
         elif isinstance(mp, MatrixProductState):
             pass
@@ -160,6 +169,9 @@ def __mul__(self, mp: Union[MatrixProductOperator, MatrixProductState]):
     @returns: a float representing the contraction of the two MatrixProduct
     """
     def __matmul__(self, mp: Union[MatrixProductOperator, MatrixProductState]):
+        min_bond = min(min(self.bond_shape), min(mp.bond_shape))
+        max_bond = max(max(self.bond_shape), max(mp.bond_shape))
+        print("min", min_bond)
         sites = []
 
         if isinstance(mp, MatrixProductState):
@@ -172,7 +184,7 @@ def __matmul__(self, mp: Union[MatrixProductOperator, MatrixProductState]):
                     self.shape[idx][3]*mp.shape[idx][2]
                 ))
                 sites.append(site)
-            return MatrixProductState.from_sites(sites)
+            return MatrixProductState.from_sites(sites) >> min_bond
         elif isinstance(mp, MatrixProductOperator):
             for idx in range(self.sites_number):
                 site = contract(self.sites[idx], [1,2,3,4], mp.sites[idx], [5,3,6,7], [1,5,2,6,4,7])
@@ -184,7 +196,7 @@ def __matmul__(self, mp: Union[MatrixProductOperator, MatrixProductState]):
                     self.shape[idx][3]*mp.shape[idx][3]
                 ))
                 sites.append(site)
-            return MatrixProductOperator.from_sites(sites)
+            return MatrixProductOperator.from_sites(sites) >> min_bond
 
     def __mod__(self, mode):
         if mode == 'L' or mode == 'left' or mode == 'l' or mode == 'left-canonical':
@@ -366,6 +378,13 @@ def compress(self, dim, mode='left', strict=False):
         n = self.sites_number-1
         parameters_number = 0
 
+        if dim >= min(self.bond_shape):
+            return self
+
+        print(self.bond_shape)
+        self.bond_shape = (dim,) * (self.sites_number-1)
+        print(self.bond_shape)
+
         if not strict:
             if mode == 'left':
                 if self.orthonormalized != 'left': self.left_orthonormalization()
@@ -381,8 +400,8 @@ def compress(self, dim, mode='left', strict=False):
 
                     W = S @ R
 
-                    print(Q.shape, L.shape, S.T.shape, R.shape)
-                    print(W.shape)
+                    # print(Q.shape, L.shape, S.T.shape, R.shape)
+                    # print(W.shape)
 
                     l1 = list(self.shape[k])
                     l2 = list(self.shape[k+1])
@@ -433,7 +452,7 @@ def compress(self, dim, mode='left', strict=False):
         else:
             that = self.copy()
 
-            print(that)
+            # print(that)
 
             if mode == 'left':
                 # if self.orthonormalized != 'left': self.left_orthonormalization()
@@ -448,8 +467,8 @@ def compress(self, dim, mode='left', strict=False):
 
                     W = S @ R
 
-                    print(Q.shape, L.shape, S.T.shape, R.shape)
-                    print(W.shape)
+                    # print(Q.shape, L.shape, S.T.shape, R.shape)
+                    # print(W.shape)
 
                     l1 = list(that.shape[k])
                     l2 = list(that.shape[k+1])
@@ -534,7 +553,7 @@ def apply(self, operator, indices):
                 # struct += output_shape
                 # print(struct)
                 T = contract(*struct)
-                print(T.shape)
+                # print(T.shape)
 
         else:
             raise Exception("MatrixProductState not decomposed")
@@ -555,7 +574,7 @@ def retrieve(self, input_indices, output_indices):
 
     def left_orthonormalization(self, bond_shape=()):
         for k in range(self.sites_number-1):
-            print(k, k+1)
+            # print(k, k+1)
             L = self.left_site_matricization(k)
             R = self.right_site_matricization(k+1)
 

tensor/matrix_product_state.py

@@ -300,6 +300,10 @@ def decompose(self, mode="left"):
     def compress(self, dim: int, mode="left", strict=False):
         n = self.sites_number-1
         parameters_number = 0
+
+        print(dim, min(self.bond_shape), self.bond_shape)
+        if dim >= min(self.bond_shape):
+            return self
 
         if not strict:
             if mode == 'left':
@@ -408,6 +412,7 @@ def apply(self, operator, index, strict=True):
         n = that.sites_number
 
         # index = n-1-index
+        # print(">",n, m, index)
 
         struct = []
         struct += [operator, [*list(range(3*n+index, 3*n+m+index)), *list(range(2*n+index, 2*n+m+index))]]
@@ -418,6 +423,7 @@ def apply(self, operator, index, strict=True):
 
         T = contract(*struct)
         # print('T', T.shape)
+        # print('op', operator.shape)
 
         # print(index, m, m+index-1, list(range(index, m+index-1)))
         for k in range(index, m+index-1):
@@ -433,15 +439,16 @@ def apply(self, operator, index, strict=True):
             R = R[:rank, :]
             # print(k, m+k, 'L', L.shape, 'Q', Q.shape, 'R', R.shape)
             # print(that.sites[k].shape)
-            # print(m+k, operator.shape[m+k], rank)
+            # print(m+k-1, operator.shape, rank)
             that.sites[k] = that.tensoricization(Q, k)
-            T = R.reshape((rank, operator.shape[m+k-1], -1))
+            T = R.reshape((rank, -1, that.shape[k+1][2]))
             # print('Tr', T.shape)
 
         # print(that.sites[m+index-1])
         # print('---')
         # print(T)
         # print(that.sites[m+index-1].shape)
+        # print(m+index-1)
         that.sites[m+index-1] = that.tensoricization(T, m+index-1)
         # print(that.sites[m+index-1].shape)
 
@@ -532,3 +539,6 @@ def left_orthogonality(self, index):
     def right_orthogonality(self, index):
         R = self.right_matricization(self.sites[index], index)
         return R @ R.T
+
+    def grad(self, index):
+        return self.sites[:index]+self.sites[index+2:]