Merge pull request #120 from WanJiawang/wjw

the ansatz evolution according to hamiltonian

quafu/synthesis/evolution.py

@@ -19,29 +19,194 @@
 import quafu.elements.element_gates as qeg
 
 
+def single_qubit_evol(pauli: str, time: float):
+    """
+    Args:
+        pauli: Pauli string (little endian convention)
+        time: Evolution time
+    """
+    reversed_pauli = pauli[::-1]
+    gates = []
+
+    for i, pauli_i in enumerate(reversed_pauli):
+        if pauli_i == "I":
+            continue
+        elif pauli_i == "X":
+            gates.append(qeg.RXGate(i, 2 * time))
+            return gates
+        elif pauli_i == "Y":
+            gates.append(qeg.RYGate(i, 2 * time))
+            return gates
+        elif pauli_i == "Z":
+            gates.append(qeg.RZGate(i, 2 * time))
+            return gates
+        else:
+            raise NotImplementedError("Pauli string not yet supported")
+
+
 def two_qubit_evol(pauli: str, time: float, cx_structure: str = "chain"):
     """
     Args:
         pauli: Pauli string (little endian convention)
         time: Evolution time
-        cx_structure: TODO
+        cx_structure: Determine the structure of CX gates, can be either "chain" for
+            next-neighbor connections or "fountain" to connect directly to the top qubit.
     """
     reversed_pauli = pauli[::-1]
     qubits = [i for i in range(len(reversed_pauli)) if reversed_pauli[i] != "I"]
     labels = np.array([reversed_pauli[i] for i in qubits])
     gates = []
 
     if all(labels == "X"):
-        pass
+        gates.append(qeg.RXXGate(qubits[0], qubits[1], 2 * time))
     elif all(labels == "Y"):
-        pass
+        gates.append(qeg.RYYGate(qubits[0], qubits[1], 2 * time))
     elif all(labels == "Z"):
         gates.append(qeg.RZZGate(qubits[0], qubits[1], 2 * time))
     else:
+        return multi_qubit_evol(pauli, time, cx_structure)
+    return gates
+
+
+def multi_qubit_evol(pauli: str, time: float, cx_structure: str = "chain"):
+    # determine whether the Pauli string consists of Pauli operators
+    if not all(pauli_char in "XYZI" for pauli_char in pauli):
         raise NotImplementedError("Pauli string not yet supported")
+    gates = []
+    # get diagonalizing clifford gate list
+    cliff, cliff_inverse = diagonalizing_clifford(pauli)
+
+    # get CX chain to reduce the evolution to the top qubit
+    if cx_structure == "chain":
+        chain, chain_inverse = cnot_chain(pauli)
+    else:
+        chain, chain_inverse = cnot_fountain(pauli)
+
+    # determine qubit to do the rotation on
+    target = None
+    # Note that all phases are removed from the pauli label and are only in the coefficients.
+    # That's because the operators we evolved have all been translated to a SparsePauliOp.
+    for i, pauli_i in enumerate(pauli[::-1]):
+        if pauli_i != "I":
+            target = i
+            break
+
+    # build the evolution as: diagonalization, reduction, 1q evolution, followed by inverses
+    gates.extend(cliff)
+    gates.extend(chain)
+    gates.append(qeg.RZGate(target, 2 * time))
+    gates.extend(chain_inverse)
+    gates.extend(cliff_inverse)
+
     return gates
 
 
+def diagonalizing_clifford(pauli: str):
+    """Get the clifford gate list to diagonalize the Pauli operator.
+
+    Args:
+        pauli: The Pauli to diagonalize.
+
+    Returns:
+        A gate list for clifford.
+    """
+    reversed_pauli = pauli[::-1]
+    gates = []
+    gates_inverse = []
+
+    for i, pauli_i in enumerate(reversed_pauli):
+        if pauli_i == "Y":
+            gates.append(qeg.SdgGate(i))
+            gates_inverse.append(qeg.SGate(i))
+        if pauli_i in ["X", "Y"]:
+            gates.append(qeg.HGate(i))
+            gates_inverse.append(qeg.HGate(i))
+    gates_inverse = gates_inverse[::-1]
+    return gates, gates_inverse
+
+
+def cnot_chain(pauli: str):
+    """CX chain.
+
+        For example, for the Pauli with the label 'XYZIX'.
+
+                           ┌───┐
+            q_0: ──────────┤ X ├
+                           └─┬─┘
+            q_1: ────────────┼──
+                      ┌───┐  │
+            q_2: ─────┤ X ├──■──
+                 ┌───┐└─┬─┘
+            q_3: ┤ X ├──■───────
+                 └─┬─┘
+            q_4: ──■────────────
+
+        Args:
+            pauli: The Pauli for which to construct the CX chain.
+
+        Returns:
+            A gate list implementing the CX chain.
+        """
+
+    gates = []
+    control, target = None, None
+
+    # iterate over the Pauli's and add CNOTs
+    for i, pauli_i in enumerate(pauli):
+        i = len(pauli) - i - 1
+        if pauli_i != "I":
+            if control is None:
+                control = i
+            else:
+                target = i
+
+        if control is not None and target is not None:
+            gates.append(qeg.CXGate(control, target))
+            control = i
+            target = None
+
+    return gates, gates[::-1]
+
+
+def cnot_fountain(pauli: str):
+    """CX chain in the fountain shape.
+
+        For example, for the Pauli with the label 'XYZIX'.
+
+                 ┌───┐┌───┐┌───┐
+            q_0: ┤ X ├┤ X ├┤ X ├
+                 └─┬─┘└─┬─┘└─┬─┘
+            q_1: ──┼────┼────┼──
+                   │    │    │
+            q_2: ──■────┼────┼──
+                        │    │
+            q_3: ───────■────┼──
+                             │
+            q_4: ────────────■──
+
+        Args:
+            pauli: The Pauli for which to construct the CX chain.
+
+        Returns:
+            A gate list implementing the CX chain.
+        """
+
+    gates = []
+    control, target = None, None
+    for i, pauli_i in enumerate(pauli[::-1]):
+        if pauli_i != "I":
+            if target is None:
+                target = i
+            else:
+                control = i
+
+        if control is not None and target is not None:
+            gates.append(qeg.CXGate(control, target))
+            control = None
+
+    return gates, gates[::-1]
+
+
 class BaseEvolution(ABC):
     """Generate evolution circuit based on operators"""
 
@@ -68,9 +233,10 @@ def evol(self, pauli: str, time: float):
         if num_non_id == 0:
             pass
         elif num_non_id == 1:
-            pass
+            return single_qubit_evol(pauli, time)
         elif num_non_id == 2:
             return two_qubit_evol(pauli, time)
         else:
-            raise NotImplementedError(f"Pauli string {pauli} not yet supported")
+            # raise NotImplementedError(f"Pauli string {pauli} not yet supported")
+            return multi_qubit_evol(pauli, time)
         return []

tests/quafu/algorithms/integration_test.py

@@ -20,7 +20,7 @@
 from quafu import simulate
 from scipy.optimize import minimize
 import heapq
-
+import matplotlib.pyplot as plt
 from quafu.algorithms.ansatz import AlterLayeredAnsatz
 
 
@@ -60,12 +60,95 @@ def test_run(self):
         should have the highest probability
         """
         num_layers = 2
-        hamlitonian = Hamiltonian.from_pauli_list(
+        print("The test for ansatz.")
+
+        # test the zero qubit evolution
+        hamiltonian__ = Hamiltonian.from_pauli_list(
+            [("IIIII", 1), ("IIIII", 1), ("IIIII", 1), ("IIIII", 1)]
+        )
+        ansatz__ = QAOAAnsatz(hamiltonian__, num_layers=num_layers)
+        ansatz__.draw_circuit()
+
+        # test the single qubit evolution
+        hamiltonian_x = Hamiltonian.from_pauli_list(
+            [("IIIIX", 1), ("IIIXI", 1), ("IXIII", 1), ("XIIII", 1)]
+        )
+        ansatz_x = QAOAAnsatz(hamiltonian_x, num_layers=num_layers)
+        ansatz_x.draw_circuit()
+        hamiltonian_y = Hamiltonian.from_pauli_list(
+            [("IIIIY", 1), ("IIIYI", 1), ("IYIII", 1), ("YIIII", 1)]
+        )
+        ansatz_y = QAOAAnsatz(hamiltonian_y, num_layers=num_layers)
+        ansatz_y.draw_circuit()
+        hamiltonian_z = Hamiltonian.from_pauli_list(
+            [("IIIIZ", 1), ("IIIZI", 1), ("IZIII", 1), ("ZIIII", 1)]
+        )
+        ansatz_z = QAOAAnsatz(hamiltonian_z, num_layers=num_layers)
+        ansatz_z.draw_circuit()
+
+        # test the two qubits evolution
+        hamiltonian_xx = Hamiltonian.from_pauli_list(
+            [("IIIXX", 1), ("IIXIX", 1), ("IXIIX", 1), ("XIIIX", 1)]
+        )
+        ansatz_xx = QAOAAnsatz(hamiltonian_xx, num_layers=num_layers)
+        ansatz_xx.draw_circuit()
+        hamiltonian_xy = Hamiltonian.from_pauli_list(
+            [("IIIYX", 1), ("IIYIX", 1), ("IYIIX", 1), ("YIIIX", 1)]
+        )
+        ansatz_xy = QAOAAnsatz(hamiltonian_xy, num_layers=num_layers)
+        ansatz_xy.draw_circuit()
+        hamiltonian_xz = Hamiltonian.from_pauli_list(
+            [("IIIXZ", 1), ("IIZIX", 1), ("IZIIX", 1), ("ZIIIX", 1)]
+        )
+        ansatz_xz = QAOAAnsatz(hamiltonian_xz, num_layers=num_layers)
+        ansatz_xz.draw_circuit()
+        hamiltonian_yx = Hamiltonian.from_pauli_list(
+            [("IIIXY", 1), ("IIXIY", 1), ("IXIIY", 1), ("XIIIY", 1)]
+        )
+        ansatz_yx = QAOAAnsatz(hamiltonian_yx, num_layers=num_layers)
+        ansatz_yx.draw_circuit()
+        hamiltonian_yy = Hamiltonian.from_pauli_list(
+            [("IIIYY", 1), ("IIYIY", 1), ("IYIIY", 1), ("YIIIY", 1)]
+        )
+        ansatz_yy = QAOAAnsatz(hamiltonian_yy, num_layers=num_layers)
+        ansatz_yy.draw_circuit()
+        hamiltonian_yz = Hamiltonian.from_pauli_list(
+            [("IIIZY", 1), ("IIZIY", 1), ("IZIIY", 1), ("ZIIIY", 1)]
+        )
+        ansatz_yz = QAOAAnsatz(hamiltonian_yz, num_layers=num_layers)
+        ansatz_yz.draw_circuit()
+        hamiltonian_zx = Hamiltonian.from_pauli_list(
+            [("IIIXZ", 1), ("IIXIZ", 1), ("IXIIZ", 1), ("XIIIZ", 1)]
+        )
+        ansatz_zx = QAOAAnsatz(hamiltonian_zx, num_layers=num_layers)
+        ansatz_zx.draw_circuit()
+        hamiltonian_zy = Hamiltonian.from_pauli_list(
+            [("IIIYZ", 1), ("IIYIZ", 1), ("IYIIZ", 1), ("YIIIZ", 1)]
+        )
+        ansatz_zy = QAOAAnsatz(hamiltonian_zy, num_layers=num_layers)
+        ansatz_zy.draw_circuit()
+        hamiltonian_zz = Hamiltonian.from_pauli_list(
+            [("IIIZZ", 1), ("IIZIZ", 1), ("IZIIZ", 1), ("ZIIIZ", 1)]
+        )
+        ansatz_zz = QAOAAnsatz(hamiltonian_zz, num_layers=num_layers)
+        ansatz_zz.draw_circuit()
+
+        # test the multiple qubits evolution
+        hamiltonian_multi = Hamiltonian.from_pauli_list(
+            [("XYZIX", 1), ("XYIZX", 1), ("XIYZX", 1), ("IXYZX", 1)]
+        )
+        ansatz_multi = QAOAAnsatz(hamiltonian_multi, num_layers=num_layers)
+        ansatz_multi.draw_circuit()
+        # ansatz_multi.plot_circuit(title='MULTI QUBITS')
+        # plt.show()
+
+        hamiltonian = Hamiltonian.from_pauli_list(
             [("IIIZZ", 1), ("IIZIZ", 1), ("IZIIZ", 1), ("ZIIIZ", 1)]
         )
         ref_mat = np.load("tests/quafu/algorithms/data/qaoa_hamiltonian.npy")
-        assert np.array_equal(ref_mat, hamlitonian.get_matrix())
-        ansatz = QAOAAnsatz(hamlitonian, num_layers=num_layers)
+        # ref_mat = np.load("data/qaoa_hamiltonian.npy")
+        assert np.array_equal(ref_mat, hamiltonian.get_matrix())
+        ansatz = QAOAAnsatz(hamiltonian, num_layers=num_layers)
         ansatz.draw_circuit()
 
         def cost_func(params, ham, estimator: Estimator):
@@ -75,7 +158,7 @@ def cost_func(params, ham, estimator: Estimator):
         est = Estimator(ansatz)
         # params = 2 * np.pi * np.random.rand(num_layers * 2)
         params = 2 * np.pi * np.random.rand(ansatz.num_parameters)
-        res = minimize(cost_func, params, args=(hamlitonian, est), method="COBYLA")
+        res = minimize(cost_func, params, args=(hamiltonian, est), method="COBYLA")
         print(res)
         ansatz.measure(list(range(5)))
         ansatz.draw_circuit()