From 3d1d739154af95792695ba1b9e9308245a88aade Mon Sep 17 00:00:00 2001 From: tinaoberoi Date: Mon, 8 Jan 2024 12:56:21 -0600 Subject: [PATCH] Add mpo implementation --- scratchpad/qtensor_MPS/mpo.py | 279 +++++++++++++++++++++++++++++++++ scratchpad/qtensor_MPS/mpo2.py | 152 ++++++++++++++++++ scratchpad/qtensor_MPS/mps.py | 133 ++++++++++------ scratchpad/qtensor_MPS/test.py | 88 +++++++++-- 4 files changed, 596 insertions(+), 56 deletions(-) create mode 100644 scratchpad/qtensor_MPS/mpo.py create mode 100644 scratchpad/qtensor_MPS/mpo2.py diff --git a/scratchpad/qtensor_MPS/mpo.py b/scratchpad/qtensor_MPS/mpo.py new file mode 100644 index 00000000..666d7790 --- /dev/null +++ b/scratchpad/qtensor_MPS/mpo.py @@ -0,0 +1,279 @@ +import tensornetwork as tn +import numpy as np +from gates import igate +from typing import List, Union, Text, Optional, Any, Type + + +class MPOLayer: + def __init__(self, tensor_name, N, physical_dim) -> None: + """ + (0) physical dim (0) physical dim (0) physical dim + | | | + @-- (1) bond dim .. (2)bond dim--@-- (1) bond dim ... (1)bond dim--@ + | | | + (2) physical dim (3) physical dim (2) physical dim + """ + self.tensor_name = tensor_name + self.N = N + self.physical_dim = physical_dim + + if N < 2: + raise ValueError("Number of tensors should be >= 2") + # [1.0 0.0 + # 0.0 0.1] + nodes = [ + tn.Node( + np.array([[[1.0, 0.0]], [[0.0, 1.0]]], dtype=np.complex64), + name=tensor_name + str(0), + ) + ] + + for i in range(N - 2): + node = tn.Node( + np.array([[[[1.0, 0.0]]], [[[0.0, 1.0]]]], dtype=np.complex64), + name=tensor_name + str(i + 1), + ) + nodes.append(node) + + nodes.append( + tn.Node( + np.array([[[1.0, 0.0]], [[0.0, 1.0]]], dtype=np.complex64), + name=tensor_name + str(N - 1), + ) + ) + + if N < 3: + tn.connect(nodes[0].get_edge(1), nodes[1].get_edge(1)) + else: + for i in range(1, N - 2): + tn.connect(nodes[i].get_edge(1), nodes[i + 1].get_edge(2)) + tn.connect(nodes[0].get_edge(1), nodes[1].get_edge(2)) + tn.connect(nodes[-1].get_edge(1), nodes[-2].get_edge(1)) + + self._nodes = nodes + + def get_mpo_nodes(self, original) -> list[tn.Node]: + if original: + return self._nodes + + nodes, edges = tn.copy(self._nodes) + return list(nodes.values()) + + def get_mpo_node(self, index, original) -> list[tn.Node]: + return self.get_mpo_nodes(original)[index] + + def construct_mpo(self, gate_function, tensor_name, N, physical_dim=1) -> "MPO": + # IIZ + gate_function = gate_function.reshape([2 * physical_dim] * N) + print(gate_function.shape) + to_split = tn.Node(gate_function, axis_names=[str(i) for i in range(N)]) + + nodes = [] + + for i in range(N - 1): + left_edges = [] + right_edges = [] + + for edge in to_split.get_all_dangling(): + if edge.name == str(i): + temp = tn.split_edge( + edge, (physical_dim, physical_dim), ["a" + str(i), "b" + str(i)] + ) + for e in temp: + left_edges.append(e) + else: + temp = tn.split_edge( + edge, + (physical_dim, physical_dim), + ["a" + str(i + 1), "b" + str(i + 1)], + ) + for e in temp: + right_edges.append(e) + + if nodes: + for edge in nodes[-1].get_all_nondangling(): + if to_split in edge.get_nodes(): + left_edges.append(edge) + + left, right, _ = tn.split_node( + to_split, left_edges, right_edges, left_name=tensor_name + str(i) + ) + nodes.append(left) + to_split = right + to_split.name = tensor_name + str(N - 1) + nodes.append(to_split) + + self._nodes = nodes + + def add_single_qubit_gate(self, gate_matrix, idx): + if idx in range(1, self.N): + gate_matrix = gate_matrix.reshape( + self.physical_dim, 1, 1, self.physical_dim + ) + else: + gate_matrix = gate_matrix.reshape(self.physical_dim, 1, self.physical_dim) + + gate_node = tn.Node( + gate_matrix, + name=self.tensor_name + str(idx + 1), + ) + + self._nodes[idx] = gate_node + + def two_qubit_svd(self, new_node, operating_qubits): + left_connected_edge = None + right_connected_edge = None + + for edge in new_node.get_all_nondangling(): + if self.tensor_name in edge.node1.name: + # Use the "node1" node by default + index = int(edge.node1.name.split(self.tensor_name)[-1]) + else: + # If "node1" is the new_mps_node, use "node2" + index = int(edge.node2.name.split(self.tensor_name)[-1]) + + if index <= operating_qubits[0]: + left_connected_edge = edge + else: + right_connected_edge = edge + + left_edges = [] + right_edges = [] + + # for edge in (left_gate_edge, left_connected_edge): + # if edge != None: + # left_edges.append(edge) + + # for edge in (right_gate_edge, right_connected_edge): + # if edge != None: + # right_edges.append(edge) + + left_edges.append(new_node.get_all_dangling()[0]) + left_edges.append(new_node.get_all_dangling()[2]) + left_edges.append(left_connected_edge) + + right_edges.append(new_node.get_all_dangling()[1]) + right_edges.append(new_node.get_all_dangling()[3]) + right_edges.append(right_connected_edge) + + u, s, vdag, _ = tn.split_node_full_svd( + new_node, left_edges=left_edges, right_edges=right_edges + ) + + new_left = u + new_right = tn.contract_between(s, vdag) + + new_left.name = self._nodes[operating_qubits[0]].name + new_right.name = self._nodes[operating_qubits[1]].name + + self._nodes[operating_qubits[0]] = new_left + self._nodes[operating_qubits[1]] = new_right + + def add_two_qubit_gate(self, gate, operating_qubits): + """ + Method to apply two qubit gates on mpo + + 0 1 + | | + gate + | | + 2 3 + + a b + | | + MPO + | | + c d + + 0 1 + | | + gate' + | | + 2 3 + """ + # transpose + gateT = tn.Node(np.conj(gate.tensor)) + + mpo_indexA = self.get_mpo_node(operating_qubits[0], True).get_all_dangling()[0] + mpo_indexB = self.get_mpo_node(operating_qubits[1], True).get_all_dangling()[0] + mpo_indexC = self.get_mpo_node(operating_qubits[0], True).get_all_dangling()[1] + mpo_indexD = self.get_mpo_node(operating_qubits[1], True).get_all_dangling()[1] + + temp_nodesA = tn.connect(mpo_indexA, gate.get_edge(2)) + temp_nodesB = tn.connect(mpo_indexB, gate.get_edge(3)) + + left_gate_edge = gate.get_edge(0) + right_gate_edge = gate.get_edge(1) + # left_gate_edgeT = gateT.get_edge(2) + # right_gate_edgeT = gateT.get_edge(3) + + new_node = tn.contract_between( + self._nodes[operating_qubits[0]], self._nodes[operating_qubits[1]] + ) + node_gate_edge = tn.flatten_edges_between(new_node, gate) + new_node = tn.contract(node_gate_edge) + + self.two_qubit_svd(new_node, operating_qubits) + + # Try connecting edges of gate 2, 3 + # Check transposition of gate, try multiindexing gate + temp_nodesC = tn.connect(mpo_indexC, gateT.get_edge(0)) + temp_nodesD = tn.connect(mpo_indexD, gateT.get_edge(1)) + + new_node = tn.contract_between( + self._nodes[operating_qubits[0]], self._nodes[operating_qubits[1]] + ) + node_gate_edge = tn.flatten_edges_between(new_node, gateT) + new_node = tn.contract(node_gate_edge) + + self.two_qubit_svd(new_node, operating_qubits) + + def mpo_mps_inner_prod(self, mps): + mpo = self.get_mpo_nodes(False) + + mps_original = mps.__copy__() + mps_original = mps_original.get_mps_nodes(False) + + mps = mps.get_mps_nodes(False) + + for wNode in mps: + wNode.set_tensor(np.conj(wNode.tensor)) + + for i in range(self.N): + tn.connect(mpo[i].get_all_dangling()[0], mps[i].get_all_dangling()[0]) + tn.connect( + mpo[i].get_all_dangling()[0], mps_original[i].get_all_dangling()[0] + ) + + for i in range(self.N - 1): + TW_i = tn.contract_between(mpo[i], mps[i]) + TW_i = tn.contract_between(TW_i, mps_original[i]) + + new_node = tn.contract_between(TW_i, mpo[i + 1]) + mpo[i + 1] = new_node + + end_node = tn.contract_between( + tn.contract_between(mpo[-1], mps[-1]), mps_original[-1] + ) + + inner_prod = np.complex128(end_node.tensor) + return inner_prod + + +class MPO: + def __init__(self, node: tn.Node, idx, physical_dim=2) -> None: + self._node = node + self._indices = idx + self._physical_dim = physical_dim + + def get_node(self, original=True) -> tn.Node: + if original: + return self._node + node_dict, _ = tn.copy([self._node]) + return node_dict[self._node] + + def is_single_qubit_mpo(self) -> True: + return len(self._indices) == 1 + + def is_two_qubit_mpo(self) -> True: + return len(self._indices) == 2 diff --git a/scratchpad/qtensor_MPS/mpo2.py b/scratchpad/qtensor_MPS/mpo2.py new file mode 100644 index 00000000..0f076970 --- /dev/null +++ b/scratchpad/qtensor_MPS/mpo2.py @@ -0,0 +1,152 @@ +import tensornetwork as tn +import numpy as np +from typing import List, Union, Text, Optional, Any, Type +from mps import MPS + + +class MPO2: + def __init__(self, nodes: List, name: Text) -> None: + temp = [] + N = len(nodes) + """ + (2) physical dim (2) physical dim (2) physical dim + | | | + 0 --- @-- (1) bond dim .. (0)bond dim--@-- (1) bond dim ... (0)bond dim--@ --- 1 + | | | + (3) physical dim (3) physical dim (3) physical dim + """ + for node in nodes: + temp.append(tn.Node(node)) + + for i in range(1, N - 2): + tn.connect(temp[i].get_edge(1), temp[i + 1].get_edge(0)) + + if N < 3: + tn.connect(temp[0].get_edge(1), temp[1].get_edge(0)) + else: + tn.connect(temp[0].get_edge(1), temp[1].get_edge(0)) + tn.connect(temp[-2].get_edge(1), temp[-1].get_edge(0)) + + self.tensors = temp + self.name = name + self.N = N + + def get_mpo_nodes(self, original) -> list[tn.Node]: + if original: + return self.tensors + + nodes, edges = tn.copy(self.tensors) + return list(nodes.values()) + + def get_mpo_node(self, index, original) -> list[tn.Node]: + return self.get_mpo_nodes(original)[index] + + def mpo_mps_prod(self, mps: MPS): + mpo = self.get_mpo_nodes(False) + + mps_original = mps.__copy__() + mps_original = mps_original.get_mps_nodes(False) + + mps = mps.get_mps_nodes(False) + + for wNode in mps: + wNode.set_tensor(np.conj(wNode.tensor)) + + # tn.connect(mpo[0].get_all_dangling()[-1], mps[0].get_all_dangling()[0]) + # tn.connect(mpo[0].get_all_dangling()[-1], mps_original[0].get_all_dangling()[0]) + + for i in range(self.N): + tn.connect(mpo[i].get_all_dangling()[-1], mps[i].get_all_dangling()[0]) + tn.connect( + mpo[i].get_all_dangling()[-1], mps_original[i].get_all_dangling()[0] + ) + # tn.connect( + # mpo[i].get_all_dangling()[0], mps_original[i].get_all_dangling()[0] + # ) + + for i in range(self.N - 1): + TW_i = tn.contract_between(mpo[i], mps[i]) + TW_i = tn.contract_between(TW_i, mps_original[i]) + + new_node = tn.contract_between(TW_i, mpo[i + 1]) + mpo[i + 1] = new_node + + end_node = tn.contract_between( + tn.contract_between(mpo[-1], mps[-1]), mps_original[-1] + ) + + inner_prod = np.complex128(end_node.tensor) + return inner_prod + + +class FiniteTFI: + """ + The famous transverse field Ising Hamiltonian. + The ground state energy of the infinite system at criticality is -4/pi. + + Convention: sigma_z=diag([-1,1]) + """ + + def __init__(self, Jx: np.ndarray, Bz: np.ndarray, name: Text = "TFI_MPO") -> None: + """ + Returns the MPO of the finite TFI model. + Args: + Jx: The Sx*Sx coupling strength between nearest neighbor lattice sites. + Bz: Magnetic field on each lattice site. + dtype: The dtype of the MPO. + backend: An optional backend. + name: A name for the MPO. + Returns: + FiniteTFI: The mpo of the infinite TFI model. + """ + dtype = np.complex64 + N = len(Bz) + sigma_x = np.array([[0, 1], [1, 0]]).astype(dtype) + sigma_z = np.diag([-1, 1]).astype(dtype) + mpo = [] + temp = np.zeros(shape=[1, 3, 2, 2], dtype=dtype) + # Bsigma_z + temp[0, 0, :, :] = Bz[0] * sigma_z + # sigma_x + temp[0, 1, :, :] = Jx[0] * sigma_x + # 11 + temp[0, 2, 0, 0] = 1.0 + temp[0, 2, 1, 1] = 1.0 + mpo.append(temp) + for n in range(1, N - 1): + temp = np.zeros(shape=[3, 3, 2, 2], dtype=dtype) + # 11 + temp[0, 0, 0, 0] = 1.0 + temp[0, 0, 1, 1] = 1.0 + # sigma_x + temp[1, 0, :, :] = sigma_x + # Bsigma_z + temp[2, 0, :, :] = Bz[n] * sigma_z + # sigma_x + temp[2, 1, :, :] = Jx[n] * sigma_x + # 11 + temp[2, 2, 0, 0] = 1.0 + temp[2, 2, 1, 1] = 1.0 + mpo.append(temp) + + temp = np.zeros([3, 1, 2, 2], dtype=dtype) + # 11 + temp[0, 0, 0, 0] = 1.0 + temp[0, 0, 1, 1] = 1.0 + # sigma_x + temp[1, 0, :, :] = sigma_x + # Bsigma_z + temp[2, 0, :, :] = Bz[-1] * sigma_z + mpo.append(temp) + + self.tensors = mpo + self.name = name + + +N = 4 +Jz = np.ones(4) +Bz = np.ones(4) +tfi = FiniteTFI(Jz, Bz) +mpo = MPO2(tfi.tensors, tfi.name) +mps1 = MPS("q", N, 2) +print(mpo.mpo_mps_prod(mps1)) diff --git a/scratchpad/qtensor_MPS/mps.py b/scratchpad/qtensor_MPS/mps.py index a8b0e499..ebfb96e0 100644 --- a/scratchpad/qtensor_MPS/mps.py +++ b/scratchpad/qtensor_MPS/mps.py @@ -1,9 +1,13 @@ import tensornetwork as tn +from mpo import MPO, MPOLayer import numpy as np + class MPS: - def __init__(self, tensor_name, N, physical_dim = 1) -> None: + def __init__(self, tensor_name, N, physical_dim=2) -> None: """ + Given: Bond dimension = 1 + (0) physical dim | (1) bond dim --@-- (2) bond dim @@ -16,14 +20,29 @@ def __init__(self, tensor_name, N, physical_dim = 1) -> None: raise ValueError("Number of tensors should be >= 2") # Initialise as |0> = [1.0 0.0 # 0.0 0.0] - nodes = [tn.Node(np.array([[1.0], *[[0.0]]*(physical_dim-1)], dtype=np.complex64), name = tensor_name + str(0))] - for i in range(N-2): - node = tn.Node(np.array([[[1.0]], *[[[0.0]]]*(physical_dim-1)], dtype=np.complex64), name = tensor_name + str(i+1)) + nodes = [ + tn.Node( + np.array([[1.0], *[[0.0]] * (physical_dim - 1)], dtype=np.complex64), + name=tensor_name + str(0), + ) + ] + for i in range(N - 2): + node = tn.Node( + np.array( + [[[1.0]], *[[[0.0]]] * (physical_dim - 1)], dtype=np.complex64 + ), + name=tensor_name + str(i + 1), + ) nodes.append(node) - nodes.append(tn.Node(np.array([[1.0], *[[0.0]]*(physical_dim-1)], dtype=np.complex64), name = tensor_name + str(N-1))) - - for i in range(1, N-2): - tn.connect(nodes[i].get_edge(2), nodes[i+1].get_edge(1)) + nodes.append( + tn.Node( + np.array([[1.0], *[[0.0]] * (physical_dim - 1)], dtype=np.complex64), + name=tensor_name + str(N - 1), + ) + ) + + for i in range(1, N - 2): + tn.connect(nodes[i].get_edge(2), nodes[i + 1].get_edge(1)) if N < 3: tn.connect(nodes[0].get_edge(1), nodes[1].get_edge(1)) @@ -34,18 +53,20 @@ def __init__(self, tensor_name, N, physical_dim = 1) -> None: self._nodes = nodes @staticmethod - def construct_mps_from_wavefunction(wavefunction, tensor_name, N, physical_dim = 1) -> 'MPS': + def construct_mps_from_wavefunction( + wavefunction, tensor_name, N, physical_dim=1 + ) -> "MPS": """ Method to create wavefunction from mps """ if wavefunction.size != physical_dim**N: raise ValueError() - - wavefunction = np.reshape(wavefunction, [physical_dim]*N) + + wavefunction = np.reshape(wavefunction, [physical_dim] * N) to_split = tn.Node(wavefunction, axis_names=[str(i) for i in range(N)]) nodes = [] - for i in range(N-1): + for i in range(N - 1): left_edges = [] right_edges = [] @@ -54,27 +75,29 @@ def construct_mps_from_wavefunction(wavefunction, tensor_name, N, physical_dim = left_edges.append(edge) else: right_edges.append(edge) - + if nodes: for edge in nodes[-1].get_all_nondangling(): if to_split in edge.get_nodes(): left_edges.append(edge) - left, right, _ = tn.split_node(to_split, left_edges, right_edges, left_name=tensor_name+str(i)) + left, right, _ = tn.split_node( + to_split, left_edges, right_edges, left_name=tensor_name + str(i) + ) nodes.append(left) to_split = right - to_split.name = tensor_name + str(N-1) + to_split.name = tensor_name + str(N - 1) nodes.append(to_split) mps = MPS(tensor_name, N, physical_dim) mps._nodes = nodes return mps - + @property def N(self): return self._N - + @property def physical_dim(self): return self._physical_dim @@ -82,13 +105,13 @@ def physical_dim(self): def get_mps_nodes(self, original) -> list[tn.Node]: if original: return self._nodes - + nodes, edges = tn.copy(self._nodes) return list(nodes.values()) - + def get_mps_node(self, index, original) -> tn.Node: return self.get_mps_nodes(original)[index] - + def get_tensors(self, original) -> list[tn.Node]: nodes = [] @@ -98,10 +121,10 @@ def get_tensors(self, original) -> list[tn.Node]: nodes, edges = tn.copy(self._nodes) return list([node.tensor for node in nodes]) - + def get_tensor(self, index, original) -> tn.Node: return self.get_tensors(original)[index] - + def is_unitary(): pass @@ -112,9 +135,9 @@ def get_wavefunction(self) -> np.array: for node in nodes: curr = tn.contract_between(curr, node) - wavefunction = np.reshape(curr.tensor, newshape=(self._physical_dim ** self._N)) + wavefunction = np.reshape(curr.tensor, newshape=(self._physical_dim**self._N)) return wavefunction - + def apply_single_qubit_gate(self, gate, index) -> None: """ Method to apply single qubit gate on mps @@ -136,7 +159,7 @@ def apply_single_qubit_gate(self, gate, index) -> None: new_node = tn.contract(temp_node, name=self._nodes[index].name) self._nodes[index] = new_node - + def apply_two_qubit_gate(self, gate, operating_qubits): """ Method to apply two qubit gates on mps @@ -153,19 +176,24 @@ def apply_two_qubit_gate(self, gate, operating_qubits): """ # reshape the gate to order-4 tensor # Assumimg operating_qubits are sorted - mps_indexA = self.get_mps_node(operating_qubits[0], True).get_all_dangling().pop() - mps_indexB = self.get_mps_node(operating_qubits[1], True).get_all_dangling().pop() + mps_indexA = ( + self.get_mps_node(operating_qubits[0], True).get_all_dangling().pop() + ) + mps_indexB = ( + self.get_mps_node(operating_qubits[1], True).get_all_dangling().pop() + ) temp_nodesA = tn.connect(mps_indexA, gate.get_edge(2)) temp_nodesB = tn.connect(mps_indexB, gate.get_edge(3)) left_gate_edge = gate.get_edge(0) right_gate_edge = gate.get_edge(1) - new_node = tn.contract_between(self._nodes[operating_qubits[0]], self._nodes[operating_qubits[1]]) + new_node = tn.contract_between( + self._nodes[operating_qubits[0]], self._nodes[operating_qubits[1]] + ) node_gate_edge = tn.flatten_edges_between(new_node, gate) new_node = tn.contract(node_gate_edge) - left_connected_edge = None right_connected_edge = None @@ -194,9 +222,7 @@ def apply_two_qubit_gate(self, gate, operating_qubits): right_edges.append(edge) u, s, vdag, _ = tn.split_node_full_svd( - new_node, - left_edges=left_edges, - right_edges=right_edges + new_node, left_edges=left_edges, right_edges=right_edges ) new_left = u @@ -208,7 +234,6 @@ def apply_two_qubit_gate(self, gate, operating_qubits): self._nodes[operating_qubits[0]] = new_left self._nodes[operating_qubits[1]] = new_right - def inner_product(self, wMPS): """ Method to calculate inner product of mps @@ -222,14 +247,14 @@ def inner_product(self, wMPS): for i in range(self.N): tn.connect(T[i].get_all_dangling().pop(), W[i].get_all_dangling().pop()) - - for i in range(self.N-1): + + for i in range(self.N - 1): TW_i = tn.contract_between(T[i], W[i]) - new_node = tn.contract_between(TW_i, W[i+1]) - W[i+1] = new_node + new_node = tn.contract_between(TW_i, W[i + 1]) + W[i + 1] = new_node return np.complex128((tn.contract_between(T[-1], W[-1])).tensor) - + def get_norm(self): """ Method to calculate norm of mps @@ -237,10 +262,9 @@ def get_norm(self): val = np.sqrt(self.inner_product(self).real) return val - def left_cannoise(self, i): nodes = [] - for i in range(i): + for i in range(i): left_edges = [] right_edges = [] @@ -249,13 +273,15 @@ def left_cannoise(self, i): left_edges.append(edge) else: right_edges.append(edge) - + if nodes: for edge in nodes[-1].get_all_nondangling(): if to_split in edge.get_nodes(): left_edges.append(edge) - left, right, _ = tn.split_node(to_split, left_edges, right_edges, left_name="q"+str(i)) + left, right, _ = tn.split_node( + to_split, left_edges, right_edges, left_name="q" + str(i) + ) nodes.append(left) to_split = right @@ -266,17 +292,32 @@ def __copy__(self): copy_mps = MPS(self.name, self._N, self._physical_dim) copy_mps._nodes = self.get_mps_nodes(original=False) return copy_mps - + def get_expectation(self, observable, idx): mps_copy = self.__copy__() mps_copy.apply_single_qubit_gate(observable, idx) return self.inner_product(mps_copy) + def apply_mpo(self, operation: MPO) -> None: + if operation.is_single_qubit_mpo(): + self.apply_single_qubit_gate(operation._node, operation._indices[0]) + else: + self.apply_two_qubit_gate(operation._node, *operation._indices) + def apply_mpo_layer(self, mpo: MPOLayer) -> None: + # Assume self.N == operation.N + mpo = mpo.get_mpo_nodes(False) + mps = self.get_mps_nodes(False) + for i in range(self.N): + tn.connect(mpo[i].get_all_dangling()[1], mps[i].get_all_dangling()[0]) + nodes = [] + for i in range(self.N): + new_node = tn.contract_between(mpo[i], mps[i]) + nodes.append(new_node) + for i in range(self.N - 1): + _ = tn.flatten_edges_between(nodes[i], nodes[i + 1]) - - - + self._nodes = nodes diff --git a/scratchpad/qtensor_MPS/test.py b/scratchpad/qtensor_MPS/test.py index 4953a95f..c88cf459 100644 --- a/scratchpad/qtensor_MPS/test.py +++ b/scratchpad/qtensor_MPS/test.py @@ -1,8 +1,9 @@ import numpy as np -import tensornetwork as tn -import xyzpy as xyz -from gates import xgate, cnot, hgate, zgate +from gates import xgate, cnot, hgate, zgate, igate from mps import MPS +from mpo import MPO, MPOLayer +import tensornetwork as tn +from constants import xmatrix, cnot_matrix def test_from_wavefunction_all_zero_state(): @@ -43,15 +44,15 @@ def test_apply_one_qubit_mps_operation_xgate(): def test_apply_twoq_cnot_two_qubits(): # In the following tests, the first qubit is always the control qubit. # Check that CNOT|10> = |11> - mps = MPS("q", 2, 2) + mps = MPS("q", 4, 2) assert np.isclose(mps.get_norm(), 1.0) - mps.apply_single_qubit_gate(xgate(), 0) - mps.apply_two_qubit_gate(cnot(), [0, 1]) - assert np.isclose(mps.get_norm(), 1.0) - assert np.allclose( - mps.get_wavefunction(), np.array([0.0, 0.0, 0.0, 1.0], dtype=np.complex64) - ) + # mps.apply_single_qubit_gate(xgate(), 0) + mps.apply_two_qubit_gate(cnot(), [1, 2]) + # assert np.isclose(mps.get_norm(), 1.0) + # assert np.allclose( + # mps.get_wavefunction(), np.array([0.0, 0.0, 0.0, 1.0], dtype=np.complex64) + # ) def test_apply_two_twoq_cnot_two_qubits(): @@ -146,3 +147,70 @@ def test_expectation_value_xgate_at_k(): expectation_array.append(mps.get_expectation(zgate(), i)) np.allclose(expectation_array, [1.0, 1.0, 1.0, -1.0, 1.0]) + + +def test_mps_mpo(): + mps00 = MPS("q", 2, 2) + mps01 = MPS("q", 2, 2) + mps10 = MPS("q", 2, 2) + mps11 = MPS("q", 2, 2) + + mps01.apply_mpo(MPO(xgate(), [1])) + mps10.apply_mpo(MPO(xgate(), [0])) + mps11.apply_mpo(MPO(xgate(), [0])) + mps11.apply_mpo(MPO(xgate(), [1])) + + allmps = (mps00, mps01, mps10, mps11) + + # Test inner products + for i in range(4): + for j in range(4): + assert np.isclose(allmps[i].inner_product(allmps[j]), i == j) + + +def test_mps_single_qubit_mpo_layer(): + mps1 = MPS("q", 3, 2) + mps1.apply_single_qubit_gate(xgate(), 1) + assert mps1.get_norm() == 1 + # expectation = mps1(original) inner prod mps1(modified) + mps2 = MPS("q", 3, 2) + mpo_layer = MPOLayer("q", 3, 2) + mpo_layer.add_single_qubit_gate(xmatrix, 1) + mps2.apply_mpo_layer(mpo_layer) + assert mps2.get_norm() == 1 + + assert (mps1.get_wavefunction() == mps2.get_wavefunction()).all() + + +def test_mpo_two_qubit(): + mpo_layer = MPOLayer("q", 4, 2) + # print("BEFORE") + # print(mpo_layer._nodes) + mpo_layer.add_two_qubit_gate(cnot(), [1, 2]) + print("AFTER") + print(mpo_layer._nodes) + # assert (mps1.get_wavefunction() == mps2.get_wavefunction()).all() + + +def test_mpo_construction_from_gate_function(): + # IZI + I = np.eye(2) + Z = np.array([[0.0, 1.0], [1.0, 0.0]], dtype=np.complex64) + H = 1 / np.sqrt(2) * np.array([[1.0, 1.0], [1.0, -1.0]], dtype=np.complex64) + gate_func = np.kron(I, I) + print(gate_func.shape) + mpo = MPOLayer("q", 2, 2) + mpo.construct_mpo(gate_func, "q", 2, 2) + mps = MPS("q", 2, 2) + print(mps.get_wavefunction()) + mps.apply_mpo_layer(mpo) + print(mps.get_wavefunction()) + + +# TODO: +# Verify for two qubit gates. - (SVD and then apply it) (apply and then SVD) +# Evolve half od circuit in mps and half in mpo. +# Evaluate Bell state, epectation value of X at a particular index +# Evaluate X, Y, Z observable +# Check for one/two/three qubits (cover edge cases) for both mps/mpo +test_mpo_construction_from_gate_function()