diff --git a/README.md b/README.md
index 070ccb33..9f60ad32 100644
--- a/README.md
+++ b/README.md
@@ -32,7 +32,7 @@ import qibo
 
 # Below shows how to set the computation_settings
 # Note that for MPS_enabled and expectation_enabled parameters the accepted inputs are boolean or a dictionary with the format shown below.
-# If computation_settings is not specified, the default setting is used in which all booleans will be False. 
+# If computation_settings is not specified, the default setting is used in which all booleans will be False.
 # This will trigger the dense vector computation of the tensornet.
 
 computation_settings = {
@@ -92,4 +92,4 @@ Multi-node is enabled by setting either the MPI or NCCL enabled flag to True in
 
 ```sh
 mpirun -n 4 -hostfile $node_list python test.py
-```
\ No newline at end of file
+```
diff --git a/setup.py b/setup.py
index 13285eb7..0f619a5c 100644
--- a/setup.py
+++ b/setup.py
@@ -1,6 +1,7 @@
-from setuptools import setup, find_packages
-import re
 import pathlib
+import re
+
+from setuptools import find_packages, setup
 
 HERE = pathlib.Path(__file__).parent.absolute()
 PACKAGE = "qibotn"
@@ -8,8 +9,8 @@
 
 # Returns the qibotn version
 def version():
-    """Gets the version from the package's __init__ file
-    if there is some problem, let it happily fail"""
+    """Gets the version from the package's __init__ file if there is some
+    problem, let it happily fail."""
     version_file = HERE / "src" / PACKAGE / "__init__.py"
     version_regex = r"^__version__ = ['\"]([^'\"]*)['\"]"
 
diff --git a/src/qibotn/MPSUtils.py b/src/qibotn/MPSUtils.py
index 4f84f672..e8068f74 100644
--- a/src/qibotn/MPSUtils.py
+++ b/src/qibotn/MPSUtils.py
@@ -1,23 +1,19 @@
 import cupy as cp
-from cuquantum.cutensornet.experimental import contract_decompose
 from cuquantum import contract
+from cuquantum.cutensornet.experimental import contract_decompose
 
 # Reference: https://github.com/NVIDIA/cuQuantum/blob/main/python/samples/cutensornet/tn_algorithms/mps_algorithms.ipynb
 
 
 def initial(num_qubits, dtype):
-    """
-    Generate the MPS with an initial state of |00...00>
-    """
+    """Generate the MPS with an initial state of |00...00>"""
     state_tensor = cp.asarray([1, 0], dtype=dtype).reshape(1, 2, 1)
     mps_tensors = [state_tensor] * num_qubits
     return mps_tensors
 
 
 def mps_site_right_swap(mps_tensors, i, **kwargs):
-    """
-    Perform the swap operation between the ith and i+1th MPS tensors.
-    """
+    """Perform the swap operation between the ith and i+1th MPS tensors."""
     # contraction followed by QR decomposition
     a, _, b = contract_decompose(
         "ipj,jqk->iqj,jpk",
@@ -30,8 +26,7 @@ def mps_site_right_swap(mps_tensors, i, **kwargs):
 
 
 def apply_gate(mps_tensors, gate, qubits, **kwargs):
-    """
-    Apply the gate operand to the MPS tensors in-place.
+    """Apply the gate operand to the MPS tensors in-place.
 
     Args:
         mps_tensors: A list of rank-3 ndarray-like tensor objects.
diff --git a/src/qibotn/QiboCircuitConvertor.py b/src/qibotn/QiboCircuitConvertor.py
index e1aabea2..f67fb8ea 100644
--- a/src/qibotn/QiboCircuitConvertor.py
+++ b/src/qibotn/QiboCircuitConvertor.py
@@ -5,9 +5,9 @@
 
 
 class QiboCircuitToEinsum:
-    """Convert a circuit to a Tensor Network (TN) representation.
-    The circuit is first processed to an intermediate form by grouping each gate
-    matrix with its corresponding qubit it is acting on to a list. It is then
+    """Convert a circuit to a Tensor Network (TN) representation. The circuit
+    is first processed to an intermediate form by grouping each gate matrix
+    with its corresponding qubit it is acting on to a list. It is then
     converted to an equivalent TN expression through the class function
     state_vector_operands() following the Einstein summation convention in the
     interleave format.
@@ -79,9 +79,8 @@ def _parse_gates_to_mode_labels_operands(
         return mode_labels, operands
 
     def op_shape_from_qubits(self, nqubits):
-        """Modify tensor to cuQuantum shape
-        (qubit_states,input_output) * qubits_involved
-        """
+        """Modify tensor to cuQuantum shape (qubit_states,input_output) *
+        qubits_involved."""
         return (2, 2) * nqubits
 
     def init_intermediate_circuit(self, circuit):
@@ -134,8 +133,7 @@ def init_inverse_circuit(self, circuit):
         self.active_qubits_inverse = np.unique(gates_qubits_inverse)
 
     def get_pauli_gates(self, pauli_map, dtype="complex128", backend=cp):
-        """
-        Populate the gates for all pauli operators.
+        """Populate the gates for all pauli operators.
 
         Args:
             pauli_map: A dictionary mapping qubits to pauli operators.
diff --git a/src/qibotn/QiboCircuitToMPS.py b/src/qibotn/QiboCircuitToMPS.py
index 816b17c0..b1d847f2 100644
--- a/src/qibotn/QiboCircuitToMPS.py
+++ b/src/qibotn/QiboCircuitToMPS.py
@@ -1,9 +1,9 @@
 import cupy as cp
 import numpy as np
-
 from cuquantum import cutensornet as cutn
+
+from qibotn.MPSUtils import apply_gate, initial
 from qibotn.QiboCircuitConvertor import QiboCircuitToEinsum
-from qibotn.MPSUtils import initial, apply_gate
 
 
 class QiboCircuitToMPS:
diff --git a/src/qibotn/backends/__init__.py b/src/qibotn/backends/__init__.py
index f9279323..e5d68de5 100644
--- a/src/qibotn/backends/__init__.py
+++ b/src/qibotn/backends/__init__.py
@@ -1,2 +1,2 @@
-from qibotn.backends.gpu import CuTensorNet
 from qibotn.backends.cpu import QuTensorNet
+from qibotn.backends.gpu import CuTensorNet
diff --git a/src/qibotn/backends/cpu.py b/src/qibotn/backends/cpu.py
index 7115b395..a85dfd9b 100644
--- a/src/qibotn/backends/cpu.py
+++ b/src/qibotn/backends/cpu.py
@@ -1,8 +1,6 @@
-import numpy as np
-
 from qibo.backends.numpy import NumpyBackend
-from qibo.states import CircuitResult
 from qibo.config import raise_error
+from qibo.states import CircuitResult
 
 
 class QuTensorNet(NumpyBackend):
@@ -60,7 +58,6 @@ def execute_circuit(
 
         Returns:
             xxx.
-
         """
 
         import qibotn.eval_qu as eval
diff --git a/src/qibotn/backends/gpu.py b/src/qibotn/backends/gpu.py
index 5777fe90..2c3f8d41 100644
--- a/src/qibotn/backends/gpu.py
+++ b/src/qibotn/backends/gpu.py
@@ -1,8 +1,7 @@
 import numpy as np
-
 from qibo.backends.numpy import NumpyBackend
-from qibo.states import CircuitResult
 from qibo.config import raise_error
+from qibo.states import CircuitResult
 
 
 class CuTensorNet(NumpyBackend):  # pragma: no cover
@@ -107,7 +106,6 @@ def execute_circuit(
 
         Returns:
             xxx.
-
         """
 
         import qibotn.eval as eval
diff --git a/src/qibotn/eval.py b/src/qibotn/eval.py
index 96fd488a..5fcb66f0 100644
--- a/src/qibotn/eval.py
+++ b/src/qibotn/eval.py
@@ -1,20 +1,22 @@
-from qibotn.QiboCircuitConvertor import QiboCircuitToEinsum
-from cuquantum import contract
-from cupy.cuda.runtime import getDeviceCount
 import cupy as cp
+from cupy.cuda.runtime import getDeviceCount
+from cuquantum import contract
 
-from qibotn.QiboCircuitToMPS import QiboCircuitToMPS
 from qibotn.mps_contraction_helper import MPSContractionHelper
+from qibotn.QiboCircuitConvertor import QiboCircuitToEinsum
+from qibotn.QiboCircuitToMPS import QiboCircuitToMPS
 
 
 def dense_vector_tn(qibo_circ, datatype):
-    """Convert qibo circuit to tensornet (TN) format and perform contraction to dense vector."""
+    """Convert qibo circuit to tensornet (TN) format and perform contraction to
+    dense vector."""
     myconvertor = QiboCircuitToEinsum(qibo_circ, dtype=datatype)
     return contract(*myconvertor.state_vector_operands())
 
 
 def expectation_pauli_tn(qibo_circ, datatype, pauli_string_pattern):
-    """Convert qibo circuit to tensornet (TN) format and perform contraction to expectation of given Pauli string."""
+    """Convert qibo circuit to tensornet (TN) format and perform contraction to
+    expectation of given Pauli string."""
     myconvertor = QiboCircuitToEinsum(qibo_circ, dtype=datatype)
     return contract(
         *myconvertor.expectation_operands(
@@ -24,14 +26,19 @@ def expectation_pauli_tn(qibo_circ, datatype, pauli_string_pattern):
 
 
 def dense_vector_tn_MPI(qibo_circ, datatype, n_samples=8):
-    """Convert qibo circuit to tensornet (TN) format and perform contraction using multi node and multi GPU through MPI.
-    The conversion is performed by QiboCircuitToEinsum(), after which it goes through 2 steps: pathfinder and execution.
-    The pathfinder looks at user defined number of samples (n_samples) iteratively to select the least costly contraction path. This is sped up with multi thread.
-    After pathfinding the optimal path is used in the actual contraction to give a dense vector representation of the TN.
+    """Convert qibo circuit to tensornet (TN) format and perform contraction
+    using multi node and multi GPU through MPI.
+
+    The conversion is performed by QiboCircuitToEinsum(), after which it
+    goes through 2 steps: pathfinder and execution. The pathfinder looks
+    at user defined number of samples (n_samples) iteratively to select
+    the least costly contraction path. This is sped up with multi
+    thread. After pathfinding the optimal path is used in the actual
+    contraction to give a dense vector representation of the TN.
     """
 
-    from mpi4py import MPI
     from cuquantum import Network
+    from mpi4py import MPI
 
     root = 0
     comm = MPI.COMM_WORLD
@@ -86,14 +93,19 @@ def dense_vector_tn_MPI(qibo_circ, datatype, n_samples=8):
 
 
 def dense_vector_tn_nccl(qibo_circ, datatype, n_samples=8):
-    """Convert qibo circuit to tensornet (TN) format and perform contraction using multi node and multi GPU through NCCL.
-    The conversion is performed by QiboCircuitToEinsum(), after which it goes through 2 steps: pathfinder and execution.
-    The pathfinder looks at user defined number of samples (n_samples) iteratively to select the least costly contraction path. This is sped up with multi thread.
-    After pathfinding the optimal path is used in the actual contraction to give a dense vector representation of the TN.
+    """Convert qibo circuit to tensornet (TN) format and perform contraction
+    using multi node and multi GPU through NCCL.
+
+    The conversion is performed by QiboCircuitToEinsum(), after which it
+    goes through 2 steps: pathfinder and execution. The pathfinder looks
+    at user defined number of samples (n_samples) iteratively to select
+    the least costly contraction path. This is sped up with multi
+    thread. After pathfinding the optimal path is used in the actual
+    contraction to give a dense vector representation of the TN.
     """
-    from mpi4py import MPI
-    from cuquantum import Network
     from cupy.cuda import nccl
+    from cuquantum import Network
+    from mpi4py import MPI
 
     root = 0
     comm_mpi = MPI.COMM_WORLD
@@ -159,15 +171,22 @@ def dense_vector_tn_nccl(qibo_circ, datatype, n_samples=8):
 
 
 def expectation_pauli_tn_nccl(qibo_circ, datatype, pauli_string_pattern, n_samples=8):
-    """Convert qibo circuit to tensornet (TN) format and perform contraction to expectation of given Pauli string using multi node and multi GPU through NCCL.
-    The conversion is performed by QiboCircuitToEinsum(), after which it goes through 2 steps: pathfinder and execution.
-    The pauli_string_pattern is used to generate the pauli string corresponding to the number of qubits of the system.
-    The pathfinder looks at user defined number of samples (n_samples) iteratively to select the least costly contraction path. This is sped up with multi thread.
-    After pathfinding the optimal path is used in the actual contraction to give an expectation value.
+    """Convert qibo circuit to tensornet (TN) format and perform contraction to
+    expectation of given Pauli string using multi node and multi GPU through
+    NCCL.
+
+    The conversion is performed by QiboCircuitToEinsum(), after which it
+    goes through 2 steps: pathfinder and execution. The
+    pauli_string_pattern is used to generate the pauli string
+    corresponding to the number of qubits of the system. The pathfinder
+    looks at user defined number of samples (n_samples) iteratively to
+    select the least costly contraction path. This is sped up with multi
+    thread. After pathfinding the optimal path is used in the actual
+    contraction to give an expectation value.
     """
-    from mpi4py import MPI
-    from cuquantum import Network
     from cupy.cuda import nccl
+    from cuquantum import Network
+    from mpi4py import MPI
 
     root = 0
     comm_mpi = MPI.COMM_WORLD
@@ -235,14 +254,21 @@ def expectation_pauli_tn_nccl(qibo_circ, datatype, pauli_string_pattern, n_sampl
 
 
 def expectation_pauli_tn_MPI(qibo_circ, datatype, pauli_string_pattern, n_samples=8):
-    """Convert qibo circuit to tensornet (TN) format and perform contraction to expectation of given Pauli string using multi node and multi GPU through MPI.
-    The conversion is performed by QiboCircuitToEinsum(), after which it goes through 2 steps: pathfinder and execution.
-    The pauli_string_pattern is used to generate the pauli string corresponding to the number of qubits of the system.
-    The pathfinder looks at user defined number of samples (n_samples) iteratively to select the least costly contraction path. This is sped up with multi thread.
-    After pathfinding the optimal path is used in the actual contraction to give an expectation value.
+    """Convert qibo circuit to tensornet (TN) format and perform contraction to
+    expectation of given Pauli string using multi node and multi GPU through
+    MPI.
+
+    The conversion is performed by QiboCircuitToEinsum(), after which it
+    goes through 2 steps: pathfinder and execution. The
+    pauli_string_pattern is used to generate the pauli string
+    corresponding to the number of qubits of the system. The pathfinder
+    looks at user defined number of samples (n_samples) iteratively to
+    select the least costly contraction path. This is sped up with multi
+    thread. After pathfinding the optimal path is used in the actual
+    contraction to give an expectation value.
     """
-    from mpi4py import MPI  # this line initializes MPI
     from cuquantum import Network
+    from mpi4py import MPI  # this line initializes MPI
 
     root = 0
     comm = MPI.COMM_WORLD
@@ -299,7 +325,8 @@ def expectation_pauli_tn_MPI(qibo_circ, datatype, pauli_string_pattern, n_sample
 
 
 def dense_vector_mps(qibo_circ, gate_algo, datatype):
-    """Convert qibo circuit to matrix product state (MPS) format and perform contraction to dense vector."""
+    """Convert qibo circuit to matrix product state (MPS) format and perform
+    contraction to dense vector."""
     myconvertor = QiboCircuitToMPS(qibo_circ, gate_algo, dtype=datatype)
     mps_helper = MPSContractionHelper(myconvertor.num_qubits)
 
@@ -309,7 +336,9 @@ def dense_vector_mps(qibo_circ, gate_algo, datatype):
 
 
 def pauli_string_gen(nqubits, pauli_string_pattern):
-    """Used internally to generate the string based on given pattern and number of qubit.
+    """Used internally to generate the string based on given pattern and number
+    of qubit.
+
     Example: pattern: "XZ", number of qubit: 7, output = XZXZXZX
     """
     if nqubits <= 0:
diff --git a/src/qibotn/eval_qu.py b/src/qibotn/eval_qu.py
index 579a42ae..7b603b50 100644
--- a/src/qibotn/eval_qu.py
+++ b/src/qibotn/eval_qu.py
@@ -3,9 +3,15 @@
 from qibo.models import Circuit as QiboCircuit
 
 
-def from_qibo(circuit: QiboCircuit, is_mps: False, psi0=None, method='svd',
-              cutoff=1e-6, cutoff_mode='abs'):
-    """Create a tensornetwork representation of the circuit"""
+def from_qibo(
+    circuit: QiboCircuit,
+    is_mps: False,
+    psi0=None,
+    method="svd",
+    cutoff=1e-6,
+    cutoff_mode="abs",
+):
+    """Create a tensornetwork representation of the circuit."""
 
     nqubits = circuit.nqubits
     gate_opt = {}
@@ -30,19 +36,17 @@ def from_qibo(circuit: QiboCircuit, is_mps: False, psi0=None, method='svd',
 
 
 def init_state_tn(nqubits, init_state_sv):
-
-    """Create a matrixproductstate directly from a dense vector"""
+    """Create a matrixproductstate directly from a dense vector."""
 
     dims = tuple(2 * np.ones(nqubits, dtype=int))
 
     return qtn.tensor_1d.MatrixProductState.from_dense(init_state_sv, dims)
 
 
-def dense_vector_tn_qu(qasm: str, initial_state, is_mps,  backend="numpy"):
-    """Evaluate QASM with Quimb
+def dense_vector_tn_qu(qasm: str, initial_state, is_mps, backend="numpy"):
+    """Evaluate QASM with Quimb.
 
     backend (quimb): numpy, cupy, jax. Passed to ``opt_einsum``.
-
     """
     circuit = QiboCircuit.from_qasm(qasm)
     if initial_state is not None:
diff --git a/src/qibotn/mps_contraction_helper.py b/src/qibotn/mps_contraction_helper.py
index 29d5e253..1c004de4 100644
--- a/src/qibotn/mps_contraction_helper.py
+++ b/src/qibotn/mps_contraction_helper.py
@@ -1,11 +1,10 @@
-from cuquantum import contract, contract_path, CircuitToEinsum, tensor
+from cuquantum import contract, contract_path
 
 # Reference: https://github.com/NVIDIA/cuQuantum/blob/main/python/samples/cutensornet/tn_algorithms/mps_algorithms.ipynb
 
 
 class MPSContractionHelper:
-    """
-    A helper class to compute various quantities for a given MPS.
+    """A helper class to compute various quantities for a given MPS.
 
     Interleaved format is used to construct the input args for `cuquantum.contract`.
     A concrete example on how the modes are populated for a 7-site MPS is provided below:
@@ -43,8 +42,8 @@ def __init__(self, num_qubits):
         ]
 
     def contract_norm(self, mps_tensors, options=None):
-        """
-        Contract the corresponding tensor network to form the norm of the MPS.
+        """Contract the corresponding tensor network to form the norm of the
+        MPS.
 
         Args:
             mps_tensors: A list of rank-3 ndarray-like tensor objects.
@@ -64,8 +63,8 @@ def contract_norm(self, mps_tensors, options=None):
         return self._contract(interleaved_inputs, options=options).real
 
     def contract_state_vector(self, mps_tensors, options=None):
-        """
-        Contract the corresponding tensor network to form the state vector representation of the MPS.
+        """Contract the corresponding tensor network to form the state vector
+        representation of the MPS.
 
         Args:
             mps_tensors: A list of rank-3 ndarray-like tensor objects.
@@ -86,8 +85,8 @@ def contract_state_vector(self, mps_tensors, options=None):
     def contract_expectation(
         self, mps_tensors, operator, qubits, options=None, normalize=False
     ):
-        """
-        Contract the corresponding tensor network to form the expectation of the MPS.
+        """Contract the corresponding tensor network to form the expectation of
+        the MPS.
 
         Args:
             mps_tensors: A list of rank-3 ndarray-like tensor objects.
diff --git a/tests/test_cuquantum_cutensor_backend.py b/tests/test_cuquantum_cutensor_backend.py
index 57b3f342..c8f1e199 100644
--- a/tests/test_cuquantum_cutensor_backend.py
+++ b/tests/test_cuquantum_cutensor_backend.py
@@ -1,8 +1,8 @@
 from timeit import default_timer as timer
 
 import config
-import numpy as np
 import cupy as cp
+import numpy as np
 import pytest
 import qibo
 from qibo.models import QFT
diff --git a/tests/test_quimb_backend.py b/tests/test_quimb_backend.py
index 81a0e2b7..15ba6524 100644
--- a/tests/test_quimb_backend.py
+++ b/tests/test_quimb_backend.py
@@ -1,5 +1,6 @@
 import copy
 import os
+
 import config
 import numpy as np
 import pytest
@@ -8,8 +9,7 @@
 
 
 def create_init_state(nqubits):
-    init_state = np.random.random(2**nqubits) + \
-        1j * np.random.random(2**nqubits)
+    init_state = np.random.random(2**nqubits) + 1j * np.random.random(2**nqubits)
     init_state = init_state / np.sqrt((np.abs(init_state) ** 2).sum())
     return init_state
 
@@ -20,10 +20,11 @@ def qibo_qft(nqubits, init_state, swaps):
     return circ_qibo, state_vec
 
 
-@pytest.mark.parametrize("nqubits, tolerance, is_mps",
-                         [(1, 1e-6, True), (2, 1e-6, False), (5, 1e-3, True), (10, 1e-3, False)])
+@pytest.mark.parametrize(
+    "nqubits, tolerance, is_mps",
+    [(1, 1e-6, True), (2, 1e-6, False), (5, 1e-3, True), (10, 1e-3, False)],
+)
 def test_eval(nqubits: int, tolerance: float, is_mps: bool):
-
     """Evaluate circuit with Quimb backend.
 
     Args:
@@ -41,20 +42,18 @@ def test_eval(nqubits: int, tolerance: float, is_mps: bool):
     init_state_tn = copy.deepcopy(init_state)
 
     # Test qibo
-    qibo.set_backend(backend=config.qibo.backend,
-                     platform=config.qibo.platform)
-   
-    qibo_circ, result_sv= qibo_qft(nqubits, init_state, swaps=True)
-    
+    qibo.set_backend(backend=config.qibo.backend, platform=config.qibo.platform)
+
+    qibo_circ, result_sv = qibo_qft(nqubits, init_state, swaps=True)
 
     # Convert to qasm for other backends
     qasm_circ = qibo_circ.to_qasm()
 
     # Test quimb
     result_tn = qibotn.eval_qu.dense_vector_tn_qu(
-            qasm_circ, init_state_tn, is_mps, backend=config.quimb.backend
-        ).flatten()
-   
+        qasm_circ, init_state_tn, is_mps, backend=config.quimb.backend
+    ).flatten()
 
-    assert np.allclose(result_sv, result_tn,
-                       atol=tolerance), "Resulting dense vectors do not match"
+    assert np.allclose(
+        result_sv, result_tn, atol=tolerance
+    ), "Resulting dense vectors do not match"