diff --git a/src/qrisp/core/quantum_array.py b/src/qrisp/core/quantum_array.py
index dd5f5956..1dd69f74 100644
--- a/src/qrisp/core/quantum_array.py
+++ b/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={},
diff --git a/src/qrisp/core/quantum_variable.py b/src/qrisp/core/quantum_variable.py
index cb352489..9dd99eee 100644
--- a/src/qrisp/core/quantum_variable.py
+++ b/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={},
diff --git a/src/qrisp/misc/utility.py b/src/qrisp/misc/utility.py
index 0c8ac890..31b4d2c4 100644
--- a/src/qrisp/misc/utility.py
+++ b/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
diff --git a/src/qrisp/qaoa/problems/maxCut.py b/src/qrisp/qaoa/problems/maxCut.py
index ad086238..cefcb6ac 100644
--- a/src/qrisp/qaoa/problems/maxCut.py
+++ b/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
 
diff --git a/src/qrisp/qaoa/qaoa_problem.py b/src/qrisp/qaoa/qaoa_problem.py
index 09496530..93299f19 100644
--- a/src/qrisp/qaoa/qaoa_problem.py
+++ b/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)
diff --git a/src/qrisp/qtypes/quantum_string.py b/src/qrisp/qtypes/quantum_string.py
index 79921ee8..d136c55a 100644
--- a/src/qrisp/qtypes/quantum_string.py
+++ b/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)
 
diff --git a/src/qrisp/simulator/bi_array_helper.py b/src/qrisp/simulator/bi_array_helper.py
index 1ecb03b1..d02225b4 100644
--- a/src/qrisp/simulator/bi_array_helper.py
+++ b/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
diff --git a/src/qrisp/simulator/bi_arrays.py b/src/qrisp/simulator/bi_arrays.py
index d29ceaae..17459a1a 100644
--- a/src/qrisp/simulator/bi_arrays.py
+++ b/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
diff --git a/src/qrisp/simulator/circuit_preprocessing.py b/src/qrisp/simulator/circuit_preprocessing.py
index 77d079e2..ce5bd21b 100644
--- a/src/qrisp/simulator/circuit_preprocessing.py
+++ b/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"])
diff --git a/src/qrisp/simulator/numerics_config.py b/src/qrisp/simulator/numerics_config.py
index 63789dc0..c5d1b09c 100644
--- a/src/qrisp/simulator/numerics_config.py
+++ b/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
diff --git a/src/qrisp/simulator/quantum_state.py b/src/qrisp/simulator/quantum_state.py
index c372bfa3..c43ae495 100644
--- a/src/qrisp/simulator/quantum_state.py
+++ b/src/qrisp/simulator/quantum_state.py
@@ -20,9 +20,12 @@
 
 from qrisp.simulator.numerics_config import xp
 from qrisp.simulator.tensor_factor import TensorFactor, multi_entangle
+from qrisp.simulator.bi_array_helper import permute_axes, invert_permutation
+import numpy as np
 
 from qrisp.misc.utility import bin_rep
 
+
 # This class describes a quantum state
 # This is achieved by another class, the TensorFactor.
 # Each qubit initially has its own tensor factor. When an entangling operation is
@@ -70,7 +73,7 @@ def apply_operation(self, operation, qubits):
         # Update the entangled factor on the involved qubits
         for i in range(len(involved_qubits)):
             self.tensor_factors[involved_qubits[i]] = entangled_factor
-
+            
         # Get the unitary of the operation
         unitary = operation.get_unitary()
 
@@ -106,8 +109,6 @@ def apply_operation(self, operation, qubits):
     def eval(self):
         [factor.unravel() for factor in self.tensor_factors]
 
-        # print(len(list(set(self.tensor_factors))))
-
         entangled_factor = multi_entangle(list(set(self.tensor_factors)))
 
         entangled_factor.unravel()
@@ -197,12 +198,14 @@ def multi_measure(self, mes_qubits, return_res_states = True):
             involved_factors.append(self.tensor_factors[mes_qubits[i]])
             involved_factors_dic[mes_qubits[i]] = self.tensor_factors[mes_qubits[i]]
             
-            
         involved_factors = list(set(involved_factors))
         
-        outcome_dic = {}
+        # involved_factors.sort(key = lambda x : len(x.qubits))
+        
+        grouped_qubits = []
         
-        outcome_dic_list = []
+        prob_array = np.ones(1)
+        ind_array = np.zeros(1, dtype = np.int64)
         
         for tf in involved_factors:
             
@@ -211,95 +214,31 @@ def multi_measure(self, mes_qubits, return_res_states = True):
                 if hash(involved_factors_dic[qb]) == hash(tf):
                     tf_mes_qubits.append(qb)
             
+            tf.unravel()
             p_list, tf_list, outcome_list = tf.multi_measure(tf_mes_qubits, return_res_tf = False)
             
+            prob_array = np.tensordot(prob_array, np.array(p_list), ((), ())).ravel()
+            grouped_qubits.extend(tf_mes_qubits)
             
-            dic_list = []
-            
-            for i in range(len(p_list)):
-                dic_list.append({"p" : p_list[i], 
-                                 "tf" : tf_list[i], 
-                                 "outcome" : bin_rep(outcome_list[i], len(tf_mes_qubits)),
-                                 "qubits" : tf_mes_qubits})
             
-            outcome_dic_list.append(dic_list)
-            outcome_dic[tf] = (p_list, tf_list, outcome_list)
-            
-        
-        outcome_states = []
-        p_list = []
-        outcome_list = []
-        
-        for dic_tuple in product(*outcome_dic_list):
-            
-            p = 1
-            if return_res_states:
-                outcome_state = self.copy()
-            outcome_value = [0]*len(mes_qubits)
-            
-            for i in range(len(dic_tuple)):
-                
-                p *= dic_tuple[i]["p"]
-                
-                for j in range(len(dic_tuple[i]["qubits"])):
-                    qb = dic_tuple[i]["qubits"][j]
-                    
-                    outcome_value[mes_qubits.index(qb)] = dic_tuple[i]["outcome"][j]
-                
-                    if return_res_states:
-                        outcome_state.tensor_factors[qb] = dic_tuple[i]["tf"]
-                    
-            
-            outcome = "".join(outcome_value)
-            outcome = int(outcome, 2)
-            outcome_list.append(outcome)
-            
-            if return_res_states:
-                outcome_states.append(outcome_state)
+            if len(grouped_qubits) < 63:
+                ind_array = tensorize_indices_jitted(np.array(outcome_list, dtype = np.int64), ind_array, len(tf_mes_qubits))
             else:
-                outcome_states.append(None)
+                if not isinstance(ind_array, list):
+                    ind_array = [int(i) for i in ind_array]
                 
-            p_list.append(p)
-            
-        
-        if not return_res_states:
-            return outcome_list, p_list
-            
-        
-        
-    # def multi_measure(self, mes_qubits, return_res_states = True):
-            
-        
-    #     tf = multi_entangle(list(set([self.tensor_factors[i] for i in mes_qubits])))
-        
-    #     p_list, tf_list, outcome_list = tf.multi_measure(mes_qubits, return_res_tf = return_res_states)
+                ind_array = tensorize_indices([int(i) for i in outcome_list], ind_array, len(tf_mes_qubits))
         
-    #     outcome_states = []
+        index_permutation = [mes_qubits.index(i) for i in grouped_qubits]
+        index_permutation = invert_permutation(np.array(index_permutation))[::-1]
         
-    #     for i in range(len(p_list)):
-            
-    #         if p_list[i] == 0:
-    #             outcome_states.append(None)
-    #             continue
-            
-    #         if return_res_states:
-    #             new_outcome_state = self.copy()
-                
-    #             for j in tf_list[i].qubits:
-    #                 new_outcome_state.tensor_factors[j] = tf_list[i]
-                    
-    #             for j in mes_qubits:
-    #                 new_outcome_state.tensor_factors[j] = TensorFactor([j])
-                    
-    #             outcome_states.append(new_outcome_state)
-                
-    #         else:
-                    
-    #             outcome_states.append(None)
+        if len(grouped_qubits) < 63:
+            ind_array = permute_axes(ind_array, index_permutation)    
+        else:
+            ind_array = permute_axes(ind_array, index_permutation, jit = False)
         
-    #     return p_list, outcome_states, outcome_list
+        return ind_array, prob_array
             
-
     def add_qubit(self):
         self.tensor_factors.append(TensorFactor([self.n]))
         self.n += 1
@@ -311,10 +250,37 @@ def disentangle(self, qb):
         disent_tf, remain_tf = tensor_factor.disentangle(qb)
         
         for i in tensor_factor.qubits:
-            if qb == i:
-                self.tensor_factors[i] = disent_tf
-            else:
-                self.tensor_factors[i] = remain_tf
-        
-        
-        
+            self.tensor_factors[i] = remain_tf
+            
+        self.tensor_factors[tensor_factor.qubits[0]] = disent_tf
+
+
+from numba import njit
+
+@njit(cache = True)
+def tensorize_indices_jitted(ind_array_0, ind_array_1, ind_0_size):
+    
+    n = len(ind_array_0)
+    m = len(ind_array_1)
+    
+    res = np.zeros((m,n), dtype = np.int64)
+    
+    for i in range(m):
+        for j in range(n):
+            res[i,j] = (ind_array_1[i]<<ind_0_size) | ind_array_0[j]
+    
+    return res.ravel()
+
+
+def tensorize_indices(ind_array_0, ind_array_1, ind_0_size):
+    
+    n = len(ind_array_0)
+    m = len(ind_array_1)
+    
+    res = []
+    
+    for i in range(m):
+        for j in range(n):
+            res.append((ind_array_1[i]<<ind_0_size) | ind_array_0[j])
+            
+    return res
\ No newline at end of file
diff --git a/src/qrisp/simulator/simulator.py b/src/qrisp/simulator/simulator.py
index 36f65b15..0a0a48ed 100644
--- a/src/qrisp/simulator/simulator.py
+++ b/src/qrisp/simulator/simulator.py
@@ -108,7 +108,6 @@ def run(qc, shots, token="", iqs=None, insert_reset=True):
         # if len(qc.qubits) < 30 or True:
             # qc, mes_list = extract_measurements(qc)
         
-            
         qc, new_mes_list = insert_multiverse_measurements(qc)
         
         mes_list = mes_list + new_mes_list
@@ -127,7 +126,6 @@ def run(qc, shots, token="", iqs=None, insert_reset=True):
             total_flops += 2 ** qc.data[i].op.num_qubits
         
         progress_bar.total = total_flops
-        
         for i in range(len(qc.data)):
             # Set alias for the instruction of this operation
             instr = qc.data[i]
@@ -166,7 +164,7 @@ def run(qc, shots, token="", iqs=None, insert_reset=True):
             mes_qubit_indices.append(qc.qubits.index(instr.qubits[0]))
             
         if len(mes_qubit_indices):
-            outcome_list, cl_prob = iqs.multi_measure(mes_qubit_indices, return_res_states = False)
+            outcome_list, cl_prob = iqs.multi_measure(mes_qubit_indices[::-1], return_res_states = False)
             mes_qubit_indices = []
             mes_clbit_indices = []
             
@@ -179,53 +177,48 @@ def run(qc, shots, token="", iqs=None, insert_reset=True):
         # and the probablities in .cl_prob. In order to ensure qiskit compatibility, we
         # reverse the bitstrings
         prob_dict = {}
-
         
-        for j in range(len(outcome_list)):
-            
-            outcome_str = bin(outcome_list[j])[2:].zfill(len(mes_list))
-            
-            try:
-                prob_dict[outcome_str] += cl_prob[j]
-            except KeyError:
-                prob_dict[outcome_str] = cl_prob[j]
-
-        #Normalize probabilities to 1 (can be slightly different due to
-        #floating point errors)
-        norm = 0
-        for v in prob_dict.values():
-            norm += v
-            
-        for k,v in prob_dict.items():
-            prob_dict[k] = prob_dict[k]/norm
+        norm = np.sum(cl_prob)
+        cl_prob = cl_prob/norm
+        
             
-
-        #If shots >= 10000, no samples will be drawn and the distribution will
+        res = {}
+        #If shots >= 1000000, no samples will be drawn and the distribution will
         #be returned instead
-        if shots >= 10000:
-            res = {}
+        if shots >= 100000:
             
-            for k, v in prob_dict.items():
-                res[k] = int(np.round(v*abs(shots)))
+            for j in range(len(outcome_list)):
                 
+                outcome_str = bin(outcome_list[j])[2:].zfill(len(mes_list))
+                
+                shot_val = int(np.round(cl_prob[j]*abs(shots)))
+                
+                try:
+                    res[outcome_str] += shot_val
+                except KeyError:
+                    res[outcome_str] = shot_val
+
         #Generate samples
         else:
             
             from numpy.random import choice
-            p_array = np.array(list(prob_dict.values()))
+            # p_array = np.array(list(prob_dict.values()))
             
-            samples = choice(len(p_array), shots, p=p_array)
+            # samples = choice(len(p_array), shots, p=p_array)
             
-            keys = list(prob_dict.keys())
+            samples = choice(len(cl_prob), shots, p=cl_prob)
             
             res = {}
             
             for s in samples:
-                if keys[s] in res:
-                    res[keys[s]] += 1
-                else:
-                    res[keys[s]] = 1
                 
+                outcome_str = bin(outcome_list[s])[2:].zfill(len(mes_list))
+                
+                if outcome_str in res:
+                    res[outcome_str] += 1
+                else:
+                    res[outcome_str] = 1
+        
         return res
 
 
diff --git a/tests/test_QAE.py b/tests/test_QAE.py
index fd4473fa..1f930cd5 100644
--- a/tests/test_QAE.py
+++ b/tests/test_QAE.py
@@ -55,6 +55,6 @@ def oracle_function(inp, tar):
 
     prec = 3 # precision
     res = QAE(input_list, state_function, oracle_function, precision=prec)
-    assert res.get_measurement() == {0.125: 0.3133, 0.875: 0.3133, 0.25: 0.1256, 0.75: 0.1256, 0.0: 0.051, 0.375: 0.0263, 0.625: 0.0263, 0.5: 0.0186}
+    assert res.get_measurement() == {0.125: 0.31334, 0.875: 0.31334, 0.25: 0.12557, 0.75: 0.12557, 0.0: 0.05096, 0.375: 0.02632, 0.625: 0.02632, 0.5: 0.01858}
 
 
diff --git a/tests/test_QAOAmaxClique.py b/tests/test_QAOAmaxClique.py
index 2d0bdde2..2c8b7689 100644
--- a/tests/test_QAOAmaxClique.py
+++ b/tests/test_QAOAmaxClique.py
@@ -102,7 +102,7 @@ def aClcostFct(state, G):
         secondQarg = QuantumVariable(secondGraph.number_of_nodes())
         secondInstance = QAOAProblem(cost_operator= maxCliqueCostOp(secondGraph), mixer= RX_mixer, cl_cost_function=maxCliqueCostfct(secondGraph)) 
         secondInstance.set_init_function(init_function=init_state)
-        theNiceQAOA2 = secondInstance.run(qarg=secondQarg, depth= 5)
+        theNiceQAOA2 = secondInstance.run(qarg=secondQarg, depth= 3, max_iter = 50, mes_kwargs = {"shots" : 100000})
         maxOne = sorted(theNiceQAOA2, key=theNiceQAOA2.get, reverse=True)[:1]
         #if no nodes connected, maxClique sol should just be the first node, i.e. "1000" (more or less zeros)
         if integ2 == 0:
@@ -116,7 +116,5 @@ def aClcostFct(state, G):
         
         if  testStr == maxOne:
             truth_list.append(1)
-    
+            
     assert sum(truth_list)/10 > 0.5
-
- 
diff --git a/tests/test_cycling_function.py b/tests/test_cycling_function.py
index 67f508e1..6ea873f6 100644
--- a/tests/test_cycling_function.py
+++ b/tests/test_cycling_function.py
@@ -50,5 +50,4 @@ def test_cycling_function():
 
     shift_amount[:] = {0: 1, -4 : 1, 1: 1}
     cyclic_shift(qa, shift_amount)
-    
-    assert qa.get_measurement() == {OutcomeArray([0, 1, 2, 3, 4, 5, 6, 7]): 0.3333, OutcomeArray([7, 0, 1, 2, 3, 4, 5, 6]): 0.3333, OutcomeArray([4, 5, 6, 7, 0, 1, 2, 3]): 0.3333}
+    assert qa.get_measurement() == {OutcomeArray([7, 0, 1, 2, 3, 4, 5, 6]): 0.33333, OutcomeArray([4, 5, 6, 7, 0, 1, 2, 3]): 0.33333, OutcomeArray([0, 1, 2, 3, 4, 5, 6, 7]): 0.33333}
diff --git a/tests/test_grovers_algorithm.py b/tests/test_grovers_algorithm.py
index b2e62a37..49d9d68e 100644
--- a/tests/test_grovers_algorithm.py
+++ b/tests/test_grovers_algorithm.py
@@ -110,8 +110,7 @@ def equation_oracle(qf_list):
         (isinstance(item, float) and item < 1) for item in mes_res.values()
     )
 
-    assert list(mes_res.keys())[0] == (0.5, -0.5)
-    assert list(mes_res.keys())[1] == (-0.5, 0.5)
+    assert set(list(mes_res.keys())[:2]) == {(0.5, -0.5), (-0.5, 0.5)}
 
     # check if the first tuple in the tuple contains expected values. Style example of the whole tuple: ((0, 1), 0.9592)
     assert all(
diff --git a/tests/test_qiskit_backend_client.py b/tests/test_qiskit_backend_client.py
index 481ea786..08413810 100644
--- a/tests/test_qiskit_backend_client.py
+++ b/tests/test_qiskit_backend_client.py
@@ -19,7 +19,7 @@
 # Created by ann81984 at 23.05.2022
 import numpy as np
 
-from qrisp.core import QuantumSession
+from qrisp import QuantumCircuit
 from qrisp.interface.backends import VirtualBackend
 from qrisp.interface.backends import VirtualQiskitBackend
 
@@ -29,46 +29,41 @@ def test_qiskit_backend_client():
     # TO-DO prevent this test from crashing regardless of functionality
     from qiskit import Aer
     # Create QuantumSession
-    qs = QuantumSession()
+    qc = QuantumCircuit()
 
-    qs.add_qubit()
-    qs.add_qubit()
-    qs.add_clbit()
+    qc.add_qubit()
+    qc.add_qubit()
+    qc.add_clbit()
 
-    qs.h(0)
+    qc.h(0)
 
-    qs.rz(np.pi / 2, 0)
+    qc.rz(np.pi / 2, 0)
 
-    qs.x(0)
-    qs.cx(0, 1)
-    qs.measure(1, 0)
-
-    qs.append(qs.to_op("composed_op"), qs.qubits, qs.clbits)
+    qc.x(0)
+    qc.cx(0, 1)
+    qc.measure(1, 0)
 
-    qs.append(qs.to_op("multi_composed_op"), qs.qubits, qs.clbits)
 
-    print(qs)
+    print(qc)
 
     ###################
 
     # Create VirtualBackend
-    def sample_run_func(qc, shots):
+    def sample_run_func(qc, shots, token = ""):
         print("Executing Circuit")
         return {"0": shots}
 
     test_virtual_backend = VirtualBackend(sample_run_func)
 
-    print(test_virtual_backend.run(qs, 100))
-    assert str(test_virtual_backend.run(qs, 100)) == "{'0': 100}"
-    assert test_virtual_backend.run(qs, 100)["0"] == 100
+    print(test_virtual_backend.run(qc, 100))
+    assert str(test_virtual_backend.run(qc, 100)) == "{'0': 100}"
+    assert test_virtual_backend.run(qc, 100)["0"] == 100
 
     ###################
 
     # Create Qiskit Backend
     test_qiskit_backend = VirtualQiskitBackend(Aer.get_backend("qasm_simulator"))
 
-    from qrisp import QuantumCircuit
-
     qc = QuantumCircuit(4, 1)
     qc.x(0)
     qc.measure(0, 0)
diff --git a/tests/test_string_test.py b/tests/test_string_test.py
index 38bdc83f..7b2ec3fe 100644
--- a/tests/test_string_test.py
+++ b/tests/test_string_test.py
@@ -41,6 +41,7 @@ def test_string_test():
 
     q_str_3 = q_str.duplicate(init=True)
 
+    
     assert (type(q_str).__name__ and type(q_str_3).__name__) == "QuantumString"
     assert str(q_str + q_str_3) == "{'hello world! hello world! ': 1.0}"
     print(q_str + q_str_3)
diff --git a/tests/test_uncomputation_example.py b/tests/test_uncomputation_example.py
index b2403d52..edc88c6c 100644
--- a/tests/test_uncomputation_example.py
+++ b/tests/test_uncomputation_example.py
@@ -71,17 +71,16 @@ def sqrt_oracle(qf):
         sqrt_oracle(qf)
         diffuser(qf)
         print(qf)
-
-    assert qf.get_measurement() == {
-        0.5: 0.9453,
-        0.0: 0.0078,
-        1.0: 0.0078,
-        1.5: 0.0078,
-        2.0: 0.0078,
-        2.5: 0.0078,
-        3.0: 0.0078,
-        3.5: 0.0078,
-    }
+    
+    print(qf)
+    assert qf.get_measurement() == {0.5: 0.94531, 
+                                    0.0: 0.00781, 
+                                    1.0: 0.00781, 
+                                    1.5: 0.00781, 
+                                    2.0: 0.00781, 
+                                    2.5: 0.00781, 
+                                    3.0: 0.00781, 
+                                    3.5: 0.00781}
 
     # ---------
     print("Test 2 passed")