Fix two qubit gate for mpo

scratchpad/qtensor_MPS/constants.py

@@ -1,8 +1,7 @@
 import numpy as np
 
-hmatrix = (
-    1 / np.sqrt(2) * np.array([[1.0, 1.0], [1.0, -1.0]], dtype=np.complex64)
-)
+pauli_matrices = ["X", "Z", "I", "Y"]
+hmatrix = 1 / np.sqrt(2) * np.array([[1.0, 1.0], [1.0, -1.0]], dtype=np.complex64)
 imatrix = np.array([[1.0, 0.0], [0.0, 1.0]], dtype=np.complex64)
 xmatrix = np.array([[0.0, 1.0], [1.0, 0.0]], dtype=np.complex64)
 ymatrix = np.array([[0.0, -1j], [1j, 0.0]], dtype=np.complex64)
@@ -28,10 +27,10 @@
 )
 swap_matrix = np.reshape(_swap_matrix, newshape=(2, 2, 2, 2))
 
-_hbar = 1.0545718*10e-34
-_sigma_z =  _hbar * 0.5 * zmatrix
-_sigma_x_pos =  _hbar * 0.5 * (xmatrix + 1j*ymatrix)
-_sigma_x_neg =  _hbar * 0.5 * (xmatrix - 1j*ymatrix)
+_hbar = 1.0545718 * 10e-34
+_sigma_z = _hbar * 0.5 * zmatrix
+_sigma_x_pos = _hbar * 0.5 * (xmatrix + 1j * ymatrix)
+_sigma_x_neg = _hbar * 0.5 * (xmatrix - 1j * ymatrix)
 
 _sigma_z_sigma_z_gate_matrix = np.kron(_sigma_z, _sigma_z)
 sigma_z_sigma_z_gate_matrix = _sigma_z_sigma_z_gate_matrix
@@ -40,10 +39,14 @@
 # TODO: What does it look if hamiltonian is e^(I) / e^(0)
 
 _sigma_x_pos_sigma_x_neg_gate_matrix = np.kron(_sigma_x_pos, _sigma_x_neg)
-sigma_x_pos_sigma_x_neg_gate_matrix = np.reshape(_sigma_x_pos_sigma_x_neg_gate_matrix, newshape=(2,2,2,2))
+sigma_x_pos_sigma_x_neg_gate_matrix = np.reshape(
+    _sigma_x_pos_sigma_x_neg_gate_matrix, newshape=(2, 2, 2, 2)
+)
 
 _sigma_x_neg_sigma_x_pos_gate_matrix = np.kron(_sigma_x_neg, _sigma_x_pos)
-sigma_x_neg_sigma_x_pos_gate_matrix = np.reshape(_sigma_x_neg_sigma_x_pos_gate_matrix, newshape=(2,2,2,2))
+sigma_x_neg_sigma_x_pos_gate_matrix = np.reshape(
+    _sigma_x_neg_sigma_x_pos_gate_matrix, newshape=(2, 2, 2, 2)
+)
 # _ising_hamiltonian_matrix = sigma_z_sigma_z_gate_matrix + 0.5 * sigma_x_pos_sigma_x_neg_gate_matrix + 0.5 * sigma_x_neg_sigma_x_pos_gate_matrix
 _ising_hamiltonian_matrix = sigma_z_sigma_z_gate_matrix
-ising_hamiltonian_matrix = np.reshape(_ising_hamiltonian_matrix, newshape=(2,2,2,2))
+ising_hamiltonian_matrix = np.reshape(_ising_hamiltonian_matrix, newshape=(2, 2, 2, 2))

scratchpad/qtensor_MPS/gates.py

@@ -6,6 +6,11 @@
 from scipy.linalg import qr
 
 
+def get_gate(gc):
+    gates = {"X": xgate, "I": igate, "Y": ygate, "Z": zgate}
+    return gates.get(gc, lambda: None)().tensor
+
+
 def xgate() -> tn.Node:
     return tn.Node(deepcopy(xmatrix), name="xgate")
 
@@ -15,7 +20,11 @@ def igate() -> tn.Node:
 
 
 def zgate() -> tn.Node:
-    return tn.Node(deepcopy(zmatrix), name="xgate")
+    return tn.Node(deepcopy(zmatrix), name="zgate")
+
+
+def ygate() -> tn.Node:
+    return tn.Node(deepcopy(ymatrix), name="ygate")
 
 
 def cnot() -> tn.Node:
@@ -29,7 +38,7 @@ def hgate() -> tn.Node:
 def sigmaRx(t) -> tn.Node:
     _hbar = 1.0545718 * 10e-34
     # theta =  -1j * t * 0.5 * _hbar * 0.5
-    theta = 1j * 0.1
+    theta = 0.1 * t
 
     gate_matrix = np.array(
         [
@@ -45,13 +54,13 @@ def sigmaRzz(t) -> tn.Node:
     _hbar = 1.0545718 * 10e-34
     # theta =  -1j * t * 0.5 * _hbar * 0.5 * _hbar * 0.5
     # J = 0.1Hz
-    theta = -1j * 0.7
+    theta = -1j * 0.7 * t
 
     gate_matrix = np.array(
         [
             [np.exp(theta / 2), 0.0, 0.0, 0.0],
-            [0.0, np.exp(theta / 2), 0.0, 0.0],
-            [0.0, 0.0, np.exp(theta / 2), 0.0],
+            [0.0, np.exp(-1 * theta / 2), 0.0, 0.0],
+            [0.0, 0.0, np.exp(-1 * theta / 2), 0.0],
             [0.0, 0.0, 0.0, np.exp(theta / 2)],
         ]
     )

scratchpad/qtensor_MPS/mpo.py

@@ -1,7 +1,9 @@
 import tensornetwork as tn
 import numpy as np
-from gates import igate
+from gates import igate, get_gate, xgate
+from constants import pauli_matrices, xmatrix
 from typing import List, Union, Text, Optional, Any, Type
+from copy import deepcopy
 
 
 class MPOLayer:
@@ -62,62 +64,87 @@ def get_mpo_nodes(self, original) -> list[tn.Node]:
     def get_mpo_node(self, index, original) -> list[tn.Node]:
         return self.get_mpo_nodes(original)[index]
 
-    def construct_mpo(self, gate_function) -> "MPO":
-        # IIZ
-        gate_function = gate_function.reshape([self.physical_dim] * 2 * self.N)
-        to_split = tn.Node(
-            gate_function,
-            axis_names=["u" + str(i) for i in range(self.N)]
-            + ["d" + str(i) for i in range(self.N)],
-        )
+    def construct_mpo(self, pauli_string) -> "MPO":
+        # check if all elements are pauli matrices
+        paulis = list(pauli_string)
+
+        if any(pstring not in pauli_matrices for pstring in paulis):
+            print("Error: Not all are pauli string")
+            # Comments: change gatefunction - paulistring
+            pauli_string = pauli_string.reshape([self.physical_dim] * 2 * self.N)
+            to_split = tn.Node(
+                pauli_string,
+                axis_names=["u" + str(i) for i in range(self.N)]
+                + ["d" + str(i) for i in range(self.N)],
+            )
 
-        nodes = []
+            nodes = []
 
-        for i in range(self.N - 1):
-            left_edges = []
-            right_edges = []
+            for i in range(self.N - 1):
+                left_edges = []
+                right_edges = []
 
-            for edge in to_split.get_all_dangling():
-                if edge.name == "u" + str(i) or edge.name == "d" + str(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():
+                for edge in to_split.get_all_dangling():
+                    if edge.name == "u" + str(i) or edge.name == "d" + str(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=self.tensor_name + str(i),
+                    # max_singular_values=1,
+                )
+                nodes.append(left)
+                to_split = right
+            to_split.name = self.tensor_name + str(self.N - 1)
+            nodes.append(to_split)
+
+            self._nodes = nodes
 
-            left, right, _ = tn.split_node(
-                to_split,
-                left_edges,
-                right_edges,
-                left_name=self.tensor_name + str(i),
-                max_singular_values=1,
-            )
-            nodes.append(left)
-            to_split = right
-        to_split.name = self.tensor_name + str(self.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)
+            if len(pauli_string) != self.N:
+                # pad pauli string
+                pauli_string += "I" * (self.N - len(pauli_string))
 
-        gate_node = tn.Node(
-            gate_matrix,
-            name=self.tensor_name + str(idx + 1),
-        )
+            nodes = self.get_mpo_nodes(False)
+            paulis = list(pauli_string)
 
-        self._nodes[idx] = gate_node
+            for i, (node, pauli) in enumerate(zip(nodes, paulis)):
+                pauli_gate = get_gate(pauli)
 
-    def two_qubit_svd(self, new_node, operating_qubits):
+                if i == 0 or i == self.N - 1:
+                    pauli_gate = pauli_gate.reshape(
+                        self.physical_dim, 1, self.physical_dim
+                    )
+                else:
+                    pauli_gate = pauli_gate.reshape(
+                        self.physical_dim, 1, 1, self.physical_dim
+                    )
+
+                node.set_tensor(pauli_gate)
+
+            self._nodes = nodes
+
+    def add_single_qubit_gate(self, gate, idx):
+        node = self._nodes[idx]
+        lst = list(node.get_all_dangling())
+        mps_index_edge = lst[0]
+        gate_edge = gate[1]
+        temp_node = tn.connect(mps_index_edge, gate_edge)
+        new_node = tn.contract(temp_node, name=self._nodes[idx].name)
+        self._nodes[idx] = new_node
+
+    def two_qubit_svd(
+        self, new_node, operating_qubits, left_gate_edge, right_gate_edge
+    ):
         left_connected_edge = None
         right_connected_edge = None
 
@@ -137,21 +164,13 @@ def two_qubit_svd(self, new_node, operating_qubits):
         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)
+        for edge in (*left_gate_edge, left_connected_edge):
+            if edge != None:
+                left_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)
+        for edge in (*right_gate_edge, right_connected_edge):
+            if edge != None:
+                right_edges.append(edge)
 
         u, s, vdag, _ = tn.split_node_full_svd(
             new_node, left_edges=left_edges, right_edges=right_edges
@@ -182,11 +201,7 @@ def add_two_qubit_gate(self, gate, operating_qubits):
             |   |
             c   d
 
-            0  1
-            |  |
-            gate'
-            |  |
-            2  3
+
         """
         # transpose
         gateT = tn.Node(np.conj(gate.tensor))
@@ -199,31 +214,38 @@ def add_two_qubit_gate(self, gate, operating_qubits):
         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_edge = [gate.get_edge(0)]
+        right_gate_edge = [gate.get_edge(1)]
+
+        left_gate_edge.append(mpo_indexC)
+        right_gate_edge.append(mpo_indexD)
         # 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]]
         )
+
+        x = self.get_mpo_node(operating_qubits[0], True).get_all_dangling()[1]
+        y = self.get_mpo_node(operating_qubits[1], True).get_all_dangling()[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)
+        self.two_qubit_svd(new_node, operating_qubits, left_gate_edge, right_gate_edge)
 
         # 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))
+        # 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)
+        # 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)
+        # self.two_qubit_svd(new_node, operating_qubits)
 
     def mpo_mps_inner_prod(self, mps):
         mpo = self.get_mpo_nodes(False)