Random Pauli applied for mirror circuit

hamiltonian-simulation/qiskit/hamiltonian_simulation_benchmark.py

@@ -46,7 +46,7 @@
 
 #def analyze_and_print_result(qc: QuantumCircuit, result, num_qubits: int,
 def analyze_and_print_result(qc, result, num_qubits: int,
-            type: str, num_shots: int, hamiltonian: str, method: int) -> tuple:
+            type: str, num_shots: int, hamiltonian: str, method: int, random_pauli_flag: bool) -> tuple:
     """
     Analyze and print the measured results. Compute the quality of the result based on operator expectation for each state.
 
@@ -63,7 +63,6 @@ def analyze_and_print_result(qc, result, num_qubits: int,
         tuple: Counts and fidelity.
     """
     counts = result.get_counts(qc)
-
     if verbose:
         print(f"For type {type} measured: {counts}")
 
@@ -79,9 +78,16 @@ def analyze_and_print_result(qc, result, num_qubits: int,
     elif method == 2 and hamiltonian == "tfim":
         correct_dist = precalculated_data[f"Exact TFIM - Qubits{num_qubits}"]
     elif method == 3 and hamiltonian == "heisenberg":
-        correct_dist = {''.join(['0' if i % 2 == 0 else '1' for i in range(num_qubits)]) if num_qubits % 2 != 0 else ''.join(['1' if i % 2 == 0 else '0' for i in range(num_qubits)]):num_shots}
+        if random_pauli_flag == True:
+            correct_dist = {''.join(['0' if i % 2 == 0 else '1' for i in range(num_qubits)]) if num_qubits % 2 != 0 else ''.join(['1' if i % 2 == 0 else '0' for i in range(num_qubits)]):num_shots}
+        else:
+            correct_dist = {''.join(['1' if i % 2 == 0 else '0' for i in range(num_qubits)]) if num_qubits % 2 != 0 else ''.join(['0' if i % 2 == 0 else '1' for i in range(num_qubits)]):num_shots}    
     elif method == 3 and hamiltonian == "tfim":
-        correct_dist = {'0' * num_qubits: num_shots // 2 + num_shots % 2, '1' * num_qubits: num_shots // 2}
+        if random_pauli_flag == True:
+            correct_dist = {'0' * num_qubits: num_shots // 2 + num_shots % 2, '1' * num_qubits: num_shots // 2}
+        else:
+            correct_dist = {'0' * num_qubits: num_shots // 2 + num_shots % 2, '1' * num_qubits: num_shots // 2}
+
     else:
         raise ValueError("Method is not 1 or 2 or 3, or hamiltonian is not tfim or heisenberg.")
 
@@ -90,15 +96,17 @@ def analyze_and_print_result(qc, result, num_qubits: int,
 
     # Use polarization fidelity rescaling
     fidelity = metrics.polarization_fidelity(counts, correct_dist)
-    print(counts,correct_dist)
+    print("count", counts)
+    print("exp:", correct_dist)
     return counts, fidelity
 
+
 ############### Benchmark Loop
 
 def run(min_qubits: int = 2, max_qubits: int = 8, max_circuits: int = 3,
         skip_qubits: int = 1, num_shots: int = 100,
         hamiltonian: str = "heisenberg", method: int = 1,
-        use_XX_YY_ZZ_gates: bool = False,
+        use_XX_YY_ZZ_gates: bool = False, random_pauli_flag = True,
         backend_id: str = None, provider_backend = None,
         hub: str = "ibm-q", group: str = "open", project: str = "main", exec_options = None,
         context = None, api = None):
@@ -148,7 +156,7 @@ def run(min_qubits: int = 2, max_qubits: int = 8, max_circuits: int = 3,
     def execution_handler(qc, result, num_qubits, type, num_shots):
         # Determine fidelity of result set
         num_qubits = int(num_qubits)
-        counts, expectation_a = analyze_and_print_result(qc, result, num_qubits, type, num_shots, hamiltonian, method)
+        counts, expectation_a = analyze_and_print_result(qc, result, num_qubits, type, num_shots, hamiltonian, method, random_pauli_flag)
         metrics.store_metric(num_qubits, type, 'fidelity', expectation_a)
 
     # Initialize execution module using the execution result handler above and specified backend_id
@@ -189,7 +197,7 @@ def execution_handler(qc, result, num_qubits, type, num_shots):
                     hamiltonian=hamiltonian,
                     w=w, hx = hx, hz = hz, 
                     use_XX_YY_ZZ_gates = use_XX_YY_ZZ_gates,
-                    method = method)
+                    method = method, random_pauli_flag = random_pauli_flag)
 
             metrics.store_metric(num_qubits, circuit_id, 'create_time', time.time() - ts)
             qc.draw()
@@ -206,7 +214,7 @@ def execution_handler(qc, result, num_qubits, type, num_shots):
     ##########
 
     # draw a sample circuit
-    kernel_draw(hamiltonian, use_XX_YY_ZZ_gates, method)
+    kernel_draw(hamiltonian, use_XX_YY_ZZ_gates, method, random_pauli_flag)
 
     # Plot metrics for all circuit sizes
     metrics.plot_metrics(f"Benchmark Results - {benchmark_name} - Qiskit")
@@ -233,6 +241,7 @@ def get_args():
     #parser.add_argument("--theta", default=0.0, help="Input Theta Value", type=float)
     parser.add_argument("--nonoise", "-non", action="store_true", help="Use Noiseless Simulator")
     parser.add_argument("--verbose", "-v", action="store_true", help="Verbose")
+    parser.add_argument("--random_pauli_flag", "-ranp", action="store_true", help="random pauli flag")
     return parser.parse_args()
 
 # if main, execute method
@@ -255,7 +264,8 @@ def get_args():
         num_shots=args.num_shots,
         hamiltonian=args.hamiltonian,
         method=args.method,
-        use_XX_YY_ZZ_gates = args.use_XX_YY_ZZ_gates,
+        random_pauli_flag=args.random_pauli_flag,
+        use_XX_YY_ZZ_gates =args.use_XX_YY_ZZ_gates,
         #theta=args.theta,
         backend_id=args.backend_id,
         exec_options = {"noise_model" : None} if args.nonoise else {},

hamiltonian-simulation/qiskit/hamiltonian_simulation_kernel.py

@@ -58,7 +58,6 @@ def initial_state(n_spins: int, initial_state: str = "checker") -> QuantumCircui
 ############## Heisenberg Circuit
 def Heisenberg(n_spins: int, K: int, t: float, tau: float, w: float, h_x: list[float], h_z: list[float],
             use_XX_YY_ZZ_gates: bool = False) -> QuantumCircuit:
-
     qr = QuantumRegister(n_spins)
     qc = QuantumCircuit(qr, name = "Heisenberg")
     # Loop over each Trotter step, adding gates to the circuit defining the Hamiltonian
@@ -94,7 +93,7 @@ def Heisenberg(n_spins: int, K: int, t: float, tau: float, w: float, h_x: list[f
 
     return qc
 
-########### TFIM circuit
+########### TFIM hamiltonian circuit
 def tfim(n_spins: int, K: int, tau: float, use_XX_YY_ZZ_gates: bool)-> QuantumCircuit:
     h = 0.2  # Strength of transverse field
     qr = QuantumRegister(n_spins)
@@ -137,24 +136,80 @@ def ResultantPauli(n_spins)-> QuantumCircuit:
         qc.barrier()
     return qc
 
-############ Quasi Inverse Heisenberg    
+############ Quasi Inverse Heisenberg 
+########### ~H P H = R ==> ~H = R H' P'  ; ~H is QuasiHamiltonian, P is Random Pauli, H is Hamiltonian, R is resultant circuit that appends on the initial state
 def QuasiHamiltonian(hamiltonian_circuit, random_pauli_oracle, res_pauli, n_spins)-> QuantumCircuit:
                 qr = QuantumRegister(n_spins)
                 qc = QuantumCircuit(qr, name = "QuasiHamiltonian")
                 hamiltonian_circuit_inverse = hamiltonian_circuit.inverse()
                 random_pauli_oracle_inverse = random_pauli_oracle.inverse()
 
-                qc.append(random_pauli_oracle_inverse,qr)
+                qc.append(random_pauli_oracle_inverse,qr)             
                 qc.append(hamiltonian_circuit_inverse ,qr)
 
-                qc.append(res_pauli,qr)
+                qc.append(res_pauli,qr) 
 
                 return qc
 
+#Inverse of Heisenberg model. mirror gates are applied.
+def HeisenbergInverse(n_spins: int, K: int, t: float, tau: float, w: float, h_x: list[float], h_z: list[float],
+            use_XX_YY_ZZ_gates: bool = False) -> QuantumCircuit:
+
+    qr = QuantumRegister(n_spins)
+    qc = QuantumCircuit(qr, name = "HeisenbergInverse")
+    # Add mirror gates for negative time simulation
+    for k in range(K): 
+        # Basic implementation of exp(-i * t * (XX + YY + ZZ)):
+        if use_XX_YY_ZZ_gates:
+            # regular inverse of XX + YY + ZZ operators on each pair of quibts in linear chain
+            # XX operator on each pair of qubits in linear chain
+            for j in range(2):
+                for i in range(j%2, n_spins - 1, 2):
+                    qc.append(zz_gate_mirror(tau).to_instruction(), [qr[i], qr[(i + 1) % n_spins]])
+
+            # YY operator on each pair of qubits in linear chain
+            for j in range(2):
+                for i in range(j%2, n_spins - 1, 2):
+                    qc.append(yy_gate_mirror(tau).to_instruction(), [qr[i], qr[(i + 1) % n_spins]])
+
+            # ZZ operation on each pair of qubits in linear chain
+            for j in range(2):
+                for i in range(j%2, n_spins - 1, 2):
+                    qc.append(xx_gate_mirror(tau).to_instruction(), [qr[i], qr[(i + 1) % n_spins]])
+
+        else:
+            # optimized Inverse of XX + YY + ZZ operator on each pair of qubits in linear chain
+            for j in range(2):
+                for i in range(j % 2, n_spins - 1, 2):
+                    qc.append(xxyyzz_opt_gate_mirror(tau).to_instruction(), [qr[i], qr[(i + 1) % n_spins]])
+        qc.barrier()
+
+        # the Pauli spin vector product
+        [qc.rz(-2 * tau * w * h_z[i], qr[i]) for i in range(n_spins)]
+        [qc.rx(-2 * tau * w * h_x[i], qr[i]) for i in range(n_spins)]
+        qc.barrier()
+    return qc
+
+#########Inverse of tfim hamiltonian
+def tfimInverse(n_spins: int, K: int, tau: float, use_XX_YY_ZZ_gates: bool)-> QuantumCircuit:
+    h = 0.2
+    qr = QuantumRegister(n_spins)
+    qc = QuantumCircuit(qr, name = "tfimInverse")
+    for k in range(K):
+        for j in range(2):
+            for i in range(j % 2, n_spins - 1, 2):
+                qc.append(zz_gate_mirror(tau).to_instruction(), [qr[i], qr[(i + 1) % n_spins]])
+        qc.barrier()
+        for i in range(n_spins):
+            qc.rx(-2 * tau * h, qr[i])
+        qc.barrier()
+    return qc
+
+
 def HamiltonianSimulation(n_spins: int, K: int, t: float,
             hamiltonian: str, w: float, hx: list[float], hz: list[float],
             use_XX_YY_ZZ_gates: bool = False,
-            method: int = 1) -> QuantumCircuit:
+            method: int = 1, random_pauli_flag: bool = True) -> QuantumCircuit:
     """
     Construct a Qiskit circuit for Hamiltonian simulation.
 
@@ -187,23 +242,29 @@ def HamiltonianSimulation(n_spins: int, K: int, t: float,
 
     hamiltonian = hamiltonian.strip().lower()
 
-
     if hamiltonian == "heisenberg": 
         init_state = "checkerboard"
         # apply initial state
-        qc.append(initial_state(n_spins, init_state), qr)
+        qc.append(initial_state(n_spins, init_state), qr) # append the initial circuit on the quantum circuit.
         qc.barrier()
-        heisenberg_circuit = Heisenberg(n_spins, K, t, tau, w, h_x, h_z, use_XX_YY_ZZ_gates)
+        heisenberg_circuit = Heisenberg(n_spins, K, t, tau, w, h_x, h_z, use_XX_YY_ZZ_gates) #append the heisenberg on the circuit
         qc.append(heisenberg_circuit, qr)
         qc.barrier()
 
         if (method == 3):
-            random_pauli_oracle = create_random_paulis(n_spins)
-            qc.append(random_pauli_oracle, qr)
-            qc.barrier()
-            res_pauli = ResultantPauli(n_spins)
-            Quasi_heisenberg = QuasiHamiltonian(heisenberg_circuit, random_pauli_oracle, res_pauli, n_spins)
-            qc.append(Quasi_heisenberg, qr)
+            if random_pauli_flag:
+                random_pauli_oracle = create_random_paulis(n_spins)
+                qc.append(random_pauli_oracle, qr)    #append the random pauli on the circuit
+                qc.barrier()
+                res_pauli = ResultantPauli(n_spins) # create a resultant pauli that we want to apply to initial state.
+                Quasi_heisenberg = QuasiHamiltonian(heisenberg_circuit, random_pauli_oracle, res_pauli, n_spins) # create a QuasiHamiltonian. 
+                qc.append(Quasi_heisenberg, qr)    #append the Quasi Hamiltonian on the circuit
+
+            else:
+                #if random_pauli_flag is False, just use traditional mirror circuit, i.e. Apply Inverse of Hamiltonian to the Hamiltonian to give Inverse.
+                heisenberg_inverse_circuit = HeisenbergInverse(n_spins, K, t, tau, w, h_x, h_z, use_XX_YY_ZZ_gates)
+                qc.append(heisenberg_inverse_circuit, qr)
+                qc.barrier()
 
 
     elif hamiltonian == "tfim":
@@ -217,12 +278,18 @@ def HamiltonianSimulation(n_spins: int, K: int, t: float,
         qc.barrier()
 
         if (method == 3):
-            random_pauli_oracle = create_random_paulis(n_spins)
-            qc.append(random_pauli_oracle, qr)
-            qc.barrier()
-            res_pauli = ResultantPauli(n_spins)
-            Quasi_tfim = QuasiHamiltonian(tfim_circuit, random_pauli_oracle, res_pauli, n_spins)
-            qc.append(Quasi_tfim, qr)
+            if random_pauli_flag == True:
+                random_pauli_oracle = create_random_paulis(n_spins)
+                qc.append(random_pauli_oracle, qr)  #append the random pauli on the circuit
+                qc.barrier()
+                res_pauli = ResultantPauli(n_spins) # create a resultant pauli that we want to apply to initial state.
+                Quasi_tfim = QuasiHamiltonian(tfim_circuit, random_pauli_oracle, res_pauli, n_spins) # create a QuasiHamiltonian.
+                qc.append(Quasi_tfim, qr) #append the Quasi Hamiltonian on the circuit
+            else:
+                #if random_pauli_flag is False, just use traditional mirror circuit, i.e. Apply Inverse of Hamiltonian to the Hamiltonian to give Inverse.
+                tfim_inverse_circuit = tfimInverse(n_spins, K, tau, use_XX_YY_ZZ_gates)
+                qc.append(tfim_inverse_circuit, qr)
+                qc.barrier()
 
     else:
         raise ValueError("Invalid Hamiltonian specification.")
@@ -461,7 +528,7 @@ def xxyyzz_opt_gate_mirror(tau: float) -> QuantumCircuit:
 ############### BV Circuit Drawer
 
 # Draw the circuits of this benchmark program
-def kernel_draw(hamiltonian: str = "heisenberg", use_XX_YY_ZZ_gates: bool = False, method: int = 1):
+def kernel_draw(hamiltonian: str = "heisenberg", use_XX_YY_ZZ_gates: bool = False, method: int = 1,random_pauli_flag: bool = True):
 
     # Print a sample circuit
     print("Sample Circuit:")