diff --git a/quantum/gate.py b/quantum/gate.py index 22e2b6f..85fd0b2 100644 --- a/quantum/gate.py +++ b/quantum/gate.py @@ -24,4 +24,11 @@ H = 1/np.sqrt(2) * np.array([ [1., 1.], [1., -1.] -]) \ No newline at end of file +]) + +CX = np.array([ + [1., 0., 0., 0.], + [0., 1., 0., 0.], + [0., 0., 0., 1.], + [0., 0., 1., 0.], +]).reshape((2,2,2,2)) \ No newline at end of file diff --git a/quantum/qbit.py b/quantum/qbit.py index b818fbe..c1e9f0d 100644 --- a/quantum/qbit.py +++ b/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) \ No newline at end of file + 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 \ No newline at end of file diff --git a/tensor/matrix_product_operator.py b/tensor/matrix_product_operator.py index ae0d130..af75912 100644 --- a/tensor/matrix_product_operator.py +++ b/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": diff --git a/tensor/matrix_product_state.py b/tensor/matrix_product_state.py index 263b1be..86be0ea 100644 --- a/tensor/matrix_product_state.py +++ b/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): diff --git a/test/test_mps.py b/test/test_mps.py index 7471654..d37095c 100644 --- a/test/test_mps.py +++ b/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() \ No newline at end of file +# test_augment() +test_apply() \ No newline at end of file diff --git a/test/test_quantum.py b/test/test_quantum.py index 2807d80..adf2679 100644 --- a/test/test_quantum.py +++ b/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 diff --git a/times.txt b/times.txt new file mode 100644 index 0000000..03633e0 Binary files /dev/null and b/times.txt differ