diff --git a/core/model.py b/core/model.py new file mode 100644 index 0000000..743e401 --- /dev/null +++ b/core/model.py @@ -0,0 +1,148 @@ +from numpy.linalg.linalg import matrix_rank +from syngular.tensor.tensor_train import MatrixProductOperator, MatrixProductState + +import numpy as np +from opt_einsum import contract + +class Model: + + def __init__(self, layers): + self.layers = layers + + def predict(self, inputs): + values = inputs + for layer in self.layers: + values = layer(values) + # print("Values shape", values._shape) + + return values + + def build(self): + for layer in self.layers: + if not layer.is_built: + layer.build(None) + layer.is_built = True + + def train(self): + pass + + def draw(self): + repr = '' + for layer in self.layers: + repr += layer.draw() + return repr + + +class Layer: + + def __init__(self): + + self.trainable_tensor_weights = [] + self.trainable_tensor_bias = [] + + self.is_built = False + + def __call__(self, inputs): + input_shape = inputs.shape + if not self.is_built: + self.build(input_shape) + self.is_built = True + else: + print("Built") + + def build(self, input_shape): + pass + + def draw(self): + repr = '' + for weight in self.trainable_tensor_weights: + mp = weight["weight"] + repr += "\t"+"| " * mp.sites_number + "\n" + repr += "\t"+("O---" * (mp.sites_number-1)) + "O" + "\n" + repr += "\t"+"| " * mp.sites_number + "\n" + return repr + + + def add_weight(self, input_shape, output_shape, bond, 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)) + + matrix_product_weight = MatrixProductOperator(weight, bond_shape=bond) + matrix_product_weight.decompose(mode='left') + + self.trainable_tensor_weights.append({'name': name, 'weight': matrix_product_weight}) + + + + def add_bias(self, size, name=None, initializer="normal"): + if name == None: + name = f'bias_{np.random.randint(0,999999)}' + + if initializer == "normal": + bias = np.random.normal(size=size) + else: + bias = np.zeros(shape=size) + + self.trainable_tensor_bias.append({name: name, bias: bias}) + +class Linear(Layer): + + def __init__(self, + input_units, output_units, + core=1, bond=None, + bias_initializer="normal", + weights_initializer="normal", + activation="relu"): + + super(Linear, self).__init__() + + self.input_units = input_units + self.output_units = output_units + + self.core = core + self.bond_dimension = bond + + self.input_core_dim = round(self.input_units**(1/self.core)) + self.output_core_dim = round(self.output_units**(1/self.core)) + + self._input_shape = (self.input_core_dim,) * self.core + self._output_shape = (self.output_core_dim,) * self.core + self._bond_shape = (self.bond_dimension,) * (self.core-1) + + self.bias_initializer = bias_initializer + self.weights_initializer = weights_initializer + + self.activation = activation + + 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="bias", initializer="normal") + + def __call__(self, inputs): + super(Linear, self).__call__(inputs) + + weight = self.trainable_tensor_weights[0]["weight"] + + print("Weight", weight) + print("input", inputs) + print("contract", MatrixProductOperator.contract(weight, inputs)) + + return MatrixProductOperator.contract(weight, inputs) + +class Output(Layer): + + def __init__(self, output_shape): + super(Output, self).__init__() + + self.output_shape = output_shape + + def __call__(self, inputs): + # print(">", inputs) + + return inputs#.reshape(self.output_shape) \ No newline at end of file diff --git a/documentation/matrix_product/MatrixProduct.md b/documentation/matrix_product/MatrixProduct.md new file mode 100644 index 0000000..3a11e61 --- /dev/null +++ b/documentation/matrix_product/MatrixProduct.md @@ -0,0 +1,46 @@ +# Matrix Product + + +Importation +```python +from syngular.tensor.tensor_train import MatrixProductState +from syngular.tensor.tensor_train import MatrixProductOperator +``` + +### Addition + +```python +tensor = np.arange(0,16).reshape((2,2,2,2)) + +X = MatrixProductOperator(tensor, bond_shape=(4,)) +Y = MatrixProductOperator(tensor, bond_shape=(4,)) + +Z = X + Y +``` + +### Multiplication + +```python +tensor = np.arange(0,16).reshape((2,2,2,2)) +tensor2 = np.arange(0,4).reshape((2,2)) + +X = MatrixProductOperator(tensor, bond_shape=(4,)) +Y = MatrixProductOperator(tensor, bond_shape=(4,)) +Z = MatrixProductOperator(tensor, bond_shape=(4,)) +S = MatrixProductState(tensor2, bond_shape=(4,)) + +W = S @ ((X @ Y) + (Y @ X)) @ S +``` + +![Diagram](tensor.jpg) + +### Compression + + +```python +tensor = np.arange(0,16).reshape((2,2,2,2)) + +X = MatrixProductOperator(tensor, bond_shape=(4,)) + +W = X >> 2 +``` diff --git a/documentation/matrix_product/tensor.jpg b/documentation/matrix_product/tensor.jpg new file mode 100644 index 0000000..cc37915 Binary files /dev/null and b/documentation/matrix_product/tensor.jpg differ diff --git a/examples/resources/images/paysage.jpg b/examples/resources/images/paysage.jpg new file mode 100644 index 0000000..309c8fa Binary files /dev/null and b/examples/resources/images/paysage.jpg differ diff --git a/examples/resources/images/paysage2.jpg b/examples/resources/images/paysage2.jpg new file mode 100644 index 0000000..e611aec Binary files /dev/null and b/examples/resources/images/paysage2.jpg differ diff --git a/examples/resources/images/test.png b/examples/resources/images/test.png new file mode 100644 index 0000000..7a9977a Binary files /dev/null and b/examples/resources/images/test.png differ diff --git a/examples/resources/images/test2.png b/examples/resources/images/test2.png new file mode 100644 index 0000000..0af5bfe Binary files /dev/null and b/examples/resources/images/test2.png differ diff --git a/layers/TensorDense.py b/layers/TensorDense.py index 9e16cd4..c11ce24 100644 --- a/layers/TensorDense.py +++ b/layers/TensorDense.py @@ -9,6 +9,8 @@ import tensornetwork as tn import numpy as np +from opt_einsum import contract + def unfold_dim(shape): return reduce(lambda x, y: x*y, shape) @@ -40,9 +42,9 @@ def build(self, tt_input_shape): self.tt_input_shape = tuple([roots] * self.cores_number) - print("reshaped", tt_input_shape) + #print("reshaped", tt_input_shape) - print("Input shape", self.tt_input_shape) + #print("Input shape", self.tt_input_shape) self.bias = tf.Variable(tf.zeros(shape=self.tt_output_shape), name="bias", trainable=True) @@ -76,31 +78,49 @@ def build(self, tt_input_shape): def call(self, inputs): def process(input, bias): - input_reshaped = tf.reshape(input,self.tt_input_shape) + # input_reshaped = tf.reshape(input,self.tt_input_shape) - last_idx = self.cores_number-1 - einsum_structure = [] + # last_idx = self.cores_number-1 + # einsum_structure = [] - cores = [tn.Node(core, backend="tensorflow").tensor for core in self.cores] - x = tn.Node(input_reshaped, backend="tensorflow") + # cores = [tn.Node(core, backend="tensorflow").tensor for core in self.cores] + # x = tn.Node(input_reshaped, backend="tensorflow") - einsum_structure = [] - einsum_structure.append(list(range(1, self.cores_number+1))) - einsum_structure.append([1, -(self.cores_number+1), 2*self.cores_number+1]) + # einsum_structure = [] + # einsum_structure.append(list(range(1, self.cores_number+1))) + # einsum_structure.append([1, -(self.cores_number+1), 2*self.cores_number+1]) - for idx in range(1, last_idx): - einsum_structure.append([idx+1, -(self.cores_number+idx+1), 2*self.cores_number+idx+1]) + # for idx in range(1, last_idx): + # einsum_structure.append([idx+1, -(self.cores_number+idx+1), 2*self.cores_number+idx+1]) - einsum_structure.append([last_idx+1, -(self.cores_number+last_idx+1), 2*self.cores_number+last_idx-1+1]) + # einsum_structure.append([last_idx+1, -(self.cores_number+last_idx+1), 2*self.cores_number+last_idx-1+1]) - print(einsum_structure) + # #print(einsum_structure) + + # result = tn.ncon( + # tensors = [x.tensor] + cores, + # network_structure = einsum_structure, + # backend = "tensorflow" + # ) + + x = tf.reshape(input,self.tt_input_shape)#tf.reshape(input, ((self.tt_input_shape[0],)*self.cores_number)) + + + struct = [] + struct += [x, list(range(1, self.cores_number+1))] + struct += [self.cores[0], [1, self.cores_number+1, 2*self.cores_number+1]] + + for idx in range(2, self.cores_number): + struct += [self.cores[idx-1]] + struct += [[idx, self.cores_number+idx, 2*self.cores_number+1, 2*self.cores_number+2]] + + struct += [self.cores[-1], [self.cores_number, 2*self.cores_number, 3*self.cores_number-1]] + struct += [list(range(self.cores_number+1, 2*self.cores_number+1))] + + + result = contract(*struct) - result = tn.ncon( - tensors = [x.tensor] + cores, - network_structure = einsum_structure, - backend = "tensorflow" - ) # print(type(self.cores[0].numpy())) # cores = list(map(lambda core: tf.convert_to_tensor(core), self.cores)) diff --git a/layers/TensorRNN.py b/layers/TensorRNN.py index 0f92ec2..92a1da1 100644 --- a/layers/TensorRNN.py +++ b/layers/TensorRNN.py @@ -84,3 +84,12 @@ def process(input): result = tf.vectorized_map(lambda vec: process(vec, self.bias), inputs) return self.activation(tf.reshape(result, (-1, self.tt_output_shape_unfold))) + +if __name__ == "__main__": + + from keras.layers import SimpleRNN + from keras import Sequential + + model = Sequential([ + SimpleRNN() + ]) \ No newline at end of file diff --git a/tensor/deep.py b/tensor/deep.py new file mode 100644 index 0000000..e045b94 --- /dev/null +++ b/tensor/deep.py @@ -0,0 +1,129 @@ +import itertools +from PIL.Image import new +from numpy.core.fromnumeric import size +from syngular.tensor.tensor_train import MatrixProductOperator, MatrixProductState + +import numpy as np + + +class Model: + + def __init__(self, layers): + self.layers = layers + + def predict(self, inputs): + values = inputs + for layer in self.layers: + values = layer(values) + + return values + + def train(self): + pass + +class Layer: + + def __init__(self) -> None: + self.trainable_weights = [] + + self.build() + + def build(self): + pass + + def train(self): + for trainable in self.trainable_weights: + pass + + +class Shaper(Layer): + + def __init__(self, shape, bond_shape=None) -> None: + super(Shaper, self).__init__() + + self.shape = shape + self.bond_shape = bond_shape + + def __call__(self, inputs): + print(inputs.shape) + outputs = np.reshape(inputs, newshape=self.shape, order="F") + + if self.bond_shape != None: + outputs = MatrixProductState(outputs, bond_shape=self.bond_shape) + outputs.decompose(mode='left') + + return outputs + + +class Linear(Layer): + + def __init__(self, unit_shape, bond_shape) -> None: + super(Linear, self).__init__() + + self.unit_shape + self.bond_shape + + self.weights = None + self.bias = None + + def build(self): + self.bias = np.zeros(self.bond_shape) + # self.weights = np.random.normal(size=) + + def __call__(self, inputs): + pass + +if __name__ == "__main__": + + import matplotlib.pyplot as plt + import matplotlib.image as mimage + + + tensor = np.arange(4**4).reshape((4,4,4,4)) + image = np.arange(32*32) + + img = mimage.imread('C:\\Antoine\\Coding\\Tensor Network\Project\\syngular\\examples\\resources\images\\paysage2.jpg') + img = img[:,:,0] / 255.0 + + + + + + model = Model([ + Shaper((20,20,30,20), (20,20,20,)) # 20800 / 240 000 + #Shaper((64,75,50), (50,50,)) # 193200 / 240 000 + ]) + + y = model.predict(img) + + + + + + img_reconstruct = y.reconstruct() + img_reconstruct = img_reconstruct.reshape((400, 600), order="F") + + + errors = [] + for i, j in zip(range(400), range(600)): + e = abs(img[i,j] - img_reconstruct[i,j]) + errors.append(e) + + # print(img_reconstruct[a,b,c,d,e]) + # errors1.append(error1) + # plt.plot(errors) + # + plt.rcParams['figure.dpi'] = 150 + + print("Shape", img_reconstruct.shape) + plt.subplot(2,1,1) + plt.imshow(img, cmap='gray') + plt.axis('off') + plt.subplot(2,1,2) + plt.imshow(img_reconstruct, cmap='gray') + plt.axis('off') + plt.show() + + + + diff --git a/tensor/tensor_train.py b/tensor/tensor_train.py index 949e580..2aeaf5d 100644 --- a/tensor/tensor_train.py +++ b/tensor/tensor_train.py @@ -1,6 +1,9 @@ +import gc import itertools +from typing import Tuple, List, Union + import numpy as np -from syngular.tensor.mpo import MatrixProductOperator +from scipy import linalg import matplotlib.pyplot as plt @@ -12,106 +15,457 @@ def __init__(self, tensor, bond_shape) -> None: self.bond_shape = bond_shape self.sites_number = len(self.bond_shape)+1 - self.sites = [] + self.sites = [None] * self.sites_number + + # print(self.bond_shape, self.sites_number) + + 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)] + self.shape = [(self.tensor_shape[i], self.bond_shape[i]) for i in range(self.sites_number-1)] - print(self.shape) + def decompose(self, mode="left"): - def decompose(self): + if mode == "left": + current_matrix = self.tensor + current_shape = self.shape[0][0] + current_rank = 1 - current_matrix = self.tensor - current_shape = self.shape[0][0] - current_rank = 1 + # print(self.shape) + for k in range(self.sites_number-1): + current_shape = self.shape[k][0] + + unfolding_matrix = np.reshape(current_matrix, newshape=(current_shape*current_rank, -1)) + rank = self.shape[k][1] + + Q, R = np.linalg.qr(unfolding_matrix, mode="complete") - for k in range(self.sites_number-1): - - unfolding_matrix = np.reshape(current_matrix, newshape=(current_shape*current_rank, -1)) - rank = self.shape[k][1] - - Q, R = np.linalg.qr(unfolding_matrix, mode="complete") + Q = Q[:,:rank] + Q = np.reshape(Q, newshape=(current_rank, current_shape, -1)) + R = R[:rank, :] + + print(f"Core {k} with {current_rank} , {current_shape}") - Q = Q[:,:rank] - Q = np.reshape(Q, newshape=(current_rank, current_shape, -1)) - R = R[:rank, :] + self.sites[k] = Q - self.sites.append(Q) + current_rank = rank + current_matrix = R - current_rank = rank - current_matrix = R + current_matrix = current_matrix[:, :, np.newaxis] + self.sites[-1] = current_matrix - current_matrix = current_matrix[:, :, np.newaxis] - self.sites.append(current_matrix) + # print([y.shape for y in self.sites]) + + elif mode == "right": + current_matrix = self.tensor.T + current_shape = self.shape[0][0] + current_rank = 1 + + self.decompose(mode="left") + + # print("sites", self.shape) + return self def retrieve(self, indices): einsum_structure = [] for idx in range(self.sites_number): + # print(indices[idx], self.sites[idx].shape) einsum_structure.append(self.sites[idx][:, indices[idx], :]) - einsum_structure.append([idx, idx+1]) + einsum_structure.append([idx, Ellipsis, idx+1]) return np.einsum(*einsum_structure) + + def update(self, sites, value): + + print(self.sites[sites[0]].shape) + print(self.sites[sites[1]].shape) + + merge_site = np.einsum(self.sites[sites[0]], [1,2,3], self.sites[sites[1]], [3,4,5], [1,2,4,5]) + + merge_site = np.einsum(merge_site, [1,2,3,4], value, [5,2,3,6], [1,5,6,4]) + + print("shape") + print(merge_site.shape) + + rank_left, rank_1, rank_2, rank_right = merge_site.shape + + unfolding_matrix = np.reshape(merge_site, newshape=(rank_left*rank_1, -1)) + + Q, R = np.linalg.qr(unfolding_matrix, mode="complete") + + Q = Q[:,:rank_left] + Q = np.reshape(Q, newshape=(rank_left, rank_1, -1)) + R = R[:rank_left, :] + R = np.reshape(R, newshape=(-1, rank_2, rank_right)) + + # print(Q.shape, R.shape) + + self.sites[sites[0]] = Q + self.sites[sites[1]] = R + + + def __repr__(self) -> str: + repr = " \n> Sites shape" + str(self.shape) + "\n" + repr += "\t"+"| " * self.sites_number + "\n" + repr += "\t"+("O---" * (self.sites_number-1)) + "O" + "\n" + return repr + + + + + + + + + + + + + + + + + + + + + + + + class MatrixProductOperator: - def __init__(self, tensor, bond_shape) -> None: - self.tensor = tensor - self.tensor_shape = tensor.shape - self.bond_shape = bond_shape + def __init__(self, tensor=(), bond_shape: Union[Union[Tuple, List], np.ndarray]=(), verbose=0) -> None: + if tensor != (): + self.tensor = tensor + self.tensor_shape = tensor.shape + self.bond_shape = tuple(bond_shape) - self.sites_number = len(self.bond_shape)+1 - self.sites = [] + self.sites_number = len(self.bond_shape)+1 + self.sites = [None] * self.sites_number - 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.bond_shape = (1,) - print(self._shape) + self.input_shape = tuple(self.tensor_shape[:self.sites_number]) + self.output_shape = tuple(self.tensor_shape[self.sites_number:]) - self.tensor = np.transpose(self.tensor, axes=sum(list(zip(range(self.sites_number), range(self.sites_number, 2*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") - def decompose(self): + if len(self.input_shape) != len(self.bond_shape)+1: + raise Exception("dimensions of bond indices do not match input dimension - 1") - current_matrix = self.tensor - current_input_shape = self._shape[0][1] - current_output_shape = self._shape[0][2] - current_rank = 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)] - for k in range(self.sites_number-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) + ] + + self.tensor = np.transpose(self.tensor, axes=sum(list(zip(range(self.sites_number), range(self.sites_number, 2*self.sites_number))), ())) + + self.verbose = verbose + self.decomposed = False + + def __add__(self, mpo): + print(mpo) + if self.decomposed and mpo.decomposed: + sites = [] + for idx in range(self.sites_number): + wing_left = idx == 0 + wing_right = idx == self.sites_number-1 + + site = np.zeros(shape=( + self.shape[idx][0] + mpo.shape[idx][0] if not wing_left else self.shape[idx][0], + self.shape[idx][1], + self.shape[idx][2], + self.shape[idx][3] + mpo.shape[idx][3] if not wing_right else self.shape[idx][3] + )) + + for inp in range(self.shape[idx][1]): + for out in range(self.shape[idx][2]): + left_matrix = self.sites[idx][:, inp, out, :] + right_matrix = mpo.sites[idx][:, inp, out, :] + if wing_left: + site[:, inp, out, :] = np.block([left_matrix, right_matrix]) + elif wing_right: + site[:, inp, out, :] = np.block([left_matrix.T, right_matrix.T]).T + else: + site[:, inp, out, :] = linalg.block_diag(left_matrix, right_matrix) + + sites.append(site) - unfolding_matrix = np.reshape(current_matrix, newshape=(current_input_shape*current_output_shape*current_rank, -1)) - rank_right = self._shape[k][3] + return MatrixProductOperator.from_sites(sites) + else: + raise Exception("Both Matrix Product Operator must be in canonical form (use .decompose()") + + def __rshift__(self, dim): + if isinstance(dim, int): + return self.compress(dim) + else: + raise Exception("dimension should be an integer") + + def __getitem__(self, key): + key_inp = key[0] + key_out = key[1] + + if len(key_inp) != self.sites_number: + raise Exception("input indices do not match the number of sites") + if len(key_out) != self.sites_number: + raise Exception("output indices do not match the number of sites") + + return self.retrieve(key_inp, key_out) + + def __repr__(self) -> str: + repr = " \n> Sites shape" + str(self.shape) + "\n" + repr += "\t"+"| " * self.sites_number + "\n" + repr += "\t"+("O---" * (self.sites_number-1)) + "O" + "\n" + repr += "\t"+"| " * self.sites_number + "\n" + return repr + + + @staticmethod + def random(input_shape, output_shape, bond_shape): + tensor = np.random.normal(size=(*input_shape, *output_shape)) + return MatrixProductOperator(tensor, bond_shape=bond_shape).decompose() + + @staticmethod + def from_sites(sites): + mp = MatrixProductOperator() + mp.sites = sites + mp.sites_number = len(sites) + + mp.decomposed = True + + mp.input_shape, mp.output_shape, mp.bond_shape = (), (), () + mp.shape = [None] * mp.sites_number + + for idx in range(mp.sites_number): + site = sites[idx] + shape = site.shape + mp.shape[idx] = shape + mp.input_shape += (shape[1],) + mp.output_shape += (shape[2],) + if idx != mp.sites_number-1: mp.bond_shape += (shape[3],) + + return mp + + @staticmethod + def empty(): + return MatrixProductOperator() + + def to_tensor(self): + tensor = np.zeros(shape=(*self.input_shape, *self.output_shape)) + + range_inp = [range(inp) for inp in self.input_shape] + range_out = [range(out) for out in self.output_shape] + + for inp in itertools.product(*range_inp): + for out in itertools.product(*range_out): + tensor[(*inp, *out)] = self[inp, out] + + return tensor + + def decompose(self, mode="left"): + if not self.decomposed: + if mode == "left": + current_matrix = self.tensor + current_input_shape = self.shape[0][1] + current_output_shape = self.shape[0][2] + current_rank = 1 + + for k in range(self.sites_number-1): + + unfolding_matrix = np.reshape(current_matrix, newshape=(current_input_shape*current_output_shape*current_rank, -1)) + rank_right = self.shape[k][3] + + Q, R = np.linalg.qr(unfolding_matrix, mode="complete") + + Q = Q[:,:rank_right] + Q = np.reshape(Q, newshape=(current_rank, current_input_shape, current_output_shape, rank_right)) + R = R[:rank_right, :] + + self.sites[k] = Q + + current_input_shape = self.shape[k][1] + current_output_shape = self.shape[k][2] + current_rank = rank_right + current_matrix = R + + current_matrix = np.reshape(current_matrix, newshape=(-1, current_input_shape, current_output_shape))[:, :, :, np.newaxis] + if self.verbose == 1: print(current_matrix.shape) + self.sites[-1] = current_matrix + + elif mode == "right": + if self.verbose == 1: print(self.sites_number) + axes = list(range(self.sites_number-1, -1, -1)) + list(range(2*self.sites_number-1, self.sites_number-1, -1)) + if self.verbose == 1: print("axes", axes) + current_matrix = np.transpose(self.tensor, axes=axes) + + self.decompose(mode="left") - Q, R = np.linalg.qr(unfolding_matrix, mode="complete") + self.sites = self.sites[::-1] + + self.decomposed = True - Q = Q[:,:rank_right] - Q = np.reshape(Q, newshape=(current_rank, current_input_shape, current_output_shape, rank_right)) - R = R[:rank_right, :] + del self.tensor + gc.collect() - self.sites.append(Q) + return self - current_input_shape = self._shape[k][1] - current_output_shape = self._shape[k][2] - current_rank = rank_right - current_matrix = R + def compress(self, dim): + for k in range(self.sites_number-1): + pass + # Q, R = np.linalg.qr(unfolding_matrix, mode="complete") - current_matrix = np.reshape(current_matrix, newshape=(-1, current_input_shape, current_output_shape))[:, :, :, np.newaxis] - print(current_matrix.shape) - self.sites.append(current_matrix) def retrieve(self, input_indices, output_indices): einsum_structure = [] + for idx in range(self.sites_number): einsum_structure.append(self.sites[idx][:, input_indices[idx], output_indices[idx], :]) einsum_structure.append([idx, idx+1]) return np.einsum(*einsum_structure) + # def copy(self): + # new_mpo = MatrixProductOperator(self.tensor, self.bond_shape) + # new_mpo.sites = self.sites + # return new_mpo + + def contract(self, mp, indices=None): + if isinstance(mp, MatrixProductState): + for idx in range(self.sites_number): + self.sites[idx] = np.einsum(self.sites[idx], [1,2,3,4], mp.sites[idx], [5,6,2,7], [5,3,6,7]) + elif isinstance(mp, MatrixProductOperator): + print("ok") + + ''' + Contraction through the input indices of the current MPO and the output indices of the passed MPO + > (0, ..., N) + + | | | | + O---O---O---O + | | | | + O---O---O---O + | | | | + ''' + + for idx in range(self.sites_number): + self.sites[idx] = np.einsum(self.sites[idx], [1,2,3,4], mp.sites[idx], [5,6,2,7], [5,3,6,7]) + return self + + + def contract(mp1, mp2, mode="left"): + sites = [None] * mp1.sites_number + if isinstance(mp2, MatrixProductState): + + if mode == 'left': + struct = [] + n = mp2.sites_number + + for idx in range(n): + struct += [mp2.sites[idx], [idx+1, n+idx+2, idx+2]] + struct += [mp1.sites[idx], [2*n+2+idx, n+idx+2, 3*n+idx+3, 2*n+2+idx+1]] + struct += [[]] + # struct += [mp2, list(range(1, self.cores_number+1))] + # struct += [self.cores[0], [1, self.cores_number+1, 2*self.cores_number+1]] + + # for idx in range(2, self.cores_number): + # struct += [self.cores[idx-1]] + # struct += [[idx, self.cores_number+idx, 2*self.cores_number+1, 2*self.cores_number+2]] + + # struct += [self.cores[-1], [self.cores_number, 2*self.cores_number, 3*self.cores_number-1]] + # struct += [list(range(self.cores_number+1, 2*self.cores_number+1))] + + + result = np.einsum(*struct) + + for idx in range(mp1.sites_number): + # print(mp1._shape[idx], mp2._shape[idx]) + print(mp1.sites[idx].shape, mp2.sites[idx].shape) + sites[idx] = np.einsum(mp1.sites[idx], [1,2,3,4], mp2.sites[idx], [5,2,6], [5,3,6]) + # p rint(sites[idx].shape) + + elif isinstance(mp2, MatrixProductOperator): + ''' + Contraction through the input indices of the current MPO and the output indices of the passed MPO + > (0, ..., N) + + | | | | + O---O---O---O + | | | | + O---O---O---O + | | | | + ''' + + if mode == 'left': + struct = [] + n = mp2.sites_number + + for idx in range(n): + struct += [[idx+1, n+idx+2, 2*n+idx+2, idx+2]] + struct += [[3*n+2+idx, 2*n+idx+2, 4*n+idx+3, 3*n+2+idx+1]] + print("STRUCT", struct) + result = np.einsum(*struct) + + # for idx in range(mp1.sites_number): + # sites[idx] = np.einsum(mp1.sites[idx], [1,2,3,4], mp2.sites[idx], [5,6,2,7]) + else: + print("Type", type(mp2)) + # print("Sites", sites) + print("Result", result) + mp = mp2.copy() + mp.sites = sites + return mp + + + + + + + + + + + + + + + + + +class MatrixProductDensityOperator: + + def __init__(self, tensor, bond_shape) -> None: + self.tensor = tensor + self.bond_shape = bond_shape + + + self.state = MatrixProductOperator(tensor, self.bond_shape) + self.state.decompose() + + self.state_conjugate = MatrixProductOperator(tensor.T, self.bond_shape) + self.state_conjugate.decompose() + + contract_indices = [self.state._shape_indices[i][2] for i in range(self.state.sites_number)] + contract_conj_indices = [self.state_conjugate._shape_indices[i][2] for i in range(self.state.sites_number)] + output_shape = [] + + self.state.copy().contract(self.state_conjugate) + # self.mpdo = #np.einsum(self.state, contract_indices, self.state_conjugate, contract_conj_indices, output_shape) + + if __name__ == "__main__": - t = np.arange(5**7).reshape((5,5,5,5,5,5,5)) - mps = MatrixProductState(t, bond_shape=(2,2,2,2,2,2)) - mps.decompose() + # t = np.arange(5**7).reshape((5,5,5,5,5,5,5)) + # mps = MatrixProductState(t, bond_shape=(2,2,2,2,2,2)) + # mps.decompose() # errors = [] @@ -131,15 +485,18 @@ def retrieve(self, input_indices, output_indices): print(tensor) mps = MatrixProductState(tensor, bond_shape=(2,2)) - mps.decompose() + mps.decompose(mode="left") + + # update = np.arange(3**4).reshape((3,3,3,3)) + # print(mps.update((1,2), update)) - print("MPS ") - print([s.shape for s in mps.sites]) - print(mps.sites) + # print("MPS ") + # print([s.shape for s in mps.sites]) + # print(mps.sites) - print(tensor[1,1,1]) - for (i,j,k) in itertools.product(range(3), range(3), range(3)): - print(mps.retrieve((i,j,k))) + # print(tensor[1,1,1]) + # for (i,j,k) in itertools.product(range(3), range(3), range(3)): + # print(mps.retrieve((i,j,k))) # errors = [] @@ -153,22 +510,29 @@ def retrieve(self, input_indices, output_indices): operator_tensor = np.arange(3**4 * 4**4).reshape((3,3,3,3,4,4,4,4)) # print(operator_tensor) - mpo = MatrixProductOperator(operator_tensor, (2,2,2)) - mpo.decompose() - print("retrieve >") - print(operator_tensor[1,1,1,1,1,1,1,1]) - print(mpo.retrieve((1,1,1,1), (1,1,1,1))) - print(operator_tensor[0,0,0,0,0,0,0,0]) - print(mpo.retrieve((0,0,0,0), (0,0,0,0))) - print(operator_tensor[1,1,1,0,1,1,1,1]) - print(mpo.retrieve((1,1,1,0), (1,1,1,1))) + mpo = MatrixProductOperator(operator_tensor, (4,4,4)) + mpo.decompose(mode="left") + # mpo2 = mpo.copy() + # mpo2.decompose(mode="left") + # print("retrieve >") + # print(operator_tensor[1,1,1,1,1,1,1,1]) + # print(mpo.retrieve((1,1,1,1), (1,1,1,1))) + # print(operator_tensor[0,0,0,0,0,0,0,0]) + # print(mpo.retrieve((0,0,0,0), (0,0,0,0))) + # print(operator_tensor[1,1,1,0,1,1,1,1]) + # print(mpo.retrieve((1,1,1,0), (1,1,1,1))) + # print(mpo) - print([s.shape for s in mpo.sites]) + # print([s.shape for s in mpo.sites]) - errors = [] + errors1, errors2 = [], [] for (a,b,c,d ,e,f,g,h) in itertools.product(range(3), range(3), range(3), range(3) , range(4), range(4), range(4), range(4)): - error = abs(operator_tensor[a,b,c,d,e,f,g,h] - mpo.retrieve((a,b,c,d), (e,f,g,h))[0][0]) - errors.append(error) - plt.plot(errors) + error1 = abs(operator_tensor[a,b,c,d,e,f,g,h] - mpo.retrieve((a,b,c,d), (e,f,g,h))[0][0]) + errors1.append(error1) + + # error2 = abs(operator_tensor[a,b,c,d,e,f,g,h] - mpo2.retrieve((a,b,c,d), (e,f,g,h))[0][0]) + # errors2.append(error2) + plt.plot(errors1) + # plt.plot(errors2) plt.show() \ No newline at end of file diff --git a/test/model.py b/test/model.py new file mode 100644 index 0000000..0f0be1a --- /dev/null +++ b/test/model.py @@ -0,0 +1,64 @@ +import itertools +from syngular.core.model import Model, Output +from syngular.core.model import Linear +from syngular.tensor.tensor_train import MatrixProductDensityOperator, MatrixProductState, MatrixProductOperator + +import numpy as np + +X = np.arange(1,17).reshape((2,2,2,2)) +Y = np.arange(18,18+16).reshape((2,2,2,2)) +Z = X + Y +print(X) +print() +print(Y) +print() +print(Z) + + +Xmp = MatrixProductOperator(X, bond_shape=(3,)).decompose() +Ymp = MatrixProductOperator(Y, bond_shape=(3,)).decompose() +Zmp = Xmp + Ymp +print(Zmp) + + +print(X[1,1,1,1]) +print(Xmp.retrieve((1,1),(1,1))) + +print(Y[1,1,1,1]) +print(Ymp.retrieve((1,1),(1,1))) + +print(Z[1,1,1,1]) +print(Zmp.retrieve((1,1),(1,1))) + +for a,b,c,d in itertools.product(range(2), range(2), range(2), range(2)): + print(Z[a,b,c,d], Zmp.retrieve((a,b),(c,d))) + +# model = Model([ +# Linear(4,4,core=2, bond=4), +# Linear(4,1,core=2,bond=2), +# Output((1,)) +# ]) +# model.build() + +# print("---") +# print(model.draw()) +# print("---") + +# X = np.random.normal(size=(2,2)) +# X = MatrixProductState(X, bond_shape=(2,)).decompose(mode='left') + +# y = model.predict(X) +# # model.predict(X) + +# print("Result") +# print(y) +# # print(y.reconstruct()) + +# y_r = np.zeros_like(y.tensor) +# print(y_r) +# print(y.tensor) + +# for k in itertools.product(range(1), range(1)): +# y_r[k] = y.retrieve(k) + +# print(y_r) \ No newline at end of file diff --git a/test/mp.py b/test/mp.py new file mode 100644 index 0000000..3168f47 --- /dev/null +++ b/test/mp.py @@ -0,0 +1,78 @@ +from syngular.tensor.tensor_train import MatrixProductState +from syngular.tensor.tensor_train import MatrixProductOperator + +import numpy as np +import time + +# rnd_mpo = MatrixProductOperator.random((4,4,4), (2,2,2), (3,3,)) +# print(rnd_mpo) +# print(rnd_mpo.to_tensor().shape) + +def test_from_sites(): + mp = MatrixProductOperator.random((4,4,4), (2,2,2), (3,3,)) + sites = mp.sites + + print("Shape", mp.shape) + + print("---------- MatrixProductOperator.from_sites(sites) ----------") + start = time.time() + MatrixProductOperator.from_sites(sites) + end = time.time() + print(f"> Execution time : {end-start:.8f}sec") + print("-------------------------------------------------------") + + + +def test_random(): + + print("---------- MatrixProductOperator.random(input_shape, output_shape, bond_shape) ----------") + start = time.time() + MatrixProductOperator.random((4,4,4), (2,2,2), (3,3,)) + end = time.time() + print(f"> Execution time : {end-start:.8f}sec") + print("-------------------------------------------------------") + + +def test_add(): + x = np.arange(1,17).reshape((2,2,2,2)) + y = np.arange(18,34).reshape((2,2,2,2)) + X = MatrixProductOperator(x, bond_shape=(3,)).decompose() + Y = MatrixProductOperator(y, bond_shape=(3,)).decompose() + + print("---------- MatrixProductOperator.__add__(a, b) ----------") + start = time.time() + Z = X + Y + Y + Y + Y + Y + Y + Y + Y + Y + end = time.time() + print(f"> Execution time : {end-start:.8f}sec") + print("-------------------------------------------------------") + + z = x+y+y+y+y+y+y+y+y+y + + print(x[1,1,1,1]) + print(X[(1,1),(1,1)]) + + print(y[1,1,1,1]) + print(Y[(1,1),(1,1)]) + + print(z[1,1,1,1]) + print(Z[(1,1),(1,1)]) + print(X) + print(Z) + print(Z.to_tensor()) + # print(Z.retrieve((1,1),(1,1))) + +def test_to_tensor(): + tensor = np.arange((2*2*2*2*2*2*2*2)*(2*2*2*2*2*2*2*2)).reshape(2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2) + mpo = MatrixProductOperator(tensor, bond_shape=(2,2,2,2,2,2,2,)).decompose() + print(mpo) + print("---------- MatrixProductOperator.to_tensor() ----------") + start = time.time() + mpo.to_tensor() + end = time.time() + print(f"> Execution time : {end-start:.8f}sec") + print("-------------------------------------------------------") + +# test_from_sites() +# test_random() +# test_to_tensor() +test_add() \ No newline at end of file