Merge pull request #591 from SRI-International/compiler-opt-changes-h…

…amsim

Compiler opt changes hamsim

hamlib/qiskit/hamlib_simulation_benchmark.py

@@ -38,6 +38,8 @@
 
 verbose = False
 
+# contains the correct bitstring for a random pauli circuit
+bitstring_dict = {}
 
 # Creates a key for distribution of initial state for method = 3.
 def key_from_initial_state(num_qubits, num_shots, init_state, random_pauli_flag):
@@ -150,19 +152,21 @@ def analyze_and_print_result(
 
         if verbose:
             print(f"... exact computation time = {round((time.time() - ts), 3)} sec")
-
-    # for method 3, compute expected distribution from the initial state
-    elif method == 3:
-
-        # check simple distribution if not inserting random Paulis
-        if not random_pauli_flag:
-            correct_dist = key_from_initial_state(num_qubits, num_shots, init_state, random_pauli_flag)
-
-        # if using random paulis, distribution is based on all ones in initial state
-        # random_pauli_flag is (for now) set to always return bitstring of 1's
-        else:
-            correct_dist = {"1" * num_qubits: num_shots}
 
+    # for method 3, compute expected distribution from the initial state
+    elif method == 3: 
+
+        # check simple distribution if not inserting random Paulis 
+        if not random_pauli_flag: 
+            correct_dist = key_from_initial_state(
+                num_qubits, num_shots, init_state, random_pauli_flag
+            )
+
+        # if using random paulis, a potentially random bitstring is collected from circuit generation
+        else: 
+            global bitstring_dict
+            correct_bitstring = bitstring_dict[qc.name]
+            correct_dist = {correct_bitstring: num_shots}
 
     else:
         raise ValueError("Method is not 1 or 2 or 3, or hamiltonian is not valid.")
@@ -220,7 +224,9 @@ def print_top_measurements(label, counts, top_n):
 def run(min_qubits: int = 2, max_qubits: int = 8, max_circuits: int = 1,
         skip_qubits: int = 1, num_shots: int = 100,
         hamiltonian: str = "TFIM", method: int = 1,
-        random_pauli_flag: bool = False, init_state: str = None,
+        random_pauli_flag: bool = False, 
+        random_init_flag: bool = False, 
+        init_state: str = None,
         K: int = None, t: float = None,
         backend_id: str = None, provider_backend = None,
         hub: str = "ibm-q", group: str = "open", project: str = "main", exec_options = None,
@@ -325,19 +331,26 @@ def execution_handler(qc, result, num_qubits, type, num_shots):
 
         #######################################################################
 
-        # Loop over only the same circuit, executing it num_circuits times
+        # in the case of random paulis, method = 3: loop over multiple random pauli circuits
+        # otherwise, loop over the same circuit, executing it num_circuits times 
         for circuit_id in range(num_circuits):
+
             ts = time.time()
-
-            # create the HamLibSimulation kernel and the associated Hamiltonian operator
-            qc, ham_op = HamiltonianSimulation(
-                    num_qubits,
-                    hamiltonian=hamiltonian, 
-                    K=K, t=t,
-                    init_state=init_state,
-                    method = method,
-                    random_pauli_flag=random_pauli_flag
-                    )
+
+            #used to store random pauli correct bitstrings
+            global bitstring_dict
+
+            # create the HamLibSimulation kernel, random pauli bitstring, and the associated Hamiltonian operator
+            qc, bs, ham_op = HamiltonianSimulation(
+                num_qubits, 
+                K=K, t=t,
+                hamiltonian=hamiltonian, 
+                init_state=init_state,
+                method = method, 
+                random_pauli_flag=random_pauli_flag, 
+                random_init_flag=random_init_flag)
+
+            bitstring_dict[qc.name] = bs
 
             metrics.store_metric(num_qubits, circuit_id, 'create_time', time.time() - ts)
 
@@ -390,6 +403,7 @@ def get_args():
     parser.add_argument("--global_rinst", "-param_rinst", default=None, help="paramater rinst")
     parser.add_argument("--num_steps", "-steps", default=None, help="Number of Trotter steps", type=int)
     parser.add_argument("--time", "-time", default=None, help="Time of evolution", type=float)
+    parser.add_argument("--random_init_flag", "-rani", action="store_true", help="random pauli flag")
     return parser.parse_args()
 
 # if main, execute method
@@ -420,6 +434,7 @@ def get_args():
         hamiltonian=args.hamiltonian,
         method=args.method,
         random_pauli_flag=args.random_pauli_flag,
+        random_init_flag=args.random_init_flag,
         init_state = args.init_state,
         K = args.num_steps,
         t = args.time,

hamlib/qiskit/hamlib_simulation_kernel.py

@@ -228,6 +228,7 @@ def create_circuit(
     method: int = 1,
     init_state: str = None,
     random_pauli_flag: bool = False,
+    random_init_flag: bool = False 
 ):
     """
     Create a quantum circuit based on the Hamiltonian data from an HDF5 file.
@@ -292,6 +293,8 @@ def create_circuit(
 
         # Append K Trotter steps of inverse, if method 3
         inv = None
+        bitstring = None
+
         if method == 3: 
 
             # if not adding random Paulis, just create simple inverse Trotter steps
@@ -305,17 +308,23 @@ def create_circuit(
             # if adding Paulis, do that here, with code from pyGSTi
             else:
                 from pygsti_mirror import convert_to_mirror_circuit
-                circuit, _ = convert_to_mirror_circuit(circuit_without_initial_state)
-
+
+                # if random init flag state is set, then discard the inital state input and use a completely (harr) random one
+                if random_init_flag:
+                    circuit, bitstring = convert_to_mirror_circuit(circuit_without_initial_state, random_pauli = True, init_state=None)
+                else: 
+                    init_state = initial_state(n_spins, init_state)
+                    circuit, bitstring = convert_to_mirror_circuit(circuit_without_initial_state, random_pauli = True, init_state=init_state)
+
         # convert_to_mirror_circuit adds its own measurement gates    
         if not random_pauli_flag:
             circuit.measure_all()
 
-        return circuit, ham_op, evo
-        
+        return circuit, bitstring, ham_op, evo
+
     else:
         # print(f"Dataset not available for n_spins = {n_spins}.")
-        return None, None, None
+        return None, None, None, None
 
 
 ############### Initial Circuit Definition
@@ -357,7 +366,8 @@ def HamiltonianSimulation(
             K: int = 5, t: float = 1.0,
             init_state = None,
             method: int = 1,
-            random_pauli_flag = False
+            random_pauli_flag = False,
+            random_init_flag = False
         ) -> QuantumCircuit:
     """
     Construct a Qiskit circuit for Hamiltonian simulation.
@@ -382,15 +392,15 @@ def HamiltonianSimulation(
     qc = QuantumCircuit(qr, cr, name=f"hamsim-{num_qubits}-{secret_int}")
 
     hamiltonian = hamiltonian.strip().lower()
-
-    # create the quantum circuit for this Hamiltonian, along with the operator and trotter evolution circuit
-    qc, ham_op, evo = create_circuit(
+    # create the quantum circuit for this Hamiltonian, along with the correct pauli bstring, the operator and trotter evolution circuit
+    qc, bitstring, ham_op, evo = create_circuit(
         n_spins=n_spins,
         time=t,
         method=method,
         init_state=init_state,
         num_trotter_steps=K,
         random_pauli_flag=random_pauli_flag,
+        random_init_flag=random_init_flag
         )
 
     # Save smaller circuit example for display
@@ -403,10 +413,12 @@ def HamiltonianSimulation(
 
     # Collapse the sub-circuits used in this benchmark (for Qiskit)
     qc2 = qc.decompose().decompose()
-
-    # return both the circuit created and the Hamiltonian operator
-    return qc2, ham_op
-
+
+    # return both the circuit created, the bitstring, and the Hamiltonian operator
+    # if random_pauli_flag is false or method isn't 3, bitstring will be None
+    return qc2, bitstring, ham_op
+
+
 
 ############### Circuit Drawer
 
@@ -442,4 +454,3 @@ def kernel_draw(hamiltonian: str = "hamlib", method: int = 1):
         print("  ... circuit too large!")
 
 
-

hamlib/qiskit/pygsti_mirror.py

@@ -3,18 +3,20 @@
 import numpy as np
 import scipy as sp
 from qiskit import transpile, QuantumCircuit
+import re
 
 # Below code created by Timothy Proctor, July 16 2024, using code developed by various members of Sandia's QPL.
+# Code has been edited to accept random pauli bitstrings by Sonny Rappaport, July 23 2024.
+# Additional functions below also by Sonny based on Tim's code. 
 
-def sample_mirror_circuits(c, num_mcs=100, randomized_state_preparation=True, rand_state=None):
+def sample_mirror_circuits(c, num_mcs=100, randomized_state_preparation=True, rand_state = None, random_pauli = True):
     if rand_state is None:
         rand_state = np.random.RandomState()
-        
+
     mirror_circuits = []
     mirror_circuits_bitstrings = []
     for j in range(num_mcs):
-        mc, bs = central_pauli_mirror_circuit(c, randomized_state_preparation=randomized_state_preparation,
-                                              rand_state=rand_state)
+        mc, bs = central_pauli_mirror_circuit(c, randomized_state_preparation=randomized_state_preparation,random_pauli = True, rand_state = rand_state)
         mirror_circuits.append(mc)
         mirror_circuits_bitstrings.append(bs)
 
@@ -36,7 +38,7 @@ def sample_reference_circuits(qubits, num_ref=100, randomized_state_preparation=
 
     return ref_circuits, ref_circuits_bitstrings
 
-def central_pauli_mirror_circuit(circ, randomized_state_preparation=True, rand_state=None):
+def central_pauli_mirror_circuit(circ, randomized_state_preparation=True, random_pauli = True, rand_state = None):
     if rand_state is None:
         rand_state = np.random.RandomState()
 
@@ -50,10 +52,13 @@ def central_pauli_mirror_circuit(circ, randomized_state_preparation=True, rand_s
     n = circuit_to_mirror.width
     d = circuit_to_mirror.depth
 
-    # SONNY'S EDIT TO TIM CODE: MAKE THE CENTRAL PAULI (for now) NOT RANDOM
-    # Note that we could rework the code to pass in a desired bitstring here later, rather than randomly generating one here.
-    #central_pauli = 2 * rand_state.randint(0, 2, 2*n)
-    central_pauli = 2 * np.ones(2*n, dtype = np.int8 )  
+    # Sonny's edit here: Rather than generating central_pauli randomly, now we can use a non-random one if desired. 
+    if random_pauli is True:
+        rand_state = np.random.RandomState()
+        central_pauli = 2 * np.random.randint(0, 2, 2*n) #old code 
+    else: 
+        central_pauli = 2 * np.ones(2*n, dtype = np.int8)  
+
     central_pauli_layer = pauli_vector_to_u3_layer(central_pauli, qubits)
     q = central_pauli.copy()
 
@@ -125,8 +130,8 @@ def haar_random_u3(q, rand_state=None):
     if rand_state is None:
         rand_state = np.random.RandomState()
 
-    a, b = 2 * np.pi * rand_state.rand(2)
-    theta = mod_2pi(2 * np.arcsin(np.sqrt(rand_state.rand(1)))[0])
+    a, b = 2 * np.pi * np.random.rand(2)
+    theta = mod_2pi(2 * np.arcsin(np.sqrt(np.random.rand(1)))[0])
     phi = mod_2pi(a - b + np.pi)
     lamb = mod_2pi(-1 * (a + b + np.pi))
     return pygsti.baseobjs.Label('Gu3', q, args=(theta, phi, lamb))
@@ -325,23 +330,36 @@ def pygsti_string(qc):
 
     return pygstr.strip()
 
-def convert_to_mirror_circuit(qc):
+def convert_to_mirror_circuit(qc, random_pauli, init_state):
     """
     Takes an arbitrary quantum circuit, transpiles it to CX and u3, and eventually returns a mirror circuit version of the quantum circuit. 
 
     The mirror circuit will have a random initial state in addition to random paulis. 
     """
 
+    # if no initial state is provided, use a ranmdom one
+    random_init_flag = True  
+
+    if init_state: 
+        random_init_flag = False
+        init_state.append(qc, range(qc.num_qubits))
+        qc = init_state
+
     qc = transpile(qc, basis_gates=["cx","u3"])
 
     pygstr = pygsti_string(qc)
 
     pygsti_circuit = pygsti.circuits.Circuit(pygstr)
 
-    mcs, bss = sample_mirror_circuits(pygsti_circuit, num_mcs=1)
-    print(bss)
+    mcs, bss = sample_mirror_circuits(pygsti_circuit, num_mcs=1, random_pauli = random_pauli, randomized_state_preparation=random_init_flag)
+
+    qasm_string = mcs[0].convert_to_openqasm()
+
+    # use regex magic to remove lines that begin with Delay gates
+    delays_removed_string = re.sub(r'^delay.*\n?', '', qasm_string, flags=re.MULTILINE) 
 
-    qiskit_circuit = QuantumCircuit.from_qasm_str(mcs[0].convert_to_openqasm())
+    qiskit_circuit = QuantumCircuit.from_qasm_str(delays_removed_string)
 
     # this circuit is made out of u3 and cx gates by pygsti default
-    return qiskit_circuit , bss[0]
+    # the bitstring is reversed to account for qiskit ordering
+    return qiskit_circuit , bss[0][::-1]