further simulator upgrades + bugfixes

src/qrisp/core/quantum_array.py

@@ -431,7 +431,7 @@ def __array_function__(self, func, types, args, kwargs):
     def get_measurement(
         self,
         backend=None,
-        shots=10000,
+        shots=100000,
         compile=True,
         compilation_kwargs={},
         subs_dic={},

src/qrisp/core/quantum_variable.py

@@ -858,7 +858,7 @@ def get_measurement(
         self,
         plot=False,
         backend=None,
-        shots=10000,
+        shots=100000,
         compile=True,
         compilation_kwargs={},
         subs_dic={},

src/qrisp/misc/utility.py

@@ -32,6 +32,7 @@ def bin_rep(n, bits):
             str(n) + " can't be represented as a " + str(bits) + " bit number"
         )
 
+    return bin(n)[2:].zfill(bits)
     zero_string = "".join(["0" for k in range(bits)])
     return (zero_string + bin(n)[2:])[-bits:]
 
@@ -639,7 +640,7 @@ def find_qs(args):
 
 
 # Function to measure multiple quantum variables at once to assess their entanglement
-def multi_measurement(qv_list, shots=10000, backend=None):
+def multi_measurement(qv_list, shots=100000, backend=None):
     """
     This functions facilitates the measurement of multiple QuantumVariables at the same
     time. This can be used if the entanglement structure between several
@@ -711,6 +712,7 @@ def multi_measurement(qv_list, shots=10000, backend=None):
     # compiled_qc = qv_list[0].qs.copy()
     # Add classical registers for the measurement results to be stored in
     cl_reg_list = []
+
 
     for var in qv_list[::-1]:
         cl_reg = []
@@ -1274,7 +1276,7 @@ def check_if_fresh(qubits, qs, ignore_q_envs = True):
     return True
 
 
-def get_measurement_from_qc(qc, qubits, backend, shots=10000):
+def get_measurement_from_qc(qc, qubits, backend, shots=100000):
     # Add classical registers for the measurement results to be stored in
     cl = []
     for i in range(len(qubits)):
@@ -1302,6 +1304,7 @@ def get_measurement_from_qc(qc, qubits, backend, shots=10000):
             new_counts_dic[new_key] = counts[key]
 
     counts = new_counts_dic
+
     # Plot result (if needed)
 
     # Normalize counts

src/qrisp/qaoa/problems/maxCut.py

@@ -21,13 +21,27 @@
 from scipy.optimize import minimize
 from sympy import Symbol
 import itertools
+from numba import njit, prange
+
+@njit(cache = True)
+def maxcut_obj(x, edge_list):
+    cut = 0
+    for i, j in edge_list:
+        # the edge is cut
+        if ((x >> i) ^ (x >>j)) & 1:
+        # if x[i] != x[j]:                          
+            cut -= 1
+    return cut
+
+@njit(parallel = True, cache = True)
+def maxcut_energy(outcome_array, count_array, edge_list):
+
+    res_array = np.zeros(len(outcome_array))    
+    for i in prange(len(outcome_array)):
+        res_array[i] = maxcut_obj(outcome_array[i], edge_list)*count_array[i]
+
+    return np.sum(res_array)
 
-def maxcut_obj(x,G):
-        cut = 0
-        for i, j in G.edges():
-            if x[i] != x[j]:                        
-                cut -= 1    
-        return cut
 
 def create_maxcut_cl_cost_function(G):
     """
@@ -48,20 +62,25 @@ def create_maxcut_cl_cost_function(G):
     """    
     def cl_cost_function(counts):
 
-        def maxcut_obj(x, G):
-            cut = 0
-            for i, j in G.edges():
-                # the edge is cut
-                if x[i] != x[j]:                          
-                    cut -= 1
-            return cut
+        edge_list = np.array(list(G.edges()), dtype = np.uint32)
+
+        counts_keys = list(counts.keys())
 
-        energy = 0
-        for meas, meas_count in counts.items():
-            obj_for_meas = maxcut_obj(meas, G)
-            energy += obj_for_meas * meas_count
+        int_list = []
+        if not isinstance(counts_keys[0], str):
+
+            for c_array in counts_keys:
+                integer = int("".join([c for c in c_array]), 2)
+                int_list.append(integer)
+        else:
+            for c_str in counts_keys:
+                integer = int(c_str, 2)
+                int_list.append(integer)
+
+        counts_array = np.array(list(counts.values()))
+        outcome_array = np.array(int_list, dtype = np.uint32)
 
-        return energy
+        return maxcut_energy(outcome_array, counts_array, edge_list)
 
     return cl_cost_function
 

src/qrisp/qaoa/qaoa_problem.py

@@ -474,6 +474,8 @@ def run(self, qarg, depth, mes_kwargs = {}, max_iter = 50, init_type = "random")
         #bound_qc = self.train_circuit(qarg, depth)
         #opt_res = bound_qc(qarg).get_measurement(**mes_kwargs)
         #return opt_res
+        if not "shots" in mes_kwargs:
+            mes_kwargs["shots"] = 5000
 
         #res_sample = self.optimization_routine(qarg, compiled_qc, symbols , depth,  mes_kwargs, max_iter)
         res_sample = self.optimization_routine(qarg, depth, mes_kwargs, max_iter)

src/qrisp/qtypes/quantum_string.py

@@ -148,6 +148,7 @@ def __iadd__(self, other):
                 self[i + len(self) - len(other)] = other[i]
 
         elif isinstance(other, QuantumChar):
+            merge(self, other)
             self.resize((len(self) + 1,), refcheck=False)
             np.ndarray.__setitem__(self, len(self) - 1, other)
 

src/qrisp/simulator/bi_array_helper.py

@@ -212,23 +212,37 @@ def permuter(index, new_index, bit_partition, permuted_bit_partition, perm):
         # bit_range = extract_bit_range(index, bit_partition[i], bit_partition[i+1])
         # insert_bit_range(new_index, bit_range, permuted_bit_partition[perm[i]])
         # Inlined version more efficient:
-        new_index |= (
-            (
-                index
-                & (
-                    (int(1 << (bit_partition[i + 1] - bit_partition[i])) - 1)
-                    << bit_partition[i]
+        if isinstance(new_index, np.ndarray):
+            new_index |= (
+                (
+                    index
+                    & (
+                        (int(1 << (bit_partition[i + 1] - bit_partition[i])) - 1)
+                        << bit_partition[i]
+                    )
                 )
-            )
-            >> bit_partition[i]
-        ) << permuted_bit_partition[perm[i]]
+                >> bit_partition[i]
+            ) << permuted_bit_partition[perm[i]]
+        else:
+            for j in range(len(index)):
+                new_index[j] |= (
+                    (
+                        index[j]
+                        & (
+                            (int(1 << (bit_partition[i + 1] - bit_partition[i])) - 1)
+                            << int(bit_partition[i])
+                        )
+                    )
+                    >> int(bit_partition[i])
+                ) << int(permuted_bit_partition[perm[i]])
+
 
 
 def invert_perm(perm):
     return [perm.index(i) for i in range(len(perm))]
 
 
-def permute_axes(index, index_bit_permutation):
+def permute_axes(index, index_bit_permutation, jit = True):
     n = len(index_bit_permutation)
 
     index_bit_permutation = -np.array(index_bit_permutation) + n - 1
@@ -261,18 +275,21 @@ def permute_axes(index, index_bit_permutation):
         else:
             i += 1
 
-    # print(perm)
-    new_index = np.zeros(index.size, dtype=index.dtype)
+    if isinstance(index, np.ndarray):
+        new_index = np.zeros(index.size, dtype=index.dtype)
+        dtype = index.dtype
+    else:
+        new_index = [0 for _ in range(len(index))]
+        dtype = np.int64
 
     bit_partition = [sum(perm_shape[:i]) for i in range(len(perm_shape) + 1)]
 
-    # print(bit_partition)
     perm_shape_permuted = [perm_shape[i] for i in invert_perm(perm)]
     permuted_bit_partition = [
         sum(perm_shape_permuted[:i]) for i in range(len(perm_shape) + 1)
     ]
 
-    if len(perm) * index.size > 2**10:
+    if len(perm) * len(index) > 2**10 and jit:
         jitted_permuter(
             index,
             new_index,
@@ -281,12 +298,13 @@ def permute_axes(index, index_bit_permutation):
             np.array(perm),
         )
     else:
+
         permuter(
             index,
             new_index,
-            np.array(bit_partition, dtype = index.dtype),
-            np.array(permuted_bit_partition, dtype = index.dtype),
-            np.array(perm, dtype = index.dtype),
+            np.array(bit_partition, dtype = dtype),
+            np.array(permuted_bit_partition, dtype = dtype),
+            np.array(perm, dtype = dtype),
         )
 
     return new_index
@@ -314,6 +332,7 @@ def bi_array_moveaxis(data_array, index_perm, f_index_array):
     # return f_index_array
     return data_array[f_index_array]
 
+
 @njit(parallel = True, cache = True)
 def dense_measurement_brute(input_array, mes_amount, outcome_index):
 
@@ -330,6 +349,10 @@ def dense_measurement_brute(input_array, mes_amount, outcome_index):
     p_array = np.abs(p_array)
     max_p_array = np.max(p_array)
 
+    indices = np.nonzero(p_array > max_p_array*cutoff_ratio)[0]
+
+    return reshaped_array[indices,:], p_array[indices], indices
+
     new_arrays = []
     p_values = []
     outcome_indices = []
@@ -343,6 +366,7 @@ def dense_measurement_brute(input_array, mes_amount, outcome_index):
 
     return new_arrays, p_values, outcome_indices
 
+
 @njit(nogil=True, cache=True)
 def dense_measurement_smart(input_array, mes_amount, outcome_index):
 
@@ -455,7 +479,7 @@ def find_unique_markers(arr):
     return unique_marker
 
 
-@njit
+@njit(cache = True)
 def find_agreements(block_a, block_b):
 
     a_pointer = 0

src/qrisp/simulator/bi_arrays.py

@@ -27,7 +27,7 @@
 )
 
 import qrisp.simulator.bi_array_helper as hlp
-from qrisp.simulator.numerics_config import float_tresh
+from qrisp.simulator.numerics_config import float_tresh, sparsification_rate, cutoff_ratio
 
 try:
     # sparse_dot_mkl seems to be only faster in situations, where the shape of the
@@ -1067,11 +1067,20 @@ def mp_wrapper():
             self.apply_swaps()
             other.apply_swaps()
 
-            res.data = np.matmul(self.data, other.data)
+            res.data = np.matmul(self.data, other.data).ravel()
 
             self.swapaxes(0, 1)
             self.reshape(original_shape_self)
             other.reshape(original_shape_other)
+
+            if np.random.random(1)[0] < sparsification_rate and res.size > 2**14:
+                temp = np.abs(res.data.ravel())
+                max_abs = np.max(temp)
+                filter_arr = temp > max_abs*cutoff_ratio
+                res.data = res.data * filter_arr
+                res.data = res.data.reshape(res_shape)
+                res.sparsity = np.sum(filter_arr)/res.size
+
 
         if res.size > multithreading_threshold:
             # Start the wrapper

src/qrisp/simulator/circuit_preprocessing.py

@@ -763,7 +763,6 @@ def insert_multiverse_measurements(qc):
                         break
             else:
                 new_data.append(Instruction(disentangler, [meas_qubit]))
-                # new_measurements.append(instr)
                 new_measurements.append((instr.qubits[0], instr.clbits[0]))
                 continue
 
@@ -778,8 +777,6 @@ def insert_multiverse_measurements(qc):
 
             mes_instr = instr.copy()
             mes_instr.qubits = [qb]
-            # new_measurements.append((qb, instr.clbits[0]))
-            # new_measurements.append(mes_instr)
 
         elif instr.op.name == "reset":
 
@@ -857,4 +854,5 @@ def circuit_preprocessor(qc):
     else:
         qc = insert_disentangling(qc)
         qc = group_qc(qc)
+
     return reorder_circuit(qc, ["measure", "reset", "disentangle"])

src/qrisp/simulator/numerics_config.py

@@ -22,5 +22,6 @@
 
 # import cupy as xp
 
-float_tresh = 1e-7
+float_tresh = 1e-5
 cutoff_ratio = 5e-4
+sparsification_rate = 0.1