From d943536652fca758a69c8a46255c7772221f6d68 Mon Sep 17 00:00:00 2001 From: Antoine Debouchage Date: Mon, 27 Dec 2021 21:59:22 +0100 Subject: [PATCH] Little change and bug in model prediction after update of weights to be fixed --- core/model.py | 46 ++++++++----- quantum/gate.py | 5 ++ quantum/qbit.py | 44 +++++++++++-- tensor/matrix_product_operator.py | 45 +++++++++---- tensor/matrix_product_state.py | 14 +++- test/model.py | 36 ----------- test/test_model.py | 82 +++++++++++++++++++++++ test/test_quantum.py | 104 +++++++++++++++++++++++++++--- 8 files changed, 293 insertions(+), 83 deletions(-) delete mode 100644 test/model.py create mode 100644 test/test_model.py diff --git a/core/model.py b/core/model.py index d003b17..7d678b9 100644 --- a/core/model.py +++ b/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): diff --git a/quantum/gate.py b/quantum/gate.py index 3aacb4c..0436232 100644 --- a/quantum/gate.py +++ b/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.], diff --git a/quantum/qbit.py b/quantum/qbit.py index 25b8721..60f9937 100644 --- a/quantum/qbit.py +++ b/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 \ No newline at end of file diff --git a/tensor/matrix_product_operator.py b/tensor/matrix_product_operator.py index 7923711..18cd920 100644 --- a/tensor/matrix_product_operator.py +++ b/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) diff --git a/tensor/matrix_product_state.py b/tensor/matrix_product_state.py index 6c4d56a..87ee2d2 100644 --- a/tensor/matrix_product_state.py +++ b/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:] \ No newline at end of file diff --git a/test/model.py b/test/model.py deleted file mode 100644 index ac7cc6b..0000000 --- a/test/model.py +++ /dev/null @@ -1,36 +0,0 @@ -import itertools -from syngular.core import Model, Output -from syngular.core import Linear -from syngular.tensor import MatrixProductState, MatrixProductOperator - -import numpy as np - -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/test_model.py b/test/test_model.py new file mode 100644 index 0000000..117df99 --- /dev/null +++ b/test/test_model.py @@ -0,0 +1,82 @@ +import itertools +from syngular.core import Model, Output +from syngular.core import Linear +from syngular.tensor import MatrixProductState, MatrixProductOperator + +import numpy as np + +w = np.arange(2**4).reshape(2,2,2,2) +W = MatrixProductOperator(w, bond_shape=(2,)) +W.decompose() + +model = Model([ + Linear(4,4,core=2, bond=4, weights_initializer=W), + Linear(4,4,core=2, bond=4, weights_initializer=W), + # Linear(4,4,core=2, bond=4, weights_initializer=W), + # Linear(4,4,core=2, bond=4, weights_initializer=W), + # Linear(4,4,core=2, bond=4, weights_initializer=W), + # Linear(4,4,core=2, bond=4, weights_initializer=W), + # Linear(4,4,core=2, bond=4, weights_initializer=W), + # Linear(4,4,core=2, bond=4, weights_initializer=W), + # Linear(4,4,core=2, bond=4, weights_initializer=W), + # Linear(4,4,core=2, bond=4, weights_initializer=W), + # Linear(4,4,core=2, bond=4, weights_initializer=W), + # Linear(4,4,core=2, bond=4, weights_initializer=W), + # Linear(4,4,core=2, bond=4, weights_initializer=W), + # Linear(4,4,core=2, bond=4, weights_initializer=W), + # Linear(4,4,core=2, bond=4, weights_initializer=W), + # Linear(4,4,core=2, bond=4, weights_initializer=W), + # Linear(4,4,core=2, bond=4, weights_initializer=W), + # Linear(4,4,core=2, bond=4, weights_initializer=W), + # Linear(4,4,core=2, bond=4, weights_initializer=W), + # Linear(4,4,core=2, bond=4, weights_initializer=W), + # Linear(4,4,core=2, bond=4, weights_initializer=W), + # Linear(4,4,core=2, bond=4, weights_initializer=W), + # Linear(4,4,core=2, bond=4, weights_initializer=W), + # Linear(4,4,core=2, bond=4, weights_initializer=W), + # Linear(4,4,core=2, bond=4, weights_initializer=W), + # Linear(4,4,core=2, bond=4, weights_initializer=W), + # Linear(4,4,core=2, bond=4, weights_initializer=W), + # Linear(4,4,core=2, bond=4, weights_initializer=W), + # Linear(4,4,core=2, bond=4, weights_initializer=W), + # Linear(4,8,core=2,bond=2), + # Output((2,)) +]) +model.build() + +print("---") +print(model.draw()) +print("---") + +# X = np.random.normal(size=(2,2)) +# X = MatrixProductState(X, bond_shape=(2,)).decompose(mode='left') + +x = np.arange(2**2).reshape(2,2) +X = MatrixProductState(x, bond_shape=(2,)).decompose() +Y = MatrixProductState(np.arange(4).reshape(2,2), bond_shape=(2,)).decompose() +# X = MatrixProductState.random((2,2), (2,)) +y = model.predict(X).to_tensor() +print(y) +print(Y.to_tensor()) + +model.train(X, Y, epochs=2) +print("X", X.to_tensor()) +y = (model.predict(X)) +# print(y) + +# print(model.layers[0].trainable_tensor_weights[0]["weight"].sites) + +# print(x.reshape((4,)) @ w.reshape(4,4)) + +# 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/test_quantum.py b/test/test_quantum.py index 7b884ac..a563ae7 100644 --- a/test/test_quantum.py +++ b/test/test_quantum.py @@ -75,34 +75,118 @@ def test_swap(): # circ = Circuit(2) # circ.add((gate.X, 0)) # circ.add((gate.SWAP, 0)) -test_swap() +# test_swap() -def test_bernstein_vazirani(size): - token = "11".zfill(size)[::-1] - print(token) +def test_bernstein_vazirani(token): + size = len(token) + # token = token.zfill(size)[::-1] + # print(token) - circ = Circuit(size=size) + circ = Circuit(size=size+1) + circ.add((gate.H, size)) + circ.add((gate.Z, size)) for i in range(size): circ.add((gate.H, i)) - + # circ.add((gate.X, 0)) + # circ.add((gate.CX, size, 0)) # circ.add(BlackBox()) for i in range(size): if token[i] == '1': print("ok", i) - circ.add((gate.CX, i)) + circ.add((gate.CX, size, i)) for i in range(size): circ.add((gate.H, i)) - circ.run() - - print(circ.current_state.state | Qbit(size).state) + # circ.run() + for stp in range(len(circ.structure)): + print('------------ step -------------') + print(circ.current_state.to_tensor()) + circ.step() print(circ.current_state.to_tensor()) + # print("ok") + + import numpy as np + print(np.argmax(circ.current_state.to_tensor())) + + + # qbit_token = Qbit.from_binary(token) + # print(circ.current_state.state | Qbit(size).state) + # print(circ.current_state.state | qbit_token.state) + + +# test_bernstein_vazirani("11") + +def test_hadamard(token): + + token = token[::-1] + size = len(token) + + qbit = Qbit(size+1) + + qbit @= (gate.H, size) + qbit @= (gate.Z, size) + + for i in range(size): qbit @= (gate.H, i) + + for i in range(len(token)): + if token[i] == "1": qbit @= (gate.CX, 2, i) + else: qbit @= (gate.I, i) + + for i in range(size): qbit @= (gate.H, i) + + print(qbit.to_tensor()) + + import numpy as np + print(np.argmax(qbit.to_tensor())) + print("Probability : ") + print( + (qbit.state | Qbit.from_binary(token[::-1]+"1").state)**2 + + (qbit.state | Qbit.from_binary(token[::-1]+"0").state)**2 + ) + + +# test_hadamard("111") + +def test_bell_state(): + start = time.time() + qbit = Qbit(2) + qbit @= (gate.H, 0) + qbit @= (gate.CX, 0, 1) ## Control 1-th qbit with 0-th qbit + end = time.time() + print(f"> Execution time : {end-start:.8f}sec") + + print(qbit.to_tensor()) + +# test_bell_state() + +def test_hzh_x(): + qbit = Qbit(1) + qbit @= gate.X + print(qbit.to_tensor()) + + qbit = Qbit(1) + qbit @= gate.H @ gate.Z @ gate.H + print(qbit.to_tensor()) + +def test_ss_z(): + qbit = Qbit(1) + qbit @= gate.X + qbit @= gate.Z + print(qbit.to_tensor()) + + qbit = Qbit(1) + qbit @= gate.X + qbit @= gate.S @ gate.S + print(qbit.to_tensor()) + +test_hzh_x() +test_ss_z() # for i in range(1): # test_bernstein_vazirani(4) # size = 1