diff --git a/aquapointer/analog/utils/__init__.py b/aquapointer/analog/utils/__init__.py
index 26ea7a2..b046d36 100644
--- a/aquapointer/analog/utils/__init__.py
+++ b/aquapointer/analog/utils/__init__.py
@@ -4,4 +4,8 @@
 # LICENSE file in the root directory of this source tree.
 
 
-from aquapointer.analog.utils import benchmark_utils, density_utils, detuning_scale_utils
+from aquapointer.analog.utils import (
+    benchmark_utils,
+    density_utils,
+    detuning_scale_utils,
+)
diff --git a/aquapointer/analog/utils/benchmark_utils.py b/aquapointer/analog/utils/benchmark_utils.py
index 6fcbe72..89979fb 100644
--- a/aquapointer/analog/utils/benchmark_utils.py
+++ b/aquapointer/analog/utils/benchmark_utils.py
@@ -14,31 +14,37 @@
 import qutip
 import collections
 
+
 def generate_binary_strings(bit_count):
     """Generates all binary strings of length `bit_count`"""
     binary_strings = []
-    def genbin(n, bs=''):
+
+    def genbin(n, bs=""):
         if len(bs) == n:
             binary_strings.append(bs)
         else:
-            genbin(n, bs + '0')
-            genbin(n, bs + '1')
+            genbin(n, bs + "0")
+            genbin(n, bs + "1")
+
     genbin(bit_count)
     return binary_strings
 
+
 def probability_from_0(bitstring, epsilon):
-    """Calculates the probability of obtaining `bitstring` from '0000..0' 
-    assuming an `epsilon` probability of measuring '1' instead of '0' """
-    n1 = bitstring.count('1')
-    n0 = bitstring.count('0')
-    return (epsilon**n1) * ((1-epsilon)**n0)
+    """Calculates the probability of obtaining `bitstring` from '0000..0'
+    assuming an `epsilon` probability of measuring '1' instead of '0'"""
+    n1 = bitstring.count("1")
+    n0 = bitstring.count("0")
+    return (epsilon**n1) * ((1 - epsilon) ** n0)
+
 
 def bitstring_projector(bitstring):
     """Generates the projector corresponding to the computational basis state
     `bitstring`"""
-    inv_bit = bitstring.replace('1', '2').replace('0', '1').replace('2', '0')
+    inv_bit = bitstring.replace("1", "2").replace("0", "1").replace("2", "0")
     return qutip.ket(inv_bit).proj()
 
+
 def flatten_samples(samples):
     """Takes a counter `samples` and returns a list where each
     entry appears as many times as its count in `samples`"""
@@ -47,6 +53,7 @@ def flatten_samples(samples):
         flat_samples.extend([key for _ in range(val)])
     return flat_samples
 
+
 def magnetization_observable(n, i):
     """Generates the magnetization operator of qubit `i`
     out of `n` qubits"""
@@ -54,6 +61,7 @@ def magnetization_observable(n, i):
     obs[i] = (qutip.sigmaz() + qutip.qeye(2)) / 2
     return qutip.tensor(obs)
 
+
 def bitstring_probability_exact(bitstring, results, t_i, t_f):
     """Calculates the exact probability of measuring `bitstring` between
     times `t_i` and `t_f` for a Pulser emulation that gives `results`
@@ -62,11 +70,13 @@ def bitstring_probability_exact(bitstring, results, t_i, t_f):
     obs = bitstring_projector(bitstring)
     return results.expect([obs])[0][t_i:t_f]
 
+
 def bitstring_probability_sampling(bitstring, samples):
     """Calculates the probability of measuring `bitstring` from a
     sampling `samples` of a state"""
     tot = np.sum(list(samples.values()))
-    return samples[bitstring]/tot 
+    return samples[bitstring] / tot
+
 
 def bitstring_probability_sampling_binned(bitstring, samples, nbins):
     """Calculates the probability of measuring `bitstring` from a
@@ -79,48 +89,52 @@ def bitstring_probability_sampling_binned(bitstring, samples, nbins):
     for i in range(nbins):
         samples_i = collections.Counter(binned_samples[i])
         tot = np.sum(list(samples_i.values()))
-        probs[i] = samples_i[bitstring]/tot 
+        probs[i] = samples_i[bitstring] / tot
     return (np.mean(probs), np.std(probs))
 
+
 def average_bitstring_probability_sampling(bitstring, results, n_samples, repetitions):
     """Calculates the average probability of measuring `bitstring` out of
     `repetitions` inependent samplings, each containing `n_samples` samples, from a Pulser
     emulation that gives `results`
     (an object pulser_simulation.simresults.CoherentResults/NoisyResults) as the output
     of an emulation"""
-    probs = [] 
+    probs = []
     for i in range(repetitions):
         samples = results.sample_final_state(n_samples)
         probs.append(bitstring_probability_sampling(bitstring, samples))
-    
+
     return (np.mean(probs), np.std(probs))
 
+
 def bitstring_distribution_exact(results):
-    """ Calculate the exact probability of every bitstring, in order
-    from '1111..1' to '000..0' """
+    """Calculate the exact probability of every bitstring, in order
+    from '1111..1' to '000..0'"""
     state = results.get_final_state()
     N = state.shape[0]
     probs = np.zeros(N)
     for i, c in enumerate(state):
-        probs[i] = np.real(c)**2 + np.imag(c)**2
+        probs[i] = np.real(c) ** 2 + np.imag(c) ** 2
     return probs
 
+
 def bitstring_distribution_sampling(samples, n):
-    """ Calculate the probability of every bitstring of length `n` from
-    sampling `samples`, in order from '1111..1' to '000..0' """
+    """Calculate the probability of every bitstring of length `n` from
+    sampling `samples`, in order from '1111..1' to '000..0'"""
     bitstrings = np.flip(generate_binary_strings(n))
     probs = np.zeros(2**n)
     tot = np.sum(list(samples.values()))
     for i, bitstring in enumerate(bitstrings):
         try:
-            probs[i] = samples[bitstring]/tot
+            probs[i] = samples[bitstring] / tot
         except KeyError:
             probs[i] = 0
     return probs
 
+
 def bitstring_distribution_sampling_binned(samples, n, nbins):
-    """ Calculate the probability of every bitstring of length `n` from
-    sampling `samples`, in order from '1111..1' to '000..0' 
+    """Calculate the probability of every bitstring of length `n` from
+    sampling `samples`, in order from '1111..1' to '000..0'
     Split the samples in `nbins` bins to get error estimate"""
     flat_samples = flatten_samples(samples)
     np.random.shuffle(flat_samples)
@@ -130,9 +144,10 @@ def bitstring_distribution_sampling_binned(samples, n, nbins):
         samples_i = collections.Counter(binned_samples[i])
         dist = bitstring_distribution_sampling(samples_i, n)
         for j, d in enumerate(dist):
-            probs[i,j] = d
+            probs[i, j] = d
     return (np.mean(probs, axis=0), np.std(probs, axis=0))
 
+
 def average_bitstring_distribution_sampling(results, n_samples, repetitions=1000):
     n = results._size
     probs = np.zeros((repetitions, 2**n))
@@ -140,9 +155,10 @@ def average_bitstring_distribution_sampling(results, n_samples, repetitions=1000
         samples = results.sample_final_state(n_samples)
         dist = bitstring_distribution_sampling(samples, n)
         for j, d in enumerate(dist):
-            probs[i,j] = d
+            probs[i, j] = d
     return (np.mean(probs, axis=0), np.std(probs, axis=0))
-    
+
+
 def magnetization_exact(k, n, results, t_i, t_f):
     """Calculate the exact magnetization of qubit `k` out
     of `n` qubits between times `t_i` and `t_f` for a Pulser emulation that
@@ -150,16 +166,18 @@ def magnetization_exact(k, n, results, t_i, t_f):
     as the output of a noiseless emulation"""
     return results.expect([magnetization_observable(n, k)])[0][t_i:t_f]
 
+
 def magnetization_sampling(k, samples):
     """Calculate the magnetization of qubit `k` from a sampling
     `samples` of a state"""
     mag = 0
     tot = 0
     for key, val in samples.items():
-        if key[k] == '1':
+        if key[k] == "1":
             mag += val
         tot += val
-    return mag/tot
+    return mag / tot
+
 
 def magnetization_sampling_binned(k, samples, nbins):
     """Calculate the magnetization of qubit `k` from a sampling
@@ -174,12 +192,13 @@ def magnetization_sampling_binned(k, samples, nbins):
         mag = 0
         tot = 0
         for key, val in samples_i.items():
-            if key[k] == '1':
+            if key[k] == "1":
                 mag += val
             tot += val
-        mags[i] = mag/tot
+        mags[i] = mag / tot
     return (np.mean(mags), np.std(mags))
 
+
 def average_magnetization_sampling(k, results, n_samples, repetitions):
     """Calculates the average magnetization of qubit `k` out of
     `repetitions` inependent samplings, each contining `n_samples` samples, from a Pulser
@@ -199,16 +218,17 @@ def is_IS(bitstring, pos, brad):
 
     for i in range(len(bitstring)):
         b1 = bitstring[i]
-        if b1=='0':
+        if b1 == "0":
             continue
-        for j in range(i+1, len(bitstring)):
+        for j in range(i + 1, len(bitstring)):
             b2 = bitstring[j]
-            if b2=='0':
+            if b2 == "0":
                 continue
-            if np.linalg.norm(pos[i]-pos[j]) < brad:
+            if np.linalg.norm(pos[i] - pos[j]) < brad:
                 return False
     return True
 
+
 def separate_IS(bitstrings, pos, brad):
     """Among all bitstrings in the list `bitstrings`, returns only those that
     are independent sets of the graph defined by nodes in position `pos` and
@@ -219,27 +239,32 @@ def separate_IS(bitstrings, pos, brad):
             res.append(bitstring)
     return res
 
+
 def generate_random_register(rows, cols, n, spacing, seed=346):
     """Generates a random register with `n` qubits as a subregister of a
     triangular lattice or size `rows` x `cols` and spacing `spacing`"""
-    if n>rows*cols:
+    if n > rows * cols:
         raise ValueError("lattice is not large enough")
     np.random.seed(seed)
     reg = Register.triangular_lattice(rows=rows, atoms_per_row=cols, spacing=spacing)
     pos = reg._coords
-    while len(pos)>n:
+    while len(pos) > n:
         k = np.random.randint(0, len(pos))
-        pos = [pos[i] for i in range(len(pos)) if i!=k]
+        pos = [pos[i] for i in range(len(pos)) if i != k]
     reg = Register.from_coordinates(pos)
     return reg, pos
 
+
 def observable(bitstring):
     """Generates the projector corresponding to the computational basis state
     `bitstring`"""
-    inv_bit = bitstring.replace('1', '2').replace('0', '1').replace('2', '0')
+    inv_bit = bitstring.replace("1", "2").replace("0", "1").replace("2", "0")
     return qutip.ket(inv_bit).proj()
 
-def pulse_landscape(register, device, independent_sets, omega, time_list, delta_list, verbose=True):
+
+def pulse_landscape(
+    register, device, independent_sets, omega, time_list, delta_list, verbose=True
+):
     """For each state in `independent_sets`, returns expectation value (and
     various other retlated stats) of its projector over a time range defined by
     `time_list` and over a range of detuning values defined by `delta_list`"""
@@ -250,18 +275,21 @@ def pulse_landscape(register, device, independent_sets, omega, time_list, delta_
     IS_dictionary = dict()
     for bitstring in independent_sets:
         IS_dictionary[bitstring] = {
-            "excitations" : bitstring.count('1'),
-            "landscape" : np.zeros((len(time_list), len(delta_list))),
-            "max_prob" : 0,
-            "min_prob" : 0,
-            "location_max" : (-1,-1),
-            "location_min" : (-1,-1),
+            "excitations": bitstring.count("1"),
+            "landscape": np.zeros((len(time_list), len(delta_list))),
+            "max_prob": 0,
+            "min_prob": 0,
+            "location_max": (-1, -1),
+            "location_min": (-1, -1),
         }
 
     T = time_list[-1]
 
     for j, d in enumerate(delta_list):
-        print(f"Simulating detuning {j+1} of {len(delta_list)}" if verbose else "", end = "\n" if verbose else "")
+        print(
+            f"Simulating detuning {j+1} of {len(delta_list)}" if verbose else "",
+            end="\n" if verbose else "",
+        )
         seq = Sequence(register, device)
         seq.declare_channel("ch", "rydberg_global")
         pulse = Pulse.ConstantPulse(T, omega, d, 0)
@@ -272,27 +300,27 @@ def pulse_landscape(register, device, independent_sets, omega, time_list, delta_
             obs = [observable(bitstring)]
             exp = res.expect(obs)
             for i, t in enumerate(time_list):
-                IS_dictionary[bitstring]["landscape"][i,j] = exp[0][t]
+                IS_dictionary[bitstring]["landscape"][i, j] = exp[0][t]
 
     for bitstring in independent_sets:
-
-        IS_dictionary[bitstring]["max_prob"] = np.max(IS_dictionary[bitstring]["landscape"])
-        IS_dictionary[bitstring]["min_prob"] = np.min(IS_dictionary[bitstring]["landscape"])
-        IS_dictionary[bitstring]["location_max"] = (
-            np.unravel_index(
-                IS_dictionary[bitstring]["landscape"].argmax(),
-                IS_dictionary[bitstring]["landscape"].shape
-            )
+        IS_dictionary[bitstring]["max_prob"] = np.max(
+            IS_dictionary[bitstring]["landscape"]
         )
-        IS_dictionary[bitstring]["location_min"] = (
-            np.unravel_index(
-                IS_dictionary[bitstring]["landscape"].argmin(),
-                IS_dictionary[bitstring]["landscape"].shape
-            )
+        IS_dictionary[bitstring]["min_prob"] = np.min(
+            IS_dictionary[bitstring]["landscape"]
+        )
+        IS_dictionary[bitstring]["location_max"] = np.unravel_index(
+            IS_dictionary[bitstring]["landscape"].argmax(),
+            IS_dictionary[bitstring]["landscape"].shape,
+        )
+        IS_dictionary[bitstring]["location_min"] = np.unravel_index(
+            IS_dictionary[bitstring]["landscape"].argmin(),
+            IS_dictionary[bitstring]["landscape"].shape,
         )
 
     return IS_dictionary
 
+
 def apply_meas_errors(samples, eps, eps_prime):
     """Copied from Pulser"""
     shots = list(samples.keys())
@@ -305,21 +333,32 @@ def apply_meas_errors(samples, eps, eps_prime):
     # Repeat flip_probs based on n_detects_list
     flip_probs_repeated = np.repeat(flip_probs, n_detects_list, axis=0)
     # Generate random matrix of shape (sum(n_detects_list), len(shot))
-    random_matrix = np.random.uniform(
-        size=(np.sum(n_detects_list), len(shot_arr[0]))
-    )
+    random_matrix = np.random.uniform(size=(np.sum(n_detects_list), len(shot_arr[0])))
     # Compare random matrix with flip probabilities
     flips = random_matrix < flip_probs_repeated
     # Perform XOR between original array and flips
     new_shots = shot_arr.repeat(n_detects_list, axis=0) ^ flips
     # Count all the new_shots
     # We are not converting to str before because tuple indexing is faster
-    detected_sample_dict: collections.Counter = collections.Counter(map(tuple, new_shots))
+    detected_sample_dict: collections.Counter = collections.Counter(
+        map(tuple, new_shots)
+    )
     return collections.Counter(
         {"".join(map(str, k)): v for k, v in detected_sample_dict.items()}
     )
 
-def pulse_landscape_SPAM(register, device, independent_sets, omega, time_list, delta_list, epsilon, epsilon_prime, verbose=True):
+
+def pulse_landscape_SPAM(
+    register,
+    device,
+    independent_sets,
+    omega,
+    time_list,
+    delta_list,
+    epsilon,
+    epsilon_prime,
+    verbose=True,
+):
     """For each state in `independent_sets`, returns expectation value (and
     various other retlated stats) of its projector over a time range defined by
     `time_list` and over a range of detuning values defined by `delta_list`.
@@ -331,18 +370,21 @@ def pulse_landscape_SPAM(register, device, independent_sets, omega, time_list, d
     IS_dictionary = dict()
     for bitstring in independent_sets:
         IS_dictionary[bitstring] = {
-            "excitations" : bitstring.count('1'),
-            "landscape" : np.zeros((len(time_list), len(delta_list))),
-            "max_prob" : 0,
-            "min_prob" : 0,
-            "location_max" : (-1,-1),
-            "location_min" : (-1,-1),
+            "excitations": bitstring.count("1"),
+            "landscape": np.zeros((len(time_list), len(delta_list))),
+            "max_prob": 0,
+            "min_prob": 0,
+            "location_max": (-1, -1),
+            "location_min": (-1, -1),
         }
 
     T = time_list[-1]
 
     for j, d in enumerate(delta_list):
-        print(f"Simulating detuning {j+1} of {len(delta_list)}" if verbose else "", end = "\n" if verbose else "")
+        print(
+            f"Simulating detuning {j+1} of {len(delta_list)}" if verbose else "",
+            end="\n" if verbose else "",
+        )
         seq = Sequence(register, device)
         seq.declare_channel("ch", "rydberg_global")
         pulse = Pulse.ConstantPulse(T, omega, d, 0)
@@ -351,29 +393,39 @@ def pulse_landscape_SPAM(register, device, independent_sets, omega, time_list, d
         res = sim.run(progress_bar=verbose)
         for bitstring in independent_sets:
             for i, t in enumerate(time_list):
-                samples = res.sample_state(t/1000, n_samples=10000)
+                samples = res.sample_state(t / 1000, n_samples=10000)
                 err_samples = apply_meas_errors(samples, epsilon, epsilon_prime)
-                IS_dictionary[bitstring]["landscape"][i,j] = err_samples[bitstring]/10000
+                IS_dictionary[bitstring]["landscape"][i, j] = (
+                    err_samples[bitstring] / 10000
+                )
     for bitstring in independent_sets:
-
-        IS_dictionary[bitstring]["max_prob"] = np.max(IS_dictionary[bitstring]["landscape"])
-        IS_dictionary[bitstring]["min_prob"] = np.min(IS_dictionary[bitstring]["landscape"])
-        IS_dictionary[bitstring]["location_max"] = (
-            np.unravel_index(
-                IS_dictionary[bitstring]["landscape"].argmax(),
-                IS_dictionary[bitstring]["landscape"].shape
-            )
+        IS_dictionary[bitstring]["max_prob"] = np.max(
+            IS_dictionary[bitstring]["landscape"]
         )
-        IS_dictionary[bitstring]["location_min"] = (
-            np.unravel_index(
-                IS_dictionary[bitstring]["landscape"].argmin(),
-                IS_dictionary[bitstring]["landscape"].shape
-            )
+        IS_dictionary[bitstring]["min_prob"] = np.min(
+            IS_dictionary[bitstring]["landscape"]
+        )
+        IS_dictionary[bitstring]["location_max"] = np.unravel_index(
+            IS_dictionary[bitstring]["landscape"].argmax(),
+            IS_dictionary[bitstring]["landscape"].shape,
+        )
+        IS_dictionary[bitstring]["location_min"] = np.unravel_index(
+            IS_dictionary[bitstring]["landscape"].argmin(),
+            IS_dictionary[bitstring]["landscape"].shape,
         )
     return IS_dictionary
 
 
-def pulse_landscape_samples(register, device, independent_sets, omega, time_list, delta_list, verbose=True, n_samples=150):
+def pulse_landscape_samples(
+    register,
+    device,
+    independent_sets,
+    omega,
+    time_list,
+    delta_list,
+    verbose=True,
+    n_samples=150,
+):
     """For each state in `independent_sets`, returns expectation value (and
     various other retlated stats) of its projector over a time range defined by
     `time_list` and over a range of detuning values defined by `delta_list`"""
@@ -384,19 +436,22 @@ def pulse_landscape_samples(register, device, independent_sets, omega, time_list
     IS_dictionary = dict()
     for bitstring in independent_sets:
         IS_dictionary[bitstring] = {
-            "excitations" : bitstring.count('1'),
-            "landscape" : np.zeros((len(time_list), len(delta_list))),
-            "samples" : np.empty((len(time_list), len(delta_list)), dtype=object),
-            "max_prob" : 0,
-            "min_prob" : 0,
-            "location_max" : (-1,-1),
-            "location_min" : (-1,-1),
+            "excitations": bitstring.count("1"),
+            "landscape": np.zeros((len(time_list), len(delta_list))),
+            "samples": np.empty((len(time_list), len(delta_list)), dtype=object),
+            "max_prob": 0,
+            "min_prob": 0,
+            "location_max": (-1, -1),
+            "location_min": (-1, -1),
         }
 
     T = time_list[-1]
 
     for j, d in enumerate(delta_list):
-        print(f"Simulating detuning {j+1} of {len(delta_list)}" if verbose else "", end = "\n" if verbose else "")
+        print(
+            f"Simulating detuning {j+1} of {len(delta_list)}" if verbose else "",
+            end="\n" if verbose else "",
+        )
         seq = Sequence(register, device)
         seq.declare_channel("ch", "rydberg_global")
         pulse = Pulse.ConstantPulse(T, omega, d, 0)
@@ -405,35 +460,43 @@ def pulse_landscape_samples(register, device, independent_sets, omega, time_list
         res = sim.run(progress_bar=verbose)
         for bitstring in independent_sets:
             for i, t in enumerate(time_list):
-                samples = res.sample_state(t/1e3, n_samples=n_samples)
+                samples = res.sample_state(t / 1e3, n_samples=n_samples)
                 val = 0
                 for key, count in samples.items():
-                    if key==bitstring:
-                        val = count/n_samples
-                IS_dictionary[bitstring]["landscape"][i,j] = val
-                IS_dictionary[bitstring]["samples"][i,j] = samples
+                    if key == bitstring:
+                        val = count / n_samples
+                IS_dictionary[bitstring]["landscape"][i, j] = val
+                IS_dictionary[bitstring]["samples"][i, j] = samples
 
     for bitstring in independent_sets:
-
-        IS_dictionary[bitstring]["max_prob"] = np.max(IS_dictionary[bitstring]["landscape"])
-        IS_dictionary[bitstring]["min_prob"] = np.min(IS_dictionary[bitstring]["landscape"])
-        IS_dictionary[bitstring]["location_max"] = (
-            np.unravel_index(
-                IS_dictionary[bitstring]["landscape"].argmax(),
-                IS_dictionary[bitstring]["landscape"].shape
-            )
+        IS_dictionary[bitstring]["max_prob"] = np.max(
+            IS_dictionary[bitstring]["landscape"]
         )
-        IS_dictionary[bitstring]["location_min"] = (
-            np.unravel_index(
-                IS_dictionary[bitstring]["landscape"].argmin(),
-                IS_dictionary[bitstring]["landscape"].shape
-            )
+        IS_dictionary[bitstring]["min_prob"] = np.min(
+            IS_dictionary[bitstring]["landscape"]
+        )
+        IS_dictionary[bitstring]["location_max"] = np.unravel_index(
+            IS_dictionary[bitstring]["landscape"].argmax(),
+            IS_dictionary[bitstring]["landscape"].shape,
+        )
+        IS_dictionary[bitstring]["location_min"] = np.unravel_index(
+            IS_dictionary[bitstring]["landscape"].argmin(),
+            IS_dictionary[bitstring]["landscape"].shape,
         )
 
     return IS_dictionary
 
 
-def pulse_landscape_selected(register, device, independent_sets, omega, shape, params_list, indices_list, verbose=True):
+def pulse_landscape_selected(
+    register,
+    device,
+    independent_sets,
+    omega,
+    shape,
+    params_list,
+    indices_list,
+    verbose=True,
+):
     """For each state in `independent_sets`, returns expectation value (and
     various other retlated stats) of its projector over a range of parameters
     defined by `params_list`"""
@@ -441,17 +504,20 @@ def pulse_landscape_selected(register, device, independent_sets, omega, shape, p
     IS_dictionary = dict()
     for bitstring in independent_sets:
         IS_dictionary[bitstring] = {
-            "excitations" : bitstring.count('1'),
-            "landscape" : np.zeros(shape),
-            "max_prob" : 0,
-            "min_prob" : 0,
-            "location_max" : (-1,-1),
-            "location_min" : (-1,-1),
+            "excitations": bitstring.count("1"),
+            "landscape": np.zeros(shape),
+            "max_prob": 0,
+            "min_prob": 0,
+            "location_max": (-1, -1),
+            "location_min": (-1, -1),
         }
 
     for k, (t, d) in enumerate(params_list):
-        print(f"Simulating parameters {k+1} of {len(params_list)}" if verbose else "", end = "\n" if verbose else "")
-        print(f"time: {t}     det: {d}" if verbose else "", end = "\n" if verbose else "")
+        print(
+            f"Simulating parameters {k+1} of {len(params_list)}" if verbose else "",
+            end="\n" if verbose else "",
+        )
+        print(f"time: {t}     det: {d}" if verbose else "", end="\n" if verbose else "")
         seq = Sequence(register, device)
         seq.declare_channel("ch", "rydberg_global")
         pulse = Pulse.ConstantPulse(t, omega, d, 0)
@@ -463,30 +529,32 @@ def pulse_landscape_selected(register, device, independent_sets, omega, shape, p
             exp = res.expect(obs)
             i = indices_list[k][0]
             j = indices_list[k][1]
-            print(f"i: {i}     j: {j}" if verbose else "", end = "\n" if verbose else "")
-            print(f"p: {exp[0][-1]}" if verbose else "", end = "\n" if verbose else "")
-            IS_dictionary[bitstring]["landscape"][i,j] = exp[0][-1]
+            print(f"i: {i}     j: {j}" if verbose else "", end="\n" if verbose else "")
+            print(f"p: {exp[0][-1]}" if verbose else "", end="\n" if verbose else "")
+            IS_dictionary[bitstring]["landscape"][i, j] = exp[0][-1]
 
     for bitstring in independent_sets:
-
-        IS_dictionary[bitstring]["max_prob"] = np.max(IS_dictionary[bitstring]["landscape"])
-        IS_dictionary[bitstring]["min_prob"] = np.min(IS_dictionary[bitstring]["landscape"])
-        IS_dictionary[bitstring]["location_max"] = (
-            np.unravel_index(
-                IS_dictionary[bitstring]["landscape"].argmax(),
-                IS_dictionary[bitstring]["landscape"].shape
-            )
+        IS_dictionary[bitstring]["max_prob"] = np.max(
+            IS_dictionary[bitstring]["landscape"]
         )
-        IS_dictionary[bitstring]["location_min"] = (
-            np.unravel_index(
-                IS_dictionary[bitstring]["landscape"].argmin(),
-                IS_dictionary[bitstring]["landscape"].shape
-            )
+        IS_dictionary[bitstring]["min_prob"] = np.min(
+            IS_dictionary[bitstring]["landscape"]
+        )
+        IS_dictionary[bitstring]["location_max"] = np.unravel_index(
+            IS_dictionary[bitstring]["landscape"].argmax(),
+            IS_dictionary[bitstring]["landscape"].shape,
+        )
+        IS_dictionary[bitstring]["location_min"] = np.unravel_index(
+            IS_dictionary[bitstring]["landscape"].argmin(),
+            IS_dictionary[bitstring]["landscape"].shape,
         )
 
     return IS_dictionary
 
-def pulse_landscape_modulation(register, device, independent_sets, omega, time_list, delta_list, verbose=True):
+
+def pulse_landscape_modulation(
+    register, device, independent_sets, omega, time_list, delta_list, verbose=True
+):
     """For each state in `independent_sets`, returns expectation value (and
     various other retlated stats) of its projector over a time range defined by
     `time_list` and over a range of detuning values defined by `delta_list`.
@@ -499,18 +567,26 @@ def pulse_landscape_modulation(register, device, independent_sets, omega, time_l
     IS_dictionary = dict()
     for bitstring in independent_sets:
         IS_dictionary[bitstring] = {
-            "excitations" : bitstring.count('1'),
-            "landscape" : np.zeros((len(time_list), len(delta_list))),
-            "max_prob" : 0,
-            "min_prob" : 0,
-            "location_max" : (-1,-1),
-            "location_min" : (-1,-1),
+            "excitations": bitstring.count("1"),
+            "landscape": np.zeros((len(time_list), len(delta_list))),
+            "max_prob": 0,
+            "min_prob": 0,
+            "location_max": (-1, -1),
+            "location_min": (-1, -1),
         }
 
     for i, t in enumerate(time_list):
-        print(f"Simulating time {i+1} of {len(time_list)}" if verbose else "", end = "\n" if verbose else "")
+        print(
+            f"Simulating time {i+1} of {len(time_list)}" if verbose else "",
+            end="\n" if verbose else "",
+        )
         for j, d in enumerate(delta_list):
-            print(f"    Simulating detuning {j+1} of {len(delta_list)}" if verbose else "", end = "\n" if verbose else "")
+            print(
+                f"    Simulating detuning {j+1} of {len(delta_list)}"
+                if verbose
+                else "",
+                end="\n" if verbose else "",
+            )
             seq = Sequence(register, device)
             seq.declare_channel("ch", "rydberg_global")
             pulse = Pulse.ConstantPulse(t, omega, d, 0)
@@ -520,28 +596,38 @@ def pulse_landscape_modulation(register, device, independent_sets, omega, time_l
             for bitstring in independent_sets:
                 obs = [observable(bitstring)]
                 exp = res.expect(obs)
-                IS_dictionary[bitstring]["landscape"][i,j] = exp[0][-1]
+                IS_dictionary[bitstring]["landscape"][i, j] = exp[0][-1]
 
     for bitstring in independent_sets:
-
-        IS_dictionary[bitstring]["max_prob"] = np.max(IS_dictionary[bitstring]["landscape"])
-        IS_dictionary[bitstring]["min_prob"] = np.min(IS_dictionary[bitstring]["landscape"])
-        IS_dictionary[bitstring]["location_max"] = (
-            np.unravel_index(
-                IS_dictionary[bitstring]["landscape"].argmax(),
-                IS_dictionary[bitstring]["landscape"].shape
-            )
+        IS_dictionary[bitstring]["max_prob"] = np.max(
+            IS_dictionary[bitstring]["landscape"]
         )
-        IS_dictionary[bitstring]["location_min"] = (
-            np.unravel_index(
-                IS_dictionary[bitstring]["landscape"].argmin(),
-                IS_dictionary[bitstring]["landscape"].shape
-            )
+        IS_dictionary[bitstring]["min_prob"] = np.min(
+            IS_dictionary[bitstring]["landscape"]
+        )
+        IS_dictionary[bitstring]["location_max"] = np.unravel_index(
+            IS_dictionary[bitstring]["landscape"].argmax(),
+            IS_dictionary[bitstring]["landscape"].shape,
+        )
+        IS_dictionary[bitstring]["location_min"] = np.unravel_index(
+            IS_dictionary[bitstring]["landscape"].argmin(),
+            IS_dictionary[bitstring]["landscape"].shape,
         )
 
     return IS_dictionary
 
-def pulse_landscape_modulation_fast(register, device, time_buffer, independent_sets, omega, time_list, delta_list, EOM=True, verbose=True):
+
+def pulse_landscape_modulation_fast(
+    register,
+    device,
+    time_buffer,
+    independent_sets,
+    omega,
+    time_list,
+    delta_list,
+    EOM=True,
+    verbose=True,
+):
     """For each state in `independent_sets`, returns expectation value (and
     various other retlated stats) of its projector over a time range defined by
     `time_list` and over a range of detuning values defined by `delta_list`.
@@ -555,19 +641,24 @@ def pulse_landscape_modulation_fast(register, device, time_buffer, independent_s
     IS_dictionary = dict()
     for bitstring in independent_sets:
         IS_dictionary[bitstring] = {
-            "excitations" : bitstring.count('1'),
-            "landscape" : np.zeros((len(time_list), len(delta_list))),
-            "max_prob" : 0,
-            "min_prob" : 0,
-            "location_max" : (-1,-1),
-            "location_min" : (-1,-1),
+            "excitations": bitstring.count("1"),
+            "landscape": np.zeros((len(time_list), len(delta_list))),
+            "max_prob": 0,
+            "min_prob": 0,
+            "location_max": (-1, -1),
+            "location_min": (-1, -1),
         }
 
-
     for j, d in enumerate(delta_list):
-        print(f"Simulating detuning {j+1} of {len(delta_list)}" if verbose else "", end = "\n" if verbose else "")
+        print(
+            f"Simulating detuning {j+1} of {len(delta_list)}" if verbose else "",
+            end="\n" if verbose else "",
+        )
 
-        print(f" Preprocessing shortcut: start" if verbose else "", end = "\n" if verbose else "")
+        print(
+            f" Preprocessing shortcut: start" if verbose else "",
+            end="\n" if verbose else "",
+        )
         T = time_list[-1]
         full_seq = Sequence(register, device)
         full_seq.declare_channel("ch", "rydberg_global")
@@ -589,15 +680,21 @@ def pulse_landscape_modulation_fast(register, device, time_buffer, independent_s
         full_res = full_sim.run(progress_bar=verbose)
         init_states = []
         for i, t in enumerate(time_list):
-            if t>=time_buffer:
-                init_states.append(full_res.get_state(t/1000))
-        print(f" Preprocessing shortcut: end" if verbose else "", end = "\n" if verbose else "")
+            if t >= time_buffer:
+                init_states.append(full_res.get_state(t / 1000))
+        print(
+            f" Preprocessing shortcut: end" if verbose else "",
+            end="\n" if verbose else "",
+        )
 
         k = 0
         for i, t in enumerate(time_list):
-            if t<time_buffer:
+            if t < time_buffer:
                 # independent simulation for each time only if you are below the time buffer
-                print(f"    Simulating time {i+1} of {len(time_list)}" if verbose else "", end = "\n" if verbose else "")
+                print(
+                    f"    Simulating time {i+1} of {len(time_list)}" if verbose else "",
+                    end="\n" if verbose else "",
+                )
                 seq = Sequence(register, device)
                 seq.declare_channel("ch", "rydberg_global")
                 if EOM:
@@ -611,10 +708,13 @@ def pulse_landscape_modulation_fast(register, device, time_buffer, independent_s
                 for bitstring in independent_sets:
                     obs = [observable(bitstring)]
                     exp = res.expect(obs)
-                    IS_dictionary[bitstring]["landscape"][i,j] = exp[0][-1]
+                    IS_dictionary[bitstring]["landscape"][i, j] = exp[0][-1]
             else:
                 # use shortcut for times above the time buffer
-                print(f"    Simulating time {i+1} of {len(time_list)}" if verbose else "", end = "\n" if verbose else "")
+                print(
+                    f"    Simulating time {i+1} of {len(time_list)}" if verbose else "",
+                    end="\n" if verbose else "",
+                )
                 seq = Sequence(register, device)
                 seq.declare_channel("ch", "rydberg_global")
                 pulse = pulse_tail
@@ -626,55 +726,71 @@ def pulse_landscape_modulation_fast(register, device, time_buffer, independent_s
                 for bitstring in independent_sets:
                     obs = [observable(bitstring)]
                     exp = res.expect(obs)
-                    IS_dictionary[bitstring]["landscape"][i,j] = exp[0][-1]
+                    IS_dictionary[bitstring]["landscape"][i, j] = exp[0][-1]
                 k += 1
 
     for bitstring in independent_sets:
-
-        IS_dictionary[bitstring]["max_prob"] = np.max(IS_dictionary[bitstring]["landscape"])
-        IS_dictionary[bitstring]["min_prob"] = np.min(IS_dictionary[bitstring]["landscape"])
-        IS_dictionary[bitstring]["location_max"] = (
-            np.unravel_index(
-                IS_dictionary[bitstring]["landscape"].argmax(),
-                IS_dictionary[bitstring]["landscape"].shape
-            )
+        IS_dictionary[bitstring]["max_prob"] = np.max(
+            IS_dictionary[bitstring]["landscape"]
         )
-        IS_dictionary[bitstring]["location_min"] = (
-            np.unravel_index(
-                IS_dictionary[bitstring]["landscape"].argmin(),
-                IS_dictionary[bitstring]["landscape"].shape
-            )
+        IS_dictionary[bitstring]["min_prob"] = np.min(
+            IS_dictionary[bitstring]["landscape"]
+        )
+        IS_dictionary[bitstring]["location_max"] = np.unravel_index(
+            IS_dictionary[bitstring]["landscape"].argmax(),
+            IS_dictionary[bitstring]["landscape"].shape,
+        )
+        IS_dictionary[bitstring]["location_min"] = np.unravel_index(
+            IS_dictionary[bitstring]["landscape"].argmin(),
+            IS_dictionary[bitstring]["landscape"].shape,
         )
 
     return IS_dictionary
 
-def plot_probability_landscape(projector_expvals, time_list, delta_list, max_p_clip=1, figsize=(8,4)):
+
+def plot_probability_landscape(
+    projector_expvals, time_list, delta_list, max_p_clip=1, figsize=(8, 4)
+):
     """Plots the projector probability over a time range defined by
     `time_list` and over a range of detuning values defined by `delta_list`"""
 
     for bitstring, data in projector_expvals.items():
         fig, ax = plt.subplots(figsize=figsize)
-        c = ax.pcolormesh(time_list, delta_list, data["landscape"].T, vmin=0, vmax=max_p_clip)
+        c = ax.pcolormesh(
+            time_list, delta_list, data["landscape"].T, vmin=0, vmax=max_p_clip
+        )
         ax.set_xlabel("time (ns)")
         ax.set_ylabel("detuning (rad/us)")
         ax.set_title(f"probability of bitstring: {bitstring}")
         fig.colorbar(c, ax=ax)
         plt.show()
         print(f"maximum probability: {data['max_prob']:.3f}")
-        print(f"probability of measuring: {1-(1-data['max_prob'])**200:.2f} (200 shots)")
-        i,j = data["location_max"]
+        print(
+            f"probability of measuring: {1-(1-data['max_prob'])**200:.2f} (200 shots)"
+        )
+        i, j = data["location_max"]
         print(f"time of max_prob: {time_list[i]}")
         print(f"detuning of max_prob: {delta_list[j]:.2f}")
         print()
         print()
 
-def plot_probability_contour(projector_expvals, time_list, delta_list, levels=8, max_p_clip=1, figsize=(8,4)):
+
+def plot_probability_contour(
+    projector_expvals, time_list, delta_list, levels=8, max_p_clip=1, figsize=(8, 4)
+):
     """Plots the projector probability over a time range defined by
     `time_list` and over a range of detuning values defined by `delta_list`"""
 
     for bitstring, data in projector_expvals.items():
         fig, ax = plt.subplots(figsize=figsize)
-        c = ax.contour(time_list, delta_list, data["landscape"].T, vmin=0, vmax=max_p_clip, levels=levels)
+        c = ax.contour(
+            time_list,
+            delta_list,
+            data["landscape"].T,
+            vmin=0,
+            vmax=max_p_clip,
+            levels=levels,
+        )
         ax.clabel(c, inline=True)
         ax.set_xlabel("time (ns)")
         ax.set_ylabel("detuning (rad/us)")
@@ -682,16 +798,19 @@ def plot_probability_contour(projector_expvals, time_list, delta_list, levels=8,
         fig.colorbar(c, ax=ax)
         plt.show()
         print(f"maximum probability: {data['max_prob']:.3f}")
-        print(f"probability of measuring: {1-(1-data['max_prob'])**200:.2f} (200 shots)")
-        i,j = data["location_max"]
+        print(
+            f"probability of measuring: {1-(1-data['max_prob'])**200:.2f} (200 shots)"
+        )
+        i, j = data["location_max"]
         print(f"time of max_prob: {time_list[i]}")
         print(f"detuning of max_prob: {delta_list[j]:.2f}")
         print()
         print()
 
+
 def stats_by_is_size(projector_expvals, n):
     """Computes projector probability statistics"""
-    for ex in range(1, n+1):
+    for ex in range(1, n + 1):
         print(f"{ex} excitations:")
         probs = []
         counter = 0
diff --git a/aquapointer/analog/utils/density_utils.py b/aquapointer/analog/utils/density_utils.py
index 3a7b3c6..3cbdf38 100644
--- a/aquapointer/analog/utils/density_utils.py
+++ b/aquapointer/analog/utils/density_utils.py
@@ -18,25 +18,27 @@ def rescale_register(pos, density, pad=0.2):
     so that they correspond as much as possible to the indices of the density.
     This means that the density value for the k-th qubit is found as:
     density[i,j] where i = int(rescaled_pos[k][1]), j = int(rescaled_pos[k][0])"""
-    reg_xlim = (np.min(pos[:,0]), np.max(pos[:,0]))
-    reg_ylim = (np.min(pos[:,1]), np.max(pos[:,1]))
-    reg_len = (reg_xlim[1]-reg_xlim[0], reg_ylim[1]-reg_ylim[0])
-    density_xlim = (0, density.shape[0]-1-pad)
-    density_ylim = (0, density.shape[1]-1-pad)
-    density_len = (density_xlim[1]-density_xlim[0], density_ylim[1]-density_ylim[0])
-    shift = (density_xlim[0]-reg_xlim[0], density_ylim[0]-reg_ylim[0])
-    scale = min(density_len[0]/reg_len[0], density_len[1]/reg_len[1])
-    return scale*(pos+shift)
+    reg_xlim = (np.min(pos[:, 0]), np.max(pos[:, 0]))
+    reg_ylim = (np.min(pos[:, 1]), np.max(pos[:, 1]))
+    reg_len = (reg_xlim[1] - reg_xlim[0], reg_ylim[1] - reg_ylim[0])
+    density_xlim = (0, density.shape[0] - 1 - pad)
+    density_ylim = (0, density.shape[1] - 1 - pad)
+    density_len = (density_xlim[1] - density_xlim[0], density_ylim[1] - density_ylim[0])
+    shift = (density_xlim[0] - reg_xlim[0], density_ylim[0] - reg_ylim[0])
+    scale = min(density_len[0] / reg_len[0], density_len[1] / reg_len[1])
+    return scale * (pos + shift)
+
 
 def density_of_qubit(k, density, rescaled_pos):
     """Returns the density value corresponding to the k-th qubit"""
     i = int(np.round(rescaled_pos[k][1]))
     j = int(np.round(rescaled_pos[k][0]))
-    return density[i,j]
+    return density[i, j]
+
 
 def move_blocks(positions, indices, movement, unit_vectors):
     """
-    `unit_vectors` is the directional unit vectors of the lattice 
+    `unit_vectors` is the directional unit vectors of the lattice
     `movement` is a tuple (nx,ny) where nx is the number of
     positions to move in direction `unit_vectors[0]`
     and similarly for ny
@@ -45,10 +47,20 @@ def move_blocks(positions, indices, movement, unit_vectors):
     for i in range(len(positions)):
         if i not in indices:
             continue
-        positions[i] += nx*unit_vectors[0] + ny*unit_vectors[1]
+        positions[i] += nx * unit_vectors[0] + ny * unit_vectors[1]
     return positions
 
-def decimate_register(full_pos, rescaled_full_pos, density, threshold, center=None, movement=None, indices=None, unit_vectors=None):
+
+def decimate_register(
+    full_pos,
+    rescaled_full_pos,
+    density,
+    threshold,
+    center=None,
+    movement=None,
+    indices=None,
+    unit_vectors=None,
+):
     """Given a full array of positions, keeps only the positions corresponding to a
     high enough value of the density. If center is given, then the register and
     the array of positions is shifted so that the central position is in `center`,
@@ -57,7 +69,7 @@ def decimate_register(full_pos, rescaled_full_pos, density, threshold, center=No
     pos = []
     rescaled_pos = []
     for k, node in enumerate(full_pos):
-        if density_of_qubit(k, density, rescaled_full_pos)/max_density > threshold:
+        if density_of_qubit(k, density, rescaled_full_pos) / max_density > threshold:
             pos.append(node)
             rescaled_pos.append(rescaled_full_pos[k])
     pos = np.array(pos)
@@ -67,16 +79,16 @@ def decimate_register(full_pos, rescaled_full_pos, density, threshold, center=No
     rescaled_pos = np.array(rescaled_pos)
     if center is not None:
         middle = select_middle_coordinate(pos, full_pos)
-        pos = pos-middle+center
+        pos = pos - middle + center
     reg = Register.from_coordinates(pos, center=False)
     return (reg, pos, rescaled_pos)
 
+
 def select_middle_coordinate(coords, full_coords):
     ideal_center = np.average(coords, axis=0)
     d = 1e10
     for c in full_coords:
-        if np.linalg.norm(c-ideal_center) < d:
-            d = np.linalg.norm(c-ideal_center)
+        if np.linalg.norm(c - ideal_center) < d:
+            d = np.linalg.norm(c - ideal_center)
             real_center = c
     return real_center
-
diff --git a/aquapointer/analog/utils/detuning_scale_utils.py b/aquapointer/analog/utils/detuning_scale_utils.py
index 4ed28be..70ef5db 100644
--- a/aquapointer/analog/utils/detuning_scale_utils.py
+++ b/aquapointer/analog/utils/detuning_scale_utils.py
@@ -21,7 +21,7 @@ def gaussian(var, m, x, y):
     distributions centered at `mean[0]`, `mean[1]`, ...
     and variance `var`
     """
-    return np.exp(-((x-m[0])**2 +(y-m[1])**2)/(2*var))/(2*np.pi*var)
+    return np.exp(-((x - m[0]) ** 2 + (y - m[1]) ** 2) / (2 * var)) / (2 * np.pi * var)
 
 
 def gaussian_mixture(shape, var, means):
@@ -29,7 +29,7 @@ def gaussian_mixture(shape, var, means):
     for i in range(len(res)):
         for j in range(len(res[0])):
             for mean in means:
-                res[j,i] += gaussian(var, mean, i, j)
+                res[j, i] += gaussian(var, mean, i, j)
     return res
 
 
diff --git a/aquapointer/digital/loaddata.py b/aquapointer/digital/loaddata.py
index a68a000..d803469 100644
--- a/aquapointer/digital/loaddata.py
+++ b/aquapointer/digital/loaddata.py
@@ -13,14 +13,16 @@
 PP_DIR = "/data/MUP1/MUP1_logfilter8_points/"
 REG_DIR = "/registers/"
 
-class LoadData:
 
+class LoadData:
     def __init__(self) -> None:
         self.d_list = [-1.0, -0.5, 0.0, 0.5, 1.0, 1.5]
         self.densities = self.load_density_slices(path=BASE_PATH + DENS_DIR)
         self.plane_points = self.load_plane_points(path=BASE_PATH + PP_DIR)
         self.register_positions = self.load_register_positions(path=BASE_PATH + REG_DIR)
-        self.rescaled_register_positions = self.load_rescaled_register_positions(path=BASE_PATH + REG_DIR)
+        self.rescaled_register_positions = self.load_rescaled_register_positions(
+            path=BASE_PATH + REG_DIR
+        )
 
     def load_density_slices(self, path: str) -> list[np.ndarray]:
         r"""The 3D-RISM density slices are saved as pickled files in the folder MUP1.
@@ -37,9 +39,9 @@ def load_density_slices(self, path: str) -> list[np.ndarray]:
         densities = []
         for d in self.d_list:
             filename = path + f"d{d}" + basename
-            with open(filename, 'rb') as file_in:
+            with open(filename, "rb") as file_in:
                 densities.append(pickle.load(file_in))
-                
+
         return densities
 
     def load_plane_points(self, path: str) -> list[np.ndarray]:
@@ -49,7 +51,7 @@ def load_plane_points(self, path: str) -> list[np.ndarray]:
 
         Args:
             path: Path to plane points files.
-        
+
         Returns:
             List of numpy arrays containing the plane points.
         """
@@ -57,9 +59,9 @@ def load_plane_points(self, path: str) -> list[np.ndarray]:
         points = []
         for d in self.d_list:
             filename = path + f"d{d}" + basename
-            with open(filename, 'rb') as file_in:
+            with open(filename, "rb") as file_in:
                 points.append(pickle.load(file_in))
-        
+
         return points
 
     def load_register_positions(self, path: str) -> list[np.ndarray]:
@@ -77,10 +79,10 @@ def load_register_positions(self, path: str) -> list[np.ndarray]:
         positions = []
         for i in range(len(self.d_list)):
             filename = path + basename + f"{i}.npy"
-            with open(filename, 'rb') as file_in:
+            with open(filename, "rb") as file_in:
                 pos = np.load(file_in)
             positions.append(pos)
-        
+
         return positions
 
     def load_rescaled_register_positions(self, path: str) -> list[np.ndarray]:
@@ -96,10 +98,9 @@ def load_rescaled_register_positions(self, path: str) -> list[np.ndarray]:
         basename = "rescaled_position_"
         rescaled_positions = []
         for i in range(len(self.d_list)):
-            filename = path + basename+f"{i}.npy"
-            with open(filename, 'rb') as file_in:
+            filename = path + basename + f"{i}.npy"
+            with open(filename, "rb") as file_in:
                 res_pos = np.load(file_in)
             rescaled_positions.append(res_pos)
-        
-        return rescaled_positions
 
+        return rescaled_positions
diff --git a/aquapointer/digital/qubo.py b/aquapointer/digital/qubo.py
index 8f7e289..67336b2 100644
--- a/aquapointer/digital/qubo.py
+++ b/aquapointer/digital/qubo.py
@@ -8,26 +8,33 @@
 from aquapointer.digital.loaddata import LoadData
 from aquapointer.digital.qubo_utils import gaussian, gaussian_mixture, gamma, Vij
 
-class Qubo:
 
+class Qubo:
     def __init__(self, loaddata: LoadData) -> None:
         self.ld = loaddata
-        self.qubo_matrices = self.get_qubo_matrices(densities=self.ld.densities, rescaled_positions=self.ld.rescaled_register_positions)
-        hamiltonians = [self.get_ising_hamiltonian(qubo=qubo) for qubo in self.qubo_matrices]
+        self.qubo_matrices = self.get_qubo_matrices(
+            densities=self.ld.densities,
+            rescaled_positions=self.ld.rescaled_register_positions,
+        )
+        hamiltonians = [
+            self.get_ising_hamiltonian(qubo=qubo) for qubo in self.qubo_matrices
+        ]
         self.qubo_hamiltonian_pairs = list(zip(self.qubo_matrices, hamiltonians))
 
-    def get_qubo_matrices(self, densities: list[np.ndarray], rescaled_positions: list[np.ndarray]) -> list[np.ndarray]:
-        r""" Given the density slices and rescaled positions of the registers,
+    def get_qubo_matrices(
+        self, densities: list[np.ndarray], rescaled_positions: list[np.ndarray]
+    ) -> list[np.ndarray]:
+        r"""Given the density slices and rescaled positions of the registers,
         one can compute the corresponding QUBO matrices.
 
         Args:
             densities: List of numpy arrays containing the 3D-RISM slices.
             rescaled_positions: List of numpy arrays containing the rescaled positions of the register.
-        
+
         Returns:
-            List of numpy arrays containing the QUBO matrices for each slice.        
-        
-        """        
+            List of numpy arrays containing the QUBO matrices for each slice.
+
+        """
         variance = 50
         amplitude = 6
         qubo_matrices = []
@@ -37,34 +44,36 @@ def get_qubo_matrices(self, densities: list[np.ndarray], rescaled_positions: lis
             # use function to calculate the one-body coefficients of the QUBO
             gamma_list = [gamma(density, pos, variance) for pos in rescaled_pos]
             qubo = np.zeros((len(rescaled_pos), len(rescaled_pos)))
-            
+
             for i in range(len(gamma_list)):
                 qubo[i][i] = -gamma_list[i]
 
             for i in range(len(rescaled_pos)):
-                for j in range(i+1, len(rescaled_pos)):
-                    qubo[i][j] = Vij(shape=density.shape, 
-                                    mean1=tuple(rescaled_pos[i]), 
-                                    mean2=tuple(rescaled_pos[j]), 
-                                    var=variance, 
-                                    amp=amplitude)
+                for j in range(i + 1, len(rescaled_pos)):
+                    qubo[i][j] = Vij(
+                        shape=density.shape,
+                        mean1=tuple(rescaled_pos[i]),
+                        mean2=tuple(rescaled_pos[j]),
+                        var=variance,
+                        amp=amplitude,
+                    )
                     qubo[j][i] = qubo[i][j]
             qubo_matrices.append(qubo)
-        
+
         return qubo_matrices
 
     def ising_energy(self, assignment: np.ndarray, qubo: np.ndarray) -> float:
-        r""" Given a binary string x and a QUBO matrix Q, computes the inner product <x, Qx>.
+        r"""Given a binary string x and a QUBO matrix Q, computes the inner product <x, Qx>.
 
-        Args:   
+        Args:
             assignment: numpy array, 0,1-valued.
             qubo: 2d numpy array.
 
         Returns:
             Float given by computing the inner product <x, Qx>.
-        
-        
-        """ 
+
+
+        """
         return np.transpose(assignment) @ qubo @ assignment
 
     def _bitfield(self, n: int, L: int) -> list[int]:
@@ -72,16 +81,16 @@ def _bitfield(self, n: int, L: int) -> list[int]:
         return [int(digit) for digit in result]
 
     def find_optimum(self, qubo: np.ndarray) -> tuple[str, float]:
-        r""" Brute-force approach to solving the QUBO problem: finding the optimal
+        r"""Brute-force approach to solving the QUBO problem: finding the optimal
         bitstring x that minimizes <x, Qx> where Q is the QUBO matrix.
 
         Args:
             qubo: 2d numpy array.
-        
+
         Returns:
-            Tuple of a bitstring and minimal energy (such that <x, Qx> is minimized).   
-        
-        """ 
+            Tuple of a bitstring and minimal energy (such that <x, Qx> is minimized).
+
+        """
         shape = qubo.shape
         L = shape[0]
 
@@ -93,12 +102,12 @@ def find_optimum(self, qubo: np.ndarray) -> tuple[str, float]:
             if energy < min_energy:
                 min_energy = energy
                 optimal_b = b
-        
-        sol = ''.join(map(str, optimal_b))        
+
+        sol = "".join(map(str, optimal_b))
         return sol, min_energy
 
     def _sparse_sigmaz_string(self, length: int, pos: list[int]) -> str:
-        r""" Given a list positions and integer length, returns a string of the
+        r"""Given a list positions and integer length, returns a string of the
         given length consisting of a Z at positions from the list of positions
         and I otherwise. This is interpreted as a sparse Pauli operator.
 
@@ -107,9 +116,9 @@ def _sparse_sigmaz_string(self, length: int, pos: list[int]) -> str:
             pos: List of integers for the positions of Z.
 
         Returns:
-            String consisting of I's and Z's.  
-                
-        """ 
+            String consisting of I's and Z's.
+
+        """
         sparse_sigmaz_str = ""
         for i in range(length):
             if i in pos:
@@ -119,34 +128,49 @@ def _sparse_sigmaz_string(self, length: int, pos: list[int]) -> str:
         return sparse_sigmaz_str
 
     def get_ising_hamiltonian(self, qubo: np.ndarray) -> SparsePauliOp:
-        r""" Given a QUBO matrix, one can associate with it a sparse Pauli operator.
+        r"""Given a QUBO matrix, one can associate with it a sparse Pauli operator.
         This is done by mapping a binary variable x -> z := (1-x)/2.
 
         Args:
             qubo: 2d numpy array.
-        
-        Returns:
-            SparsePauliOp corresponding to the QUBO matrix.       
-        
-        """ 
-        #the constant term (coefficient in front of II...I)
-        coeff_id = 0.5*np.sum([qubo[i][i] for i in range(len(qubo))])+0.5*np.sum([np.sum([qubo[i][j] for j in range(i+1,len(qubo))]) for i in range(len(qubo))])
-
-        #the linear terms (coefficient in front of I ... I sigma^z_i I ... I)
-        coeff_linear = [0.5*np.sum([qubo[i][j] for j in range(len(qubo))]) for i in range(len(qubo))]
 
-        #the quadratic terms (coefficient in front of I ... I sigma^z_i I ... I sigma^z_j I ... I)
-        coeff_quadratic = 0.25*qubo
-
-        #creat the list of sparse pauli operators and coefficients
+        Returns:
+            SparsePauliOp corresponding to the QUBO matrix.
+
+        """
+        # the constant term (coefficient in front of II...I)
+        coeff_id = 0.5 * np.sum([qubo[i][i] for i in range(len(qubo))]) + 0.5 * np.sum(
+            [
+                np.sum([qubo[i][j] for j in range(i + 1, len(qubo))])
+                for i in range(len(qubo))
+            ]
+        )
+
+        # the linear terms (coefficient in front of I ... I sigma^z_i I ... I)
+        coeff_linear = [
+            0.5 * np.sum([qubo[i][j] for j in range(len(qubo))])
+            for i in range(len(qubo))
+        ]
+
+        # the quadratic terms (coefficient in front of I ... I sigma^z_i I ... I sigma^z_j I ... I)
+        coeff_quadratic = 0.25 * qubo
+
+        # creat the list of sparse pauli operators and coefficients
         sparse_list = [(self._sparse_sigmaz_string(len(qubo), []), coeff_id)]
 
         for i in range(len(qubo)):
-            sparse_list.append((self._sparse_sigmaz_string(len(qubo), [i]), coeff_linear[i]))
+            sparse_list.append(
+                (self._sparse_sigmaz_string(len(qubo), [i]), coeff_linear[i])
+            )
             for j in range(len(qubo)):
                 if i != j:
-                    sparse_list.append((self._sparse_sigmaz_string(len(qubo), [i, j]), coeff_quadratic[i][j]))
+                    sparse_list.append(
+                        (
+                            self._sparse_sigmaz_string(len(qubo), [i, j]),
+                            coeff_quadratic[i][j],
+                        )
+                    )
 
         hamiltionian = SparsePauliOp.from_list(sparse_list)
 
-        return hamiltionian
\ No newline at end of file
+        return hamiltionian
diff --git a/aquapointer/digital/qubo_utils.py b/aquapointer/digital/qubo_utils.py
index 51981d2..43382c4 100644
--- a/aquapointer/digital/qubo_utils.py
+++ b/aquapointer/digital/qubo_utils.py
@@ -8,6 +8,7 @@
 from qiskit.quantum_info import SparsePauliOp
 from numba import njit
 
+
 @njit
 def gaussian(var, m, x, y):
     """
@@ -16,7 +17,8 @@ def gaussian(var, m, x, y):
     and variance `var`
     """
     res = 0
-    return np.exp(-((x-m[0])**2 +(y-m[1])**2)/(2*var))/(2*np.pi*var)
+    return np.exp(-((x - m[0]) ** 2 + (y - m[1]) ** 2) / (2 * var)) / (2 * np.pi * var)
+
 
 @njit
 def gaussian_mixture(shape, var, means):
@@ -24,41 +26,52 @@ def gaussian_mixture(shape, var, means):
     for i in range(len(res)):
         for j in range(len(res[0])):
             for mean in means:
-                res[j,i] += gaussian(var, mean, i, j)
+                res[j, i] += gaussian(var, mean, i, j)
     return res
 
+
 @njit
 def gamma(density, m, var, amp=1):
-    Nm = amp*gaussian_mixture(density.shape, var, [m])
+    Nm = amp * gaussian_mixture(density.shape, var, [m])
     res = 0
     for i in range(len(Nm[0])):
         for j in range(len(Nm)):
-            res += 2*Nm[i,j]*density[i,j]
-            res -= Nm[i,j]*Nm[i,j]
+            res += 2 * Nm[i, j] * density[i, j]
+            res -= Nm[i, j] * Nm[i, j]
     return res
 
+
 @njit
-def Vij(shape: tuple[int, int], mean1: tuple[float, float], mean2: tuple[float, float], var: float, amp: float) -> float:
-    Nm1 = amp*gaussian_mixture(shape, var, [mean1])
-    Nm2 = amp*gaussian_mixture(shape, var, [mean2])
+def Vij(
+    shape: tuple[int, int],
+    mean1: tuple[float, float],
+    mean2: tuple[float, float],
+    var: float,
+    amp: float,
+) -> float:
+    Nm1 = amp * gaussian_mixture(shape, var, [mean1])
+    Nm2 = amp * gaussian_mixture(shape, var, [mean2])
     res = 0
     for i in range(len(Nm1[0])):
         for j in range(len(Nm1)):
-            res += Nm1[i,j]*Nm2[i,j]
+            res += Nm1[i, j] * Nm2[i, j]
     return res
 
-def get_qubo_matrices(densities: list[np.ndarray], rescaled_positions: list[np.ndarray]) -> list[np.ndarray]:
-    r""" Given the density slices and rescaled positions of the registers,
+
+def get_qubo_matrices(
+    densities: list[np.ndarray], rescaled_positions: list[np.ndarray]
+) -> list[np.ndarray]:
+    r"""Given the density slices and rescaled positions of the registers,
     one can compute the corresponding QUBO matrices.
 
     Args:
         densities: List of numpy arrays containing the 3D-RISM slices.
         rescaled_positions: List of numpy arrays containing the rescaled positions of the register.
-    
+
     Returns:
-        List of numpy arrays containing the QUBO matrices for each slice.        
-    
-    """   
+        List of numpy arrays containing the QUBO matrices for each slice.
+
+    """
     variance = 50
     amplitude = 6
     qubo_matrices = []
@@ -68,53 +81,58 @@ def get_qubo_matrices(densities: list[np.ndarray], rescaled_positions: list[np.n
         # use function to calculate the one-body coefficients of the QUBO
         gamma_list = [gamma(density, pos, variance) for pos in rescaled_pos]
         qubo = np.zeros((len(rescaled_pos), len(rescaled_pos)))
-        
+
         for i in range(len(gamma_list)):
             qubo[i][i] = -gamma_list[i]
 
         for i in range(len(rescaled_pos)):
-            for j in range(i+1, len(rescaled_pos)):
-                #compute V_ij (the quadratic coefficients in the QUBO)
-                #if i!=j:
-                qubo[i][j] = Vij(shape=density.shape, 
-                                mean1=tuple(rescaled_pos[i]), 
-                                mean2=tuple(rescaled_pos[j]), 
-                                var=variance, 
-                                amp=amplitude)
+            for j in range(i + 1, len(rescaled_pos)):
+                # compute V_ij (the quadratic coefficients in the QUBO)
+                # if i!=j:
+                qubo[i][j] = Vij(
+                    shape=density.shape,
+                    mean1=tuple(rescaled_pos[i]),
+                    mean2=tuple(rescaled_pos[j]),
+                    var=variance,
+                    amp=amplitude,
+                )
                 qubo[j][i] = qubo[i][j]
         qubo_matrices.append(qubo)
-    
+
     return qubo_matrices
 
+
 def ising_energy(assignment: np.ndarray, qubo: np.ndarray) -> float:
-    r""" Given a binary string x and a QUBO matrix Q, computes the inner product <x, Qx>.
+    r"""Given a binary string x and a QUBO matrix Q, computes the inner product <x, Qx>.
 
-    Args:   
+    Args:
         assignment: numpy array, 0,1-valued.
         qubo: 2d numpy array.
 
     Returns:
         Float given by computing the inner product <x, Qx>.
-    
-    
-    """ 
+
+
+    """
     return np.transpose(assignment) @ qubo @ assignment
 
+
 def bitfield(n: int, L: int) -> list[int]:
     result = np.binary_repr(n, L)
     return [int(digit) for digit in result]
 
+
 def find_optimum(qubo: np.ndarray) -> tuple[str, float]:
-    r""" Brute-force approach to solving the QUBO problem: finding the optimal
+    r"""Brute-force approach to solving the QUBO problem: finding the optimal
     bitstring x that minimizes <x, Qx> where Q is the QUBO matrix.
 
     Args:
         qubo: 2d numpy array.
-    
+
     Returns:
-        Tuple of a bitstring and minimal energy (such that <x, Qx> is minimized).   
-    
-    """ 
+        Tuple of a bitstring and minimal energy (such that <x, Qx> is minimized).
+
+    """
     shape = qubo.shape
     L = shape[0]
 
@@ -126,7 +144,7 @@ def find_optimum(qubo: np.ndarray) -> tuple[str, float]:
         if energy < min_energy:
             min_energy = energy
             optimal_b = b
-    sol = ''.join(map(str, optimal_b))
+    sol = "".join(map(str, optimal_b))
     return sol, min_energy
 
 def sparse_sigmaz_string(length: int, pos: list[int]) -> str:
@@ -181,4 +199,4 @@ def get_ising_hamiltonian(qubo: np.ndarray) -> SparsePauliOp:
 
     hamiltionian = SparsePauliOp.from_list(sparse_list)
 
-    return hamiltionian
\ No newline at end of file
+    return hamiltionian
diff --git a/aquapointer/digital/vqe.py b/aquapointer/digital/vqe.py
index e1a9e74..a119a98 100644
--- a/aquapointer/digital/vqe.py
+++ b/aquapointer/digital/vqe.py
@@ -10,8 +10,16 @@
 from scipy.optimize import minimize
 from aquapointer.digital.qubo_utils import ising_energy
 
+
 class VQE:
-    def __init__(self, qubo: np.ndarray, ansatz: QuantumCircuit, sampler: BackendSampler, params: np.ndarray, prob_opt_sol: bool) -> None:
+    def __init__(
+        self,
+        qubo: np.ndarray,
+        ansatz: QuantumCircuit,
+        sampler: BackendSampler,
+        params: np.ndarray,
+        prob_opt_sol: bool,
+    ) -> None:
         self.qubo = qubo
         self.ansatz = ansatz
         self.sampler = sampler
@@ -19,29 +27,38 @@ def __init__(self, qubo: np.ndarray, ansatz: QuantumCircuit, sampler: BackendSam
         if params.any():
             self.params = params
         else:
-            self.params = np.array([np.random.random()]*self.ansatz.num_parameters)
+            self.params = np.array([np.random.random()] * self.ansatz.num_parameters)
 
         self.r = 0.1
         self.prob_opt_sol = prob_opt_sol
         self.history = []
 
     def run(self, alpha: float, maxiter: int, method="COBYLA"):
-        r""" Runs the minization.
+        r"""Runs the minization.
 
         Args:
             alpha: Confidence level.
             maxiter: Maximum number of iterations.
             method: Method for updating parameters.
-        
+
         Returns:
             Result from running scipy.optimize.minimize.
         """
-        res = minimize(self.cvar_energy, self.params, args=(alpha, ), method=method, tol=1e-8, options={"maxiter": maxiter})
+        res = minimize(
+            self.cvar_energy,
+            self.params,
+            args=(alpha,),
+            method=method,
+            tol=1e-8,
+            options={"maxiter": maxiter},
+        )
         self.params = res.x
         return res
-           
-    def _compute_cvar(self, probabilities: np.ndarray, values: np.ndarray, confidence_level: float) -> float:
-        r""" Compute Conditional Value at Risk (CVaR) for given probabilities, values, and confidence level.
+
+    def _compute_cvar(
+        self, probabilities: np.ndarray, values: np.ndarray, confidence_level: float
+    ) -> float:
+        r"""Compute Conditional Value at Risk (CVaR) for given probabilities, values, and confidence level.
 
         Args:
             probabilities: List or array of probabilities
@@ -62,21 +79,21 @@ def _compute_cvar(self, probabilities: np.ndarray, values: np.ndarray, confidenc
         exceed_index = np.argmax(cumulative_prob >= confidence_level)
 
         # Calculate CVaR
-        cvar_values = sorted_values[:exceed_index + 1]
-        cvar_probabilities = sorted_probabilities[:exceed_index + 1]
+        cvar_values = sorted_values[: exceed_index + 1]
+        cvar_probabilities = sorted_probabilities[: exceed_index + 1]
 
         cvar = np.sum(cvar_values * cvar_probabilities) / np.sum(cvar_probabilities)
 
         return cvar
-    
+
     def cvar_energy(self, params: np.ndarray, alpha: float) -> float:
-        r""" Function that takes parameters to bind to the ansatz and confidence level
+        r"""Function that takes parameters to bind to the ansatz and confidence level
         alpha, to compute the cvar energy (by sampling the ansatz and computing cvar).
-        
+
         Args:
             params: numpy array of parameters for the ansatz.
             alpha: Confidence level.
-        
+
         Returns:
             CVaR energy.
         """
@@ -87,13 +104,13 @@ def cvar_energy(self, params: np.ndarray, alpha: float) -> float:
         # Sample ansatz
         samp_dist = self.sampler.run(qc, shots=int(1e4)).result().quasi_dists[0]
 
-        samp_dist_binary=samp_dist.binary_probabilities()
+        samp_dist_binary = samp_dist.binary_probabilities()
 
         correct_dist = {}
         for key in samp_dist_binary.keys():
             reverse_key = key[::-1]
-            keynot = [(int(b)+1)%2 for b in reverse_key]
-            correct_dist[''.join(map(str, keynot))] = samp_dist_binary[key]
+            keynot = [(int(b) + 1) % 2 for b in reverse_key]
+            correct_dist["".join(map(str, keynot))] = samp_dist_binary[key]
 
         prob_energy = []
         bitstrings = []
@@ -128,7 +145,9 @@ def cvar_energy(self, params: np.ndarray, alpha: float) -> float:
         cvar_energy = self._compute_cvar(prob_energy[:, 0], prob_energy[:, 1], alpha)
 
         top_opt_prob = np.sum(sorted_probs[eps_rel_energies])
-        avg_top_energies = np.mean(sorted_probs[eps_rel_energies]*sorted_values[eps_rel_energies])
+        avg_top_energies = np.mean(
+            sorted_probs[eps_rel_energies] * sorted_values[eps_rel_energies]
+        )
 
         # save intermediate optimal bitsting and energy to self.history
         if self.prob_opt_sol:
@@ -136,6 +155,4 @@ def cvar_energy(self, params: np.ndarray, alpha: float) -> float:
         else:
             self.history.append([top_opt_prob, avg_top_energies])
 
-
-
-        return cvar_energy
\ No newline at end of file
+        return cvar_energy