diff --git a/doc/source/api-reference/qibo.rst b/doc/source/api-reference/qibo.rst
index 1b254088b3..d877c4c28d 100644
--- a/doc/source/api-reference/qibo.rst
+++ b/doc/source/api-reference/qibo.rst
@@ -473,7 +473,8 @@ factor results in the mitigated Pauli expectation value :math:`\langle O\rangle_
 .. math::
    \langle O\rangle_{ideal} = \frac{\langle O\rangle_{noisy}}{\lambda}
 
-.. autofunction:: qibo.models.error_mitigation.apply_randomized_readout_mitigation
+This process can be implemented with the aforementioned
+:func:`qibo.models.error_mitigation.apply_randomized_readout_mitigation`.
 
 
 Zero Noise Extrapolation (ZNE)
@@ -2382,7 +2383,7 @@ Hellinger fidelity
 Hellinger shot error
 """"""""""""""""""""
 
-.. autofunction:: qibo.quantum_info.hellinger_fidelity
+.. autofunction:: qibo.quantum_info.hellinger_shot_error
 
 
 Haar integral
diff --git a/doc/source/code-examples/examples.rst b/doc/source/code-examples/examples.rst
index 1896cffb51..081e8fc71d 100644
--- a/doc/source/code-examples/examples.rst
+++ b/doc/source/code-examples/examples.rst
@@ -347,8 +347,6 @@ For example, we can draw the QFT circuit for 5-qubits:
         # new plot function based on matplotlib
         from qibo.ui import plot_circuit
 
-        %matplotlib inline
-
         # create a 5-qubits QFT circuit
         c = QFT(5)
         c.add(gates.M(qubit) for qubit in range(2))
diff --git a/examples/shor/functions.py b/examples/shor/functions.py
index f4b2b36646..f36f2b507a 100644
--- a/examples/shor/functions.py
+++ b/examples/shor/functions.py
@@ -393,8 +393,9 @@ def quantum_order_finding_semiclassical(N, a):
         circuit.add(gates.U1(q_reg, -angle))
         circuit.add(gates.H(q_reg))
         results.append(circuit.add(gates.M(q_reg, collapse=True)))
+    circuit.add(gates.M(q_reg))
 
-    circuit()  # execute
+    circuit(nshots=1)  # execute
     s = sum(int(r.symbols[0].outcome()) * (2**i) for i, r in enumerate(results))
     print(f"The quantum circuit measures s = {s}.\n")
     return s
diff --git a/src/qibo/backends/abstract.py b/src/qibo/backends/abstract.py
index 5a14c24025..abe1f20d83 100644
--- a/src/qibo/backends/abstract.py
+++ b/src/qibo/backends/abstract.py
@@ -193,7 +193,7 @@ def execute_circuit(
     def execute_circuits(
         self, circuits, initial_states=None, nshots=None
     ):  # pragma: no cover
-        """Execute multiple :class:`qibo.models.circuit.Circuit`s in parallel."""
+        """Execute multiple :class:`qibo.models.circuit.Circuit` in parallel."""
         raise_error(NotImplementedError)
 
     @abc.abstractmethod
diff --git a/src/qibo/backends/clifford.py b/src/qibo/backends/clifford.py
index 1ba071cd40..aee8750e34 100644
--- a/src/qibo/backends/clifford.py
+++ b/src/qibo/backends/clifford.py
@@ -90,13 +90,13 @@ def calculate_frequencies(self, samples):
         return collections.Counter(dict(zip(res, counts)))
 
     def zero_state(self, nqubits: int):
-        """Construct the zero state |00...00>.
+        """Construct the zero state :math`\\ket{00...00}`.
 
         Args:
-            nqubits (int): Number of qubits.
+            nqubits (int): number of qubits.
 
         Returns:
-            (ndarray): Symplectic matrix for the zero state.
+            ndarray: Symplectic matrix for the zero state.
         """
         identity = self.np.eye(nqubits)
         symplectic_matrix = self.np.zeros(
diff --git a/src/qibo/gates/channels.py b/src/qibo/gates/channels.py
index bde8c50273..be43678318 100644
--- a/src/qibo/gates/channels.py
+++ b/src/qibo/gates/channels.py
@@ -48,8 +48,11 @@ def to_choi(self, nqubits: Optional[int] = None, order: str = "row", backend=Non
         of the Kraus channel :math:`\\{K_{\\alpha}\\}_{\\alpha}`.
 
         .. math::
-            \\mathcal{E} = \\sum_{\\alpha} \\, |K_{\\alpha}\\rangle\\rangle
-                \\langle\\langle K_{\\alpha}|
+            \\mathcal{E} = \\sum_{\\alpha} \\, |K_{\\alpha})(K_{\\alpha}| \\, ,
+
+        where :math:`|K_{\\alpha})` is the vectorization of the Kraus operator
+        :math:`K_{\\alpha}`.
+        For a definition of vectorization, see :func:`qibo.quantum_info.vectorization`.
 
         Args:
             nqubits (int, optional): total number of qubits to be considered
diff --git a/src/qibo/gates/measurements.py b/src/qibo/gates/measurements.py
index 7e1559e9d9..64a7a98e50 100644
--- a/src/qibo/gates/measurements.py
+++ b/src/qibo/gates/measurements.py
@@ -17,27 +17,26 @@ class M(Gate):
             If the qubits to measure are held in an iterable (eg. list) the ``*``
             operator can be used, for example ``gates.M(*[0, 1, 4])`` or
             ``gates.M(*range(5))``.
-        register_name (str): Optional name of the register to distinguish it
+        register_name (str, optional): Optional name of the register to distinguish it
             from other registers when used in circuits.
         collapse (bool): Collapse the state vector after the measurement is
             performed. Can be used only for single shot measurements.
             If ``True`` the collapsed state vector is returned. If ``False``
             the measurement result is returned.
-        basis (:class:`qibo.gates.Gate`, str, list): Basis to measure.
+        basis (:class:`qibo.gates.Gate` or str or list, optional): Basis to measure.
             Can be either:
             - a qibo gate
             - the string representing the gate
             - a callable that accepts a qubit, for example: ``lambda q: gates.RX(q, 0.2)``
-            - a list of the above, if a different basis will be used for each
-              measurement qubit.
-            Default is Z.
-        p0 (dict): Optional bitflip probability map. Can be:
+            - a list of the above, if a different basis will be used for each measurement qubit.
+            Defaults is to :class:`qibo.gates.Z`.
+        p0 (dict, optional): bitflip probability map. Can be:
             A dictionary that maps each measured qubit to the probability
             that it is flipped, a list or tuple that has the same length
             as the tuple of measured qubits or a single float number.
             If a single float is given the same probability will be used
             for all qubits.
-        p1 (dict): Optional bitflip probability map for asymmetric bitflips.
+        p1 (dict, optional): bitflip probability map for asymmetric bitflips.
             Same as ``p0`` but controls the 1->0 bitflip probability.
             If ``p1`` is ``None`` then ``p0`` will be used both for 0->1 and
             1->0 bitflips.
diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index eea1071399..3b16f8ac9d 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -658,8 +658,8 @@ def apply_randomized_readout_mitigation(
         nshots (int, optional): number of shots. Defaults to :math:`10000`.
         ncircuits (int, optional): number of randomized circuits. Each of them uses
             ``int(nshots / ncircuits)`` shots. Defaults to 10.
-        qubit_map (list, optional): the qubit map. If None, a list of range of circuit's qubits is used.
-            Defaults to ``None``.
+        qubit_map (list, optional): the qubit map. If ``None``, a list of range of circuit's
+            qubits is used. Defaults to ``None``.
         seed (int or :class:`numpy.random.Generator`, optional): Either a generator of random
             numbers or a fixed seed to initialize a generator. If ``None``, initializes
             a generator with a random seed. Default: ``None``.
@@ -667,15 +667,15 @@ def apply_randomized_readout_mitigation(
             in the execution. If ``None``, it uses :class:`qibo.backends.GlobalBackend`.
             Defaults to ``None``.
 
-    Return:
+    Returns:
         :class:`qibo.measurements.CircuitResult`: the state of the input circuit with
             mitigated frequencies.
 
 
-    Reference:
+    References:
         1. Ewout van den Berg, Zlatko K. Minev et al,
         *Model-free readout-error mitigation for quantum expectation values*.
-           `arXiv:2012.09738 [quant-ph] <https://arxiv.org/abs/2012.09738>`_.
+        `arXiv:2012.09738 [quant-ph] <https://arxiv.org/abs/2012.09738>`_.
     """
     from qibo import Circuit  # pylint: disable=import-outside-toplevel
     from qibo.quantum_info import (  # pylint: disable=import-outside-toplevel
diff --git a/src/qibo/quantum_info/basis.py b/src/qibo/quantum_info/basis.py
index 5ea541bca0..e3bf042a60 100644
--- a/src/qibo/quantum_info/basis.py
+++ b/src/qibo/quantum_info/basis.py
@@ -28,15 +28,15 @@ def pauli_basis(
         vectorize (bool, optional): If ``False``, returns a nested array with
             all Pauli matrices. If ``True``, retuns an array where every
             row is a vectorized Pauli matrix. Defaults to ``False``.
-        sparse (bool, optional) If ``True``, retuns Pauli basis in a sparse
-            representation. Default is ``False``.
+        sparse (bool, optional): If ``True``, retuns Pauli basis in a sparse
+            representation. Defaults to ``False``.
         order (str, optional): If ``"row"``, vectorization of Pauli basis is
             performed row-wise. If ``"column"``, vectorization is performed
             column-wise. If ``"system"``, system-wise vectorization is
             performed. If ``vectorization=False``, then ``order=None`` is
-            forced. Default is ``None``.
+            forced. Defaults to ``None``.
         pauli_order (str, optional): corresponds to the order of 4 single-qubit
-            Pauli elements. Default is "IXYZ".
+            Pauli elements. Defaults to ``"IXYZ"``.
         backend (:class:`qibo.backends.abstract.Backend`, optional): backend
             to be used in the execution. If ``None``, it uses
             :class:`qibo.backends.GlobalBackend`. Defaults to ``None``.
@@ -148,15 +148,12 @@ def comp_basis_to_pauli(
     The unitary :math:`U` is given by
 
     .. math::
-        U = \\sum_{k = 0}^{d^{2} - 1} \\, \\ketbra{k}{P_{k}} \\,\\, ,
+        U = \\sum_{k = 0}^{d^{2} - 1} \\, |k)(P_{k}| \\,\\, ,
 
-    where :math:`\\ket{P_{k}}` is the system-vectorization of the :math:`k`-th
-    Pauli operator :math:`P_{k}`, and :math:`\\ket{k}` is the computational
-    basis element.
-
-    When converting a state :math:`\\ket{\\rho}` to its Pauli-Liouville
-    representation :math:`\\ket{\\rho'}`, one should use ``order="system"``
-    in :func:`vectorization`.
+    where :math:`|P_{k})` is the vectorization of the :math:`k`-th
+    Pauli operator :math:`P_{k}`, and :math:`|k)` is the vectorization
+    of the :math:`k`-th computational basis element.
+    For a definition of vectorization, see :func:`qibo.quantum_info.vectorization`.
 
     Example:
         .. code-block:: python
@@ -174,13 +171,13 @@ def comp_basis_to_pauli(
         normalize (bool, optional): If ``True``, converts to the
             Pauli basis. Defaults to ``False``.
         sparse (bool, optional): If ``True``, returns unitary matrix in
-            sparse representation. Default is ``False``.
+            sparse representation. Defaults to ``False``.
         order (str, optional): If ``"row"``, vectorization of Pauli basis is
             performed row-wise. If ``"column"``, vectorization is performed
             column-wise. If ``"system"``, system-wise vectorization is
-            performed. Default is ``"row"``.
+            performed. Defaults to ``"row"``.
         pauli_order (str, optional): corresponds to the order of 4 single-qubit
-            Pauli elements. Default is "IXYZ".
+            Pauli elements. Defaults to ``"IXYZ"``.
         backend (:class:`qibo.backends.abstract.Backend`, optional): backend to be
             used in the execution. If ``None``, it uses
             :class:`qibo.backends.GlobalBackend`. Defaults to ``None``.
@@ -237,20 +234,25 @@ def pauli_to_comp_basis(
     The unitary :math:`U` is given by
 
     .. math::
-        U = \\sum_{k = 0}^{d^{2} - 1} \\, \\ketbra{P_{k}}{b_{k}} \\, .
+        U = \\sum_{k = 0}^{d^{2} - 1} \\, |P_{k})(b_{k}| \\, ,
+
+    where :math:`|P_{k})` is the vectorization of the :math:`k`-th
+    Pauli operator :math:`P_{k}`, and :math:`|k)` is the vectorization
+    of the :math:`k`-th computational basis element.
+    For a definition of vectorization, see :func:`qibo.quantum_info.vectorization`.
 
     Args:
         nqubits (int): number of qubits.
         normalize (bool, optional): If ``True``, converts to the
             Pauli basis. Defaults to ``False``.
         sparse (bool, optional): If ``True``, returns unitary matrix in
-            sparse representation. Default is ``False``.
+            sparse representation. Defaults to ``False``.
         order (str, optional): If ``"row"``, vectorization of Pauli basis is
             performed row-wise. If ``"column"``, vectorization is performed
             column-wise. If ``"system"``, system-wise vectorization is
-            performed. Default is ``"row"``.
+            performed. Defaults to ``"row"``.
         pauli_order (str, optional): corresponds to the order of 4 single-qubit
-            Pauli elements. Default is "IXYZ".
+            Pauli elements. Defaults to ``"IXYZ"``.
         backend (:class:`qibo.backends.abstract.Backend`, optional): backend to be
             used in the execution. If ``None``, it uses
             :class:`qibo.backends.GlobalBackend`. Defaults to ``None``.
diff --git a/src/qibo/quantum_info/clifford.py b/src/qibo/quantum_info/clifford.py
index 000aa6bcd0..f4ec2d6bce 100644
--- a/src/qibo/quantum_info/clifford.py
+++ b/src/qibo/quantum_info/clifford.py
@@ -297,16 +297,16 @@ def frequencies(self, binary: bool = True, registers: bool = False):
             of times each measured value/bitstring appears.
 
             If ``binary`` is ``True``
-                the keys of the `Counter` are in binary form, as strings of
-                :math:`0`s and :math`1`s.
+                the keys of the :class:`collections.Counter` are in binary form,
+                as strings of :math:`0` and :math`1`.
             If ``binary`` is ``False``
-                the keys of the ``Counter`` are integers.
+                the keys of the :class:`collections.Counter` are integers.
             If ``registers`` is ``True``
-                a `dict` of `Counter` s is returned where keys are the name of
-                each register.
+                a `dict` of :class:`collections.Counter` is returned where keys are
+                the name of each register.
             If ``registers`` is ``False``
-                a single ``Counter`` is returned which contains samples from all
-                the measured qubits, independently of their registers.
+                a single :class:`collections.Counter` is returned which contains samples
+                from all the measured qubits, independently of their registers.
         """
         measured_qubits = self.measurement_gate.target_qubits
         freq = self._backend.calculate_frequencies(self.samples(False))
diff --git a/src/qibo/quantum_info/entropies.py b/src/qibo/quantum_info/entropies.py
index 05bf3f823f..c9c889b603 100644
--- a/src/qibo/quantum_info/entropies.py
+++ b/src/qibo/quantum_info/entropies.py
@@ -78,7 +78,7 @@ def classical_relative_entropy(prob_dist_p, prob_dist_q, base: float = 2, backen
 
     For probabilities :math:`\\mathbf{p}` and :math:`\\mathbf{q}`, it is defined as
 
-    ..math::
+    .. math::
         D(\\mathbf{p} \\, \\| \\, \\mathbf{q}) = \\sum_{x} \\, \\mathbf{p}(x) \\,
             \\log\\left( \\frac{\\mathbf{p}(x)}{\\mathbf{q}(x)} \\right) \\, .
 
@@ -680,7 +680,7 @@ def relative_renyi_entropy(
     This is known as the `min-relative entropy <https://arxiv.org/abs/1310.7178>`_.
 
     .. note::
-        Function raises ``NotImplementedError`` when ``target`` :math:`sigma`
+        Function raises ``NotImplementedError`` when ``target`` :math:`\\sigma`
         is a pure state and :math:`\\alpha > 1`. This is due to the fact that
         it is not possible to calculate :math:`\\sigma^{1 - \\alpha}` when
         :math:`\\alpha > 1` and :math:`\\sigma` is a projector, i.e. a singular matrix.
diff --git a/src/qibo/quantum_info/superoperator_transformations.py b/src/qibo/quantum_info/superoperator_transformations.py
index 73bce85f97..158dd73ac4 100644
--- a/src/qibo/quantum_info/superoperator_transformations.py
+++ b/src/qibo/quantum_info/superoperator_transformations.py
@@ -14,18 +14,17 @@
 
 
 def vectorization(state, order: str = "row", backend=None):
-    """Returns state :math:`\\rho` in its Liouville
-    representation :math:`|\\rho\\rangle\\rangle`.
+    """Returns state :math:`\\rho` in its Liouville representation :math:`|\\rho)`.
 
     If ``order="row"``, then:
 
     .. math::
-        |\\rho\\rangle\\rangle = \\sum_{k, l} \\, \\rho_{kl} \\, \\ket{k} \\otimes \\ket{l}
+        |\\rho) = \\sum_{k, l} \\, \\rho_{kl} \\, \\ket{k} \\otimes \\ket{l}
 
     If ``order="column"``, then:
 
     .. math::
-        |\\rho\\rangle\\rangle = \\sum_{k, l} \\, \\rho_{kl} \\, \\ket{l} \\otimes \\ket{k}
+        |\\rho) = \\sum_{k, l} \\, \\rho_{kl} \\, \\ket{l} \\otimes \\ket{k}
 
     Args:
         state: state vector or density matrix.
@@ -88,12 +87,12 @@ def vectorization(state, order: str = "row", backend=None):
 
 def unvectorization(state, order: str = "row", backend=None):
     """Returns state :math:`\\rho` from its Liouville
-    representation :math:`|\\rho\\rangle\\rangle`. This operation is
+    representation :math:`|\\rho)`. This operation is
     the inverse function of :func:`vectorization`, i.e.
 
     .. math::
         \\begin{align}
-            \\rho &= \\text{unvectorization}(|\\rho\\rangle\\rangle) \\nonumber \\\\
+            \\rho &= \\text{unvectorization}(|\\rho)) \\nonumber \\\\
             &= \\text{unvectorization}(\\text{vectorization}(\\rho)) \\nonumber
         \\end{align}
 
@@ -152,9 +151,9 @@ def to_choi(channel, order: str = "row", backend=None):
     """Converts quantum ``channel`` :math:`U` to its Choi representation :math:`\\Lambda`.
 
     .. math::
-        \\Lambda = | U \\rangle\\rangle \\langle\\langle U | \\, ,
+        \\Lambda = | U ) ( U | \\, ,
 
-    where :math:`| \\cdot \\rangle\\rangle` is the :func:`qibo.quantum_info.vectorization`
+    where :math:`| \\cdot )` is the :func:`qibo.quantum_info.vectorization`
     operation.
 
     Args:
@@ -376,7 +375,9 @@ def choi_to_kraus(
 
     .. math::
         \\Lambda = \\sum_{\\alpha} \\, \\lambda_{\\alpha}^{2} \\,
-            |\\tilde{K}_{\\alpha}\\rangle\\rangle \\langle\\langle \\tilde{K}_{\\alpha}| \\, .
+            |\\tilde{K}_{\\alpha})(\\tilde{K}_{\\alpha}| \\, .
+
+    where :math:`|\\cdot)` is the :func:`qibo.quantum_info.vectorization` operation.
 
     This is the spectral decomposition of :math:`\\Lambda`, Hence, the set
     :math:`\\{\\lambda_{\\alpha}, \\, \\tilde{K}_{\\alpha}\\}_{\\alpha}`
@@ -385,7 +386,7 @@ def choi_to_kraus(
 
     .. math::
         K_{\\alpha} = \\lambda_{\\alpha} \\,
-            \\text{unvectorization}(|\\tilde{K}_{\\alpha}\\rangle\\rangle) \\, .
+            \\text{unvectorization}(|\\tilde{K}_{\\alpha})) \\, .
 
     If :math:`\\mathcal{E}` is not CP, then spectral composition is replaced by
     a singular value decomposition (SVD), i.e.
@@ -638,7 +639,11 @@ def kraus_to_choi(kraus_ops, order: str = "row", backend=None):
     of quantum channel to its Choi representation :math:`\\Lambda`.
 
     .. math::
-        \\Lambda = \\sum_{\\alpha} \\, |K_{\\alpha}\\rangle\\rangle \\langle\\langle K_{\\alpha}|
+        \\Lambda = \\sum_{\\alpha} \\, |K_{\\alpha})( K_{\\alpha}|
+
+    where :math:`|K_{\\alpha})` is the vectorization of the Kraus operator
+    :math:`K_{\\alpha}`.
+    For a definition of vectorization, see :func:`qibo.quantum_info.vectorization`.
 
     Args:
         kraus_ops (list): List of Kraus operators as pairs ``(qubits, Ak)``
@@ -757,10 +762,11 @@ def kraus_to_chi(
     of quantum channel to  its :math:`\\chi`-matrix representation.
 
     .. math::
-        \\chi = \\sum_{\\alpha} \\, |c_{\\alpha}\\rangle\\rangle \\langle\\langle c_{\\alpha}|,
+        \\chi = \\sum_{\\alpha} \\, |c_{\\alpha})( c_{\\alpha}|,
 
-    where :math:`|c_{\\alpha}\\rangle\\rangle \\cong |K_{\\alpha}\\rangle\\rangle`
-    in Pauli-Liouville basis.
+    where :math:`|c_{\\alpha}) \\cong |K_{\\alpha})` in Pauli-Liouville basis,
+    and :math:`| \\cdot )` is the :func:`qibo.quantum_info.vectorization`
+    operation.
 
     Args:
         kraus_ops (list): List of Kraus operators as pairs ``(qubits, Ak)``
@@ -2130,8 +2136,6 @@ def _reshuffling(super_op, order: str = "row", backend=None):
     Returns:
         ndarray: Choi (Liouville) representation of the quantum channel.
     """
-    super_op = backend.cast(super_op)
-
     if not isinstance(order, str):
         raise_error(TypeError, f"order must be type str, but it is type {type(order)}.")
 
@@ -2150,6 +2154,8 @@ def _reshuffling(super_op, order: str = "row", backend=None):
 
     backend = _check_backend(backend)
 
+    super_op = backend.cast(super_op, dtype=super_op.dtype)
+
     dim = np.sqrt(super_op.shape[0])
 
     if (
diff --git a/src/qibo/result.py b/src/qibo/result.py
index b2fa8a95fd..e92707c3b2 100644
--- a/src/qibo/result.py
+++ b/src/qibo/result.py
@@ -191,25 +191,27 @@ def frequencies(self, binary: bool = True, registers: bool = False):
         """Returns the frequencies of measured samples.
 
         Args:
-            binary (bool, optional): Return frequency keys in binary or decimal form.
+            binary (bool, optional): If ``True``, returns frequency keys in binary form.
+                If ``False``, returns them in decimal form. Defaults to ``True``.
             registers (bool, optional): Group frequencies according to registers.
+                Defaults to ``False``.
 
         Returns:
-            A `collections.Counter` where the keys are the observed values
+            A :class:`collections.Counter` where the keys are the observed values
             and the values the corresponding frequencies, that is the number
             of times each measured value/bitstring appears.
 
             If ``binary`` is ``True``
-                the keys of the `Counter` are in binary form, as strings of
-                :math:`0`s and :math`1`s.
+                the keys of the :class:`collections.Counter` are in binary form,
+                as strings of :math:`0` and :math`1`.
             If ``binary`` is ``False``
-                the keys of the ``Counter`` are integers.
+                the keys of the :class:`collections.Counter` are integers.
             If ``registers`` is ``True``
-                a `dict` of `Counter` s is returned where keys are the name of
-                each register.
+                a `dict` of :class:`collections.Counter` is returned where keys are
+                the name of each register.
             If ``registers`` is ``False``
-                a single ``Counter`` is returned which contains samples from all
-                the measured qubits, independently of their registers.
+                a single :class:`collections.Counter` is returned which contains samples
+                from all the measured qubits, independently of their registers.
         """
         qubits = self.measurement_gate.qubits
 
diff --git a/src/qibo/transpiler/placer.py b/src/qibo/transpiler/placer.py
index 62d2a3e97c..7944ea81f2 100644
--- a/src/qibo/transpiler/placer.py
+++ b/src/qibo/transpiler/placer.py
@@ -54,7 +54,7 @@ def assert_mapping_consistency(layout: dict, connectivity: nx.Graph = None):
     ref_keys = (
         ["q" + str(i) for i in nodes] if isinstance(physical_qubits[0], str) else nodes
     )
-    if physical_qubits != ref_keys:
+    if sorted(physical_qubits) != sorted(ref_keys):
         raise_error(
             PlacementError,
             "Some physical qubits in the layout may be missing or duplicated.",
@@ -411,7 +411,6 @@ def __init__(
     ):
         self.connectivity = connectivity
         self.routing_algorithm = routing_algorithm
-        self.routing_algorithm.connectivity = connectivity
         self.depth = depth
 
     def __call__(self, circuit: Circuit):
@@ -425,6 +424,7 @@ def __call__(self, circuit: Circuit):
         """
         initial_placer = Trivial(self.connectivity)
         initial_placement = initial_placer(circuit=circuit)
+        self.routing_algorithm.connectivity = self.connectivity
         new_circuit = self._assemble_circuit(circuit)
         final_placement = self._routing_step(initial_placement, new_circuit)
 
diff --git a/src/qibo/transpiler/router.py b/src/qibo/transpiler/router.py
index e3fecb2a4c..b1aad8fc69 100644
--- a/src/qibo/transpiler/router.py
+++ b/src/qibo/transpiler/router.py
@@ -34,13 +34,34 @@ def assert_connectivity(connectivity: nx.Graph, circuit: Circuit):
         if len(gate.qubits) > 2 and not isinstance(gate, gates.M):
             raise_error(ConnectivityError, f"{gate.name} acts on more than two qubits.")
         if len(gate.qubits) == 2:
-            if (gate.qubits[0], gate.qubits[1]) not in connectivity.edges:
+            # physical_qubits = tuple(sorted((circuit.wire_names[gate.qubits[0]], circuit.wire_names[gate.qubits[1]])))
+            physical_qubits = tuple(sorted(gate.qubits))  # for q_i naming
+            if physical_qubits not in connectivity.edges:
                 raise_error(
                     ConnectivityError,
-                    f"Circuit does not respect connectivity. {gate.name} acts on {gate.qubits}.",
+                    f"The circuit does not respect the connectivity. {gate.name} acts on {physical_qubits} but only the following qubits are directly connected: {connectivity.edges}.",
                 )
 
 
+def _relabel_connectivity(connectivity, layout):
+    """Relabels the connectivity graph using the passed layout.
+
+    Args:
+        connectivity (nx.Graph): input connectivity.
+        layout (dict): input qubit layout.
+    Returns:
+        (dict) the updated connectivity.
+    """
+    node_mapping = {}
+    layout = dict(
+        sorted(layout.items(), key=lambda item: int(item[0][1:]))
+    )  # for q_i naming
+    for i, node in enumerate(list(layout.keys())):
+        node_mapping[int(node[1:])] = i  # for q_i naming
+    new_connectivity = nx.relabel_nodes(connectivity, node_mapping)
+    return new_connectivity
+
+
 class StarConnectivityRouter(Router):
     """Transforms an arbitrary circuit to one that can be executed on hardware.
 
@@ -164,7 +185,7 @@ class CircuitMap:
     Also implements the initial two-qubit block decompositions.
 
     Args:
-        initial_layout (dict): initial logical-to-physical qubit mapping.
+        initial_layout (dict): initial physical to logical qubit mapping.
         circuit (:class:`qibo.models.circuit.Circuit`): circuit to be routed.
         blocks (:class:`qibo.transpiler.blocks.CircuitBlocks`, optional): circuit
             block representation. If ``None``, the blocks will be computed from the circuit.
@@ -173,35 +194,86 @@ class CircuitMap:
 
     def __init__(
         self,
-        initial_layout: dict,
-        circuit: Circuit,
+        initial_layout: Optional[dict] = None,
+        circuit: Optional[Circuit] = None,
         blocks: Optional[CircuitBlocks] = None,
+        temp: Optional[bool] = False,
     ):
+        self._p2l, self._l2p = [], []
+
+        self._temporary = temp
+        if self._temporary:
+            return
+        elif circuit is None:
+            raise_error(ValueError, "Circuit must be provided.")
+
         if blocks is not None:
             self.circuit_blocks = blocks
         else:
             self.circuit_blocks = CircuitBlocks(circuit, index_names=True)
-        # Order the initial layout based on the hardware qubit names
-        # to avoid problems in custom layouts
-        self.initial_layout = dict(sorted(initial_layout.items()))
-        self._graph_qubits_names = [int(key[1:]) for key in self.initial_layout.keys()]
-        self._circuit_logical = list(range(len(self.initial_layout)))
-        self._physical_logical = list(self.initial_layout.values())
+
+        self._nqubits = circuit.nqubits
         self._routed_blocks = CircuitBlocks(Circuit(circuit.nqubits))
         self._swaps = 0
 
-    def set_circuit_logical(self, circuit_logical_map: list):
-        """Sets the current circuit to logical qubit mapping.
+        if initial_layout is None:
+            return
 
-        Method works in-place.
+        self.wire_names = list(initial_layout.keys())
+        self.physical_to_logical = list(initial_layout.values())
+
+    @property
+    def physical_to_logical(self):
+        """Returns the physical to logical qubit mapping."""
+        return self._p2l
+
+    @property
+    def logical_to_physical(self):
+        """Returns the logical to physical qubit mapping."""
+        return self._l2p
+
+    @physical_to_logical.setter
+    def physical_to_logical(self, p2l_map: list):
+        """Sets the physical to logical qubit mapping and updates the logical to physical mapping.
 
         Args:
-            circuit_logical_map (list): logical mapping.
+            p2l_map (list): physical to logical mapping.
         """
-        self._circuit_logical = circuit_logical_map
+        self._p2l = p2l_map.copy()
+        self._l2p = [0] * len(self._p2l)
+        for i, l in enumerate(self._p2l):
+            self._l2p[l] = i
+
+    @logical_to_physical.setter
+    def logical_to_physical(self, l2p_map: list):
+        """Sets the logical to physical qubit mapping and updates the physical to logical mapping.
 
-    def blocks_qubits_pairs(self):
-        """Returns a list containing the qubit pairs of each block."""
+        Args:
+            l2p_map (list): logical to physical mapping.
+        """
+        self._l2p = l2p_map.copy()
+        self._p2l = [0] * len(self._l2p)
+        for i, p in enumerate(self._l2p):
+            self._p2l[p] = i
+
+    def _update_mappings_swap(self, logical_swap: tuple, physical_swap: tuple):
+        """Updates the qubit mappings after applying a SWAP gate.
+
+        Args:
+            logical_swap (tuple[int]): the indices of the logical qubits to be swapped.
+            physical_swap (tuple[int]): the indices of the corresponding physical qubits to be swapped.
+        """
+        self._p2l[physical_swap[0]], self._p2l[physical_swap[1]] = (
+            logical_swap[1],
+            logical_swap[0],
+        )
+        self._l2p[logical_swap[0]], self._l2p[logical_swap[1]] = (
+            physical_swap[1],
+            physical_swap[0],
+        )
+
+    def blocks_logical_qubits_pairs(self):
+        """Returns a list containing the logical qubit pairs of each block."""
         return [block.qubits for block in self.circuit_blocks()]
 
     def execute_block(self, block: Block):
@@ -213,9 +285,7 @@ def execute_block(self, block: Block):
         Args:
             block (:class:`qibo.transpiler.blocks.Block`): block to be removed.
         """
-        self._routed_blocks.add_block(
-            block.on_qubits(self.get_physical_qubits(block, index=True))
-        )
+        self._routed_blocks.add_block(block.on_qubits(self.get_physical_qubits(block)))
         self.circuit_blocks.remove_block(block)
 
     def routed_circuit(self, circuit_kwargs: Optional[dict] = None):
@@ -230,16 +300,12 @@ def routed_circuit(self, circuit_kwargs: Optional[dict] = None):
         return self._routed_blocks.circuit(circuit_kwargs=circuit_kwargs)
 
     def final_layout(self):
-        """Returns the final physical-circuit qubits mapping."""
-        unsorted_dict = {
-            "q" + str(self.circuit_to_physical(i)): i
-            for i in range(len(self._circuit_logical))
-        }
+        """Returns the final physical-logical qubits mapping."""
 
-        return dict(sorted(unsorted_dict.items()))
+        return {self.wire_names[i]: self._p2l[i] for i in range(self._nqubits)}
 
-    def update(self, swap: tuple):
-        """Updates the logical-physical qubit mapping after applying a ``SWAP``
+    def update(self, logical_swap: tuple):
+        """Updates the qubit mapping after applying a ``SWAP``
 
         Adds the :class:`qibo.gates.gates.SWAP` gate to the routed blocks.
         Method works in-place.
@@ -247,113 +313,50 @@ def update(self, swap: tuple):
         Args:
             swap (tuple): tuple containing the logical qubits to be swapped.
         """
-        physical_swap = self.logical_to_physical(swap, index=True)
-        self._routed_blocks.add_block(
-            Block(qubits=physical_swap, gates=[gates.SWAP(*physical_swap)])
-        )
-        self._swaps += 1
-        idx_0, idx_1 = self._circuit_logical.index(
-            swap[0]
-        ), self._circuit_logical.index(swap[1])
-        self._circuit_logical[idx_0], self._circuit_logical[idx_1] = swap[1], swap[0]
+
+        physical_swap = self.logical_pair_to_physical(logical_swap)
+        if not self._temporary:
+            self._routed_blocks.add_block(
+                Block(qubits=physical_swap, gates=[gates.SWAP(*physical_swap)])
+            )
+            self._swaps += 1
+
+        self._update_mappings_swap(logical_swap, physical_swap)
 
     def undo(self):
         """Undo the last swap. Method works in-place."""
         last_swap_block = self._routed_blocks.return_last_block()
-        swap = tuple(self.physical_to_logical(q) for q in last_swap_block.qubits)
+        physical_swap = last_swap_block.qubits
+        logical_swap = self._p2l[physical_swap[0]], self._p2l[physical_swap[1]]
         self._routed_blocks.remove_block(last_swap_block)
         self._swaps -= 1
 
-        idx_0, idx_1 = self._circuit_logical.index(
-            swap[0]
-        ), self._circuit_logical.index(swap[1])
-        self._circuit_logical[idx_0], self._circuit_logical[idx_1] = swap[1], swap[0]
-
-    def get_logical_qubits(self, block: Block):
-        """Returns the current logical qubits where a block is acting on.
-
-        Args:
-            block (:class:`qibo.transpiler.blocks.Block`): block to be analysed.
+        self._update_mappings_swap(logical_swap, physical_swap)
 
-        Returns:
-            tuple: logical qubits where a block is acting on.
-        """
-        return self.circuit_to_logical(block.qubits)
-
-    def get_physical_qubits(self, block: Union[int, Block], index: bool = False):
+    def get_physical_qubits(self, block: Union[int, Block]):
         """Returns the physical qubits where a block is acting on.
 
         Args:
             block (int or :class:`qibo.transpiler.blocks.Block`): block to be analysed.
-            index (bool, optional): If ``True``, qubits are returned as indices of
-                the connectivity nodes. Defaults to ``False``.
 
         Returns:
-            tuple: physical qubits where a block is acting on.
-
+            tuple: physical qubit numbers where a block is acting on.
         """
         if isinstance(block, int):
             block = self.circuit_blocks.search_by_index(block)
 
-        return self.logical_to_physical(self.get_logical_qubits(block), index=index)
-
-    def logical_to_physical(self, logical_qubits: tuple, index: bool = False):
-        """Returns the physical qubits associated to the logical qubits.
-
-        Args:
-            logical_qubits (tuple): physical qubits.
-            index (bool, optional): If ``True``, qubits are returned as indices of
-                `the connectivity nodes. Defaults to ``False``.
-
-        Returns:
-            tuple: physical qubits associated to the logical qubits.
-        """
-        if not index:
-            return tuple(
-                self._graph_qubits_names[
-                    self._physical_logical.index(logical_qubits[i])
-                ]
-                for i in range(2)
-            )
-
-        return tuple(self._physical_logical.index(logical_qubits[i]) for i in range(2))
-
-    def circuit_to_logical(self, circuit_qubits: tuple):
-        """Returns the current logical qubits associated to the initial circuit qubits.
-
-        Args:
-            circuit_qubits (tuple): circuit qubits.
-
-        Returns:
-            tuple: logical qubits.
-        """
-        return tuple(self._circuit_logical[circuit_qubits[i]] for i in range(2))
-
-    def circuit_to_physical(self, circuit_qubit: int):
-        """Returns the current physical qubit associated to an initial circuit qubit.
-
-        Args:
-            circuit_qubit (int): circuit qubit.
-
-        Returns:
-            int: physical qubit.
-        """
-        return self._graph_qubits_names[
-            self._physical_logical.index(self._circuit_logical[circuit_qubit])
-        ]
+        return tuple(self._l2p[q] for q in block.qubits)
 
-    def physical_to_logical(self, physical_qubit: int):
-        """Returns current logical qubit associated to a physical qubit (connectivity graph node).
+    def logical_pair_to_physical(self, logical_qubits: tuple):
+        """Returns the physical qubits associated to the logical qubit pair.
 
         Args:
-            physical_qubit (int): physical qubit.
+            logical_qubits (tuple): logical qubit pair.
 
         Returns:
-            int: logical qubit.
+            tuple: physical qubit numbers associated to the logical qubit pair.
         """
-        physical_qubit_index = self._graph_qubits_names.index(physical_qubit)
-
-        return self._physical_logical[physical_qubit_index]
+        return self._l2p[logical_qubits[0]], self._l2p[logical_qubits[1]]
 
 
 class ShortestPaths(Router):
@@ -368,9 +371,10 @@ class ShortestPaths(Router):
     def __init__(self, connectivity: nx.Graph, seed: Optional[int] = None):
         self.connectivity = connectivity
         self._front_layer = None
-        self.circuit = None
+        self.circuit_map = None
         self._dag = None
         self._final_measurements = None
+        self._node_mapping_inv = None
         if seed is None:
             seed = 42
         random.seed(seed)
@@ -378,7 +382,7 @@ def __init__(self, connectivity: nx.Graph, seed: Optional[int] = None):
     @property
     def added_swaps(self):
         """Returns the number of SWAP gates added to the circuit during routing."""
-        return self.circuit._swaps
+        return self.circuit_map._swaps
 
     def __call__(self, circuit: Circuit, initial_layout: dict):
         """Circuit connectivity matching.
@@ -402,13 +406,13 @@ def __call__(self, circuit: Circuit, initial_layout: dict):
 
         circuit_kwargs = circuit.init_kwargs
         circuit_kwargs["wire_names"] = list(initial_layout.keys())
-        routed_circuit = self.circuit.routed_circuit(circuit_kwargs=circuit_kwargs)
+        routed_circuit = self.circuit_map.routed_circuit(circuit_kwargs=circuit_kwargs)
         if self._final_measurements is not None:
             routed_circuit = self._append_final_measurements(
                 routed_circuit=routed_circuit
             )
 
-        return routed_circuit, self.circuit.final_layout()
+        return routed_circuit, self.circuit_map.final_layout()
 
     def _find_new_mapping(self):
         """Find new qubit mapping. Mapping is found by looking for the shortest path.
@@ -426,13 +430,13 @@ def _find_new_mapping(self):
             if candidate[1] == best_cost
         ]
         best_candidate = random.choice(best_candidates)
-        self._add_swaps(best_candidate, self.circuit)
+        self._add_swaps(best_candidate, self.circuit_map)
 
     def _candidates(self):
         """Returns all possible shortest paths in a ``list`` that contains
         the new mapping and a second ``list`` containing the path meeting point.
         """
-        target_qubits = self.circuit.get_physical_qubits(self._front_layer[0])
+        target_qubits = self.circuit_map.get_physical_qubits(self._front_layer[0])
         path_list = list(
             nx.all_shortest_paths(
                 self.connectivity, source=target_qubits[0], target=target_qubits[1]
@@ -459,22 +463,20 @@ def _add_swaps(candidate: tuple, circuitmap: CircuitMap):
         meeting_point = candidate[1]
         forward = path[0 : meeting_point + 1]
         backward = list(reversed(path[meeting_point + 1 :]))
-        if len(forward) > 1:
-            for f1, f2 in zip(forward[:-1], forward[1:]):
-                circuitmap.update(
-                    (
-                        circuitmap.physical_to_logical(f1),
-                        circuitmap.physical_to_logical(f2),
-                    )
+        for f in forward[1:]:
+            circuitmap.update(
+                (
+                    circuitmap.physical_to_logical[f],
+                    circuitmap.physical_to_logical[forward[0]],
                 )
-        if len(backward) > 1:
-            for b1, b2 in zip(backward[:-1], backward[1:]):
-                circuitmap.update(
-                    (
-                        circuitmap.physical_to_logical(b1),
-                        circuitmap.physical_to_logical(b2),
-                    )
+            )
+        for b in backward[1:]:
+            circuitmap.update(
+                (
+                    circuitmap.physical_to_logical[b],
+                    circuitmap.physical_to_logical[backward[0]],
                 )
+            )
 
     def _compute_cost(self, candidate: tuple):
         """Greedy algorithm that decides which path to take and how qubits should be walked.
@@ -488,11 +490,11 @@ def _compute_cost(self, candidate: tuple):
             (list, int): best path to move qubits and qubit meeting point in the path.
         """
         temporary_circuit = CircuitMap(
-            initial_layout=self.circuit.initial_layout,
-            circuit=Circuit(len(self.circuit.initial_layout)),
-            blocks=deepcopy(self.circuit.circuit_blocks),
+            circuit=Circuit(self.circuit_map._nqubits),
+            blocks=deepcopy(self.circuit_map.circuit_blocks),
         )
-        temporary_circuit.set_circuit_logical(deepcopy(self.circuit._circuit_logical))
+
+        temporary_circuit.physical_to_logical = self.circuit_map.physical_to_logical
         self._add_swaps(candidate, temporary_circuit)
         temporary_dag = deepcopy(self._dag)
         successive_executed_gates = 0
@@ -525,7 +527,7 @@ def _compute_cost(self, candidate: tuple):
         return -successive_executed_gates
 
     def _check_execution(self):
-        """Checks if some blocks in the front layer can be executed in the current configuration.
+        """Check if some blocks in the front layer can be executed in the current configuration.
 
         Returns:
             (list): executable blocks if there are, ``None`` otherwise.
@@ -533,8 +535,8 @@ def _check_execution(self):
         executable_blocks = []
         for block in self._front_layer:
             if (
-                self.circuit.get_physical_qubits(block) in self.connectivity.edges
-                or not self.circuit.circuit_blocks.search_by_index(block).entangled
+                self.circuit_map.get_physical_qubits(block) in self.connectivity.edges
+                or not self.circuit_map.circuit_blocks.search_by_index(block).entangled
             ):
                 executable_blocks.append(block)
         if len(executable_blocks) == 0:
@@ -554,8 +556,8 @@ def _execute_blocks(self, blocklist: list):
             blocklist (list): list of blocks.
         """
         for block_id in blocklist:
-            block = self.circuit.circuit_blocks.search_by_index(block_id)
-            self.circuit.execute_block(block)
+            block = self.circuit_map.circuit_blocks.search_by_index(block_id)
+            self.circuit_map.execute_block(block)
             self._dag.remove_node(block_id)
         self._update_front_layer()
 
@@ -581,10 +583,13 @@ def _preprocessing(self, circuit: Circuit, initial_layout: dict):
             circuit (:class:`qibo.models.circuit.Circuit`): circuit to be preprocessed.
             initial_layout (dict): initial physical-to-logical qubit mapping.
         """
+
+        self.connectivity = _relabel_connectivity(self.connectivity, initial_layout)
+
         copied_circuit = circuit.copy(deep=True)
         self._final_measurements = self._detach_final_measurements(copied_circuit)
-        self.circuit = CircuitMap(initial_layout, copied_circuit)
-        self._dag = _create_dag(self.circuit.blocks_qubits_pairs())
+        self.circuit_map = CircuitMap(initial_layout, copied_circuit)
+        self._dag = _create_dag(self.circuit_map.blocks_logical_qubits_pairs())
         self._update_front_layer()
 
     def _detach_final_measurements(self, circuit: Circuit):
@@ -613,8 +618,8 @@ def _append_final_measurements(self, routed_circuit: Circuit):
         conserving the measurement register."""
         for measurement in self._final_measurements:
             original_qubits = measurement.qubits
-            routed_qubits = (
-                self.circuit.circuit_to_physical(qubit) for qubit in original_qubits
+            routed_qubits = list(
+                self.circuit_map.logical_to_physical[qubit] for qubit in original_qubits
             )
             routed_circuit.add(
                 measurement.on_qubits(dict(zip(original_qubits, routed_qubits)))
@@ -668,7 +673,7 @@ def __init__(
         self._dist_matrix = None
         self._dag = None
         self._front_layer = None
-        self.circuit = None
+        self.circuit_map = None
         self._memory_map = None
         self._final_measurements = None
         self._temp_added_swaps = []
@@ -700,25 +705,25 @@ def __call__(self, circuit: Circuit, initial_layout: dict):
                 len(self._temp_added_swaps) > self.swap_threshold * longest_path
             ):  # threshold is arbitrary
                 while self._temp_added_swaps:
-                    swap = self._temp_added_swaps.pop()
-                    self.circuit.undo()
+                    self._temp_added_swaps.pop()
+                    self.circuit_map.undo()
                 self._temp_added_swaps = []
                 self._shortest_path_routing()
 
         circuit_kwargs = circuit.init_kwargs
         circuit_kwargs["wire_names"] = list(initial_layout.keys())
-        routed_circuit = self.circuit.routed_circuit(circuit_kwargs=circuit_kwargs)
+        routed_circuit = self.circuit_map.routed_circuit(circuit_kwargs=circuit_kwargs)
         if self._final_measurements is not None:
             routed_circuit = self._append_final_measurements(
                 routed_circuit=routed_circuit
             )
 
-        return routed_circuit, self.circuit.final_layout()
+        return routed_circuit, self.circuit_map.final_layout()
 
     @property
     def added_swaps(self):
         """Returns the number of SWAP gates added to the circuit during routing."""
-        return self.circuit._swaps
+        return self.circuit_map._swaps
 
     def _preprocessing(self, circuit: Circuit, initial_layout: dict):
         """The following objects will be initialised:
@@ -735,11 +740,14 @@ def _preprocessing(self, circuit: Circuit, initial_layout: dict):
             circuit (:class:`qibo.models.circuit.Circuit`): circuit to be preprocessed.
             initial_layout (dict): initial physical-to-logical qubit mapping.
         """
+
+        self.connectivity = _relabel_connectivity(self.connectivity, initial_layout)
+
         copied_circuit = circuit.copy(deep=True)
         self._final_measurements = self._detach_final_measurements(copied_circuit)
-        self.circuit = CircuitMap(initial_layout, copied_circuit)
+        self.circuit_map = CircuitMap(initial_layout, copied_circuit)
         self._dist_matrix = nx.floyd_warshall_numpy(self.connectivity)
-        self._dag = _create_dag(self.circuit.blocks_qubits_pairs())
+        self._dag = _create_dag(self.circuit_map.blocks_logical_qubits_pairs())
         self._memory_map = []
         self._update_dag_layers()
         self._update_front_layer()
@@ -770,7 +778,7 @@ def _append_final_measurements(self, routed_circuit: Circuit):
         for measurement in self._final_measurements:
             original_qubits = measurement.qubits
             routed_qubits = list(
-                self.circuit.circuit_to_physical(qubit) for qubit in original_qubits
+                self.circuit_map.logical_to_physical[qubit] for qubit in original_qubits
             )
             routed_circuit.add(
                 measurement.on_qubits(dict(zip(original_qubits, routed_qubits)))
@@ -793,14 +801,32 @@ def _update_front_layer(self):
         """
         self._front_layer = self._get_dag_layer(0)
 
-    def _get_dag_layer(self, n_layer):
-        """Return the :math:`n`-topological layer of the dag."""
+    def _get_dag_layer(self, n_layer, qubits=False):
+        """Return the :math:`n`-topological layer of the dag.
+
+        Args:
+            n_layer (int): layer number.
+            qubits (bool, optional): if ``True``, return the target qubits of the blocks in the layer.
+                If ``False``, return the block numbers. Defaults to ``False``.
+
+        Returns:
+            (list): list of block numbers or target qubits.
+        """
+
+        if qubits:
+            return [
+                node[1]["qubits"]
+                for node in self._dag.nodes(data=True)
+                if node[1]["layer"] == n_layer
+            ]
+
         return [node[0] for node in self._dag.nodes(data="layer") if node[1] == n_layer]
 
     def _find_new_mapping(self):
         """Find the new best mapping by adding one swap."""
         candidates_evaluation = {}
-        self._memory_map.append(deepcopy(self.circuit._circuit_logical))
+
+        self._memory_map.append(self.circuit_map.physical_to_logical.copy())
         for candidate in self._swap_candidates():
             candidates_evaluation[candidate] = self._compute_cost(candidate)
 
@@ -810,31 +836,28 @@ def _find_new_mapping(self):
         ]
         best_candidate = random.choice(best_candidates)
 
-        for qubit in self.circuit.logical_to_physical(best_candidate, index=True):
+        for qubit in self.circuit_map.logical_pair_to_physical(best_candidate):
             self._delta_register[qubit] += self.delta
-        self.circuit.update(best_candidate)
+        self.circuit_map.update(best_candidate)
         self._temp_added_swaps.append(best_candidate)
 
     def _compute_cost(self, candidate: int):
         """Compute the cost associated to a possible SWAP candidate."""
-        temporary_circuit = CircuitMap(
-            initial_layout=self.circuit.initial_layout,
-            circuit=Circuit(len(self.circuit.initial_layout)),
-            blocks=self.circuit.circuit_blocks,
-        )
-        temporary_circuit.set_circuit_logical(deepcopy(self.circuit._circuit_logical))
+
+        temporary_circuit = CircuitMap(temp=True)
+        temporary_circuit.physical_to_logical = self.circuit_map.physical_to_logical
         temporary_circuit.update(candidate)
 
-        if temporary_circuit._circuit_logical in self._memory_map:
+        if temporary_circuit.physical_to_logical in self._memory_map:
             return float("inf")
 
         tot_distance = 0.0
         weight = 1.0
         for layer in range(self.lookahead + 1):
-            layer_gates = self._get_dag_layer(layer)
+            layer_gates = self._get_dag_layer(layer, qubits=True)
             avg_layer_distance = 0.0
-            for gate in layer_gates:
-                qubits = temporary_circuit.get_physical_qubits(gate, index=True)
+            for lq_pair in layer_gates:
+                qubits = temporary_circuit.logical_pair_to_physical(lq_pair)
                 avg_layer_distance += (
                     max(self._delta_register[i] for i in qubits)
                     * (self._dist_matrix[qubits[0], qubits[1]] - 1.0)
@@ -856,13 +879,13 @@ def _swap_candidates(self):
         """
         candidates = []
         for block in self._front_layer:
-            for qubit in self.circuit.get_physical_qubits(block):
+            for qubit in self.circuit_map.get_physical_qubits(block):
                 for connected in self.connectivity.neighbors(qubit):
                     candidate = tuple(
                         sorted(
                             (
-                                self.circuit.physical_to_logical(qubit),
-                                self.circuit.physical_to_logical(connected),
+                                self.circuit_map.physical_to_logical[qubit],
+                                self.circuit_map.physical_to_logical[connected],
                             )
                         )
                     )
@@ -880,8 +903,8 @@ def _check_execution(self):
         executable_blocks = []
         for block in self._front_layer:
             if (
-                self.circuit.get_physical_qubits(block) in self.connectivity.edges
-                or not self.circuit.circuit_blocks.search_by_index(block).entangled
+                self.circuit_map.get_physical_qubits(block) in self.connectivity.edges
+                or not self.circuit_map.circuit_blocks.search_by_index(block).entangled
             ):
                 executable_blocks.append(block)
 
@@ -903,8 +926,8 @@ def _execute_blocks(self, blocklist: list):
             blocklist (list): list of blocks.
         """
         for block_id in blocklist:
-            block = self.circuit.circuit_blocks.search_by_index(block_id)
-            self.circuit.execute_block(block)
+            block = self.circuit_map.circuit_blocks.search_by_index(block_id)
+            self.circuit_map.execute_block(block)
             self._dag.remove_node(block_id)
         self._update_dag_layers()
         self._update_front_layer()
@@ -922,7 +945,7 @@ def _shortest_path_routing(self):
         shortest_path_qubits = None
 
         for block in self._front_layer:
-            q1, q2 = self.circuit.get_physical_qubits(block)
+            q1, q2 = self.circuit_map.get_physical_qubits(block)
             distance = self._dist_matrix[q1, q2]
 
             if distance < min_distance:
@@ -933,14 +956,10 @@ def _shortest_path_routing(self):
             self.connectivity, shortest_path_qubits[0], shortest_path_qubits[1]
         )
 
-        # Q1 is moved
-        shortest_path = [
-            self.circuit.physical_to_logical(q) for q in shortest_path[:-1]
-        ]
-        swaps = list(zip(shortest_path[:-1], shortest_path[1:]))
-
-        for swap in swaps:
-            self.circuit.update(swap)
+        # move q1
+        q1 = self.circuit_map.physical_to_logical[shortest_path[0]]
+        for q2 in shortest_path[1:-1]:
+            self.circuit_map.update((q1, self.circuit_map.physical_to_logical[q2]))
 
 
 def _create_dag(gates_qubits_pairs: list):
@@ -957,6 +976,10 @@ def _create_dag(gates_qubits_pairs: list):
     """
     dag = nx.DiGraph()
     dag.add_nodes_from(range(len(gates_qubits_pairs)))
+
+    for i in range(len(gates_qubits_pairs)):
+        dag.nodes[i]["qubits"] = gates_qubits_pairs[i]
+
     # Find all successors
     connectivity_list = []
     for idx, gate in enumerate(gates_qubits_pairs):
@@ -985,7 +1008,7 @@ def _remove_redundant_connections(dag: nx.DiGraph):
         (:class:`networkx.DiGraph`): reduced dag.
     """
     new_dag = nx.DiGraph()
-    new_dag.add_nodes_from(range(dag.number_of_nodes()))
+    new_dag.add_nodes_from(dag.nodes(data=True))
     transitive_reduction = nx.transitive_reduction(dag)
     new_dag.add_edges_from(transitive_reduction.edges)
 
diff --git a/src/qibo/ui/mpldrawer.py b/src/qibo/ui/mpldrawer.py
index 778a429fe2..439541d38d 100644
--- a/src/qibo/ui/mpldrawer.py
+++ b/src/qibo/ui/mpldrawer.py
@@ -227,20 +227,6 @@ def _draw_controls(ax, i, gate, labels, gate_grid, wire_grid, plot_params, measu
             wire_grid[max_wire],
             plot_params,
         )
-        ismeasured = False
-        for index in control_indices:
-            if measured.get(index, 1000) < i:
-                ismeasured = True
-        if ismeasured:
-            dy = 0.04  # TODO: put in plot_params
-            _line(
-                ax,
-                gate_grid[i] + dy,
-                gate_grid[i] + dy,
-                wire_grid[min_wire],
-                wire_grid[max_wire],
-                plot_params,
-            )
 
         for ci in control_indices:
             x = gate_grid[i]
@@ -301,12 +287,8 @@ def _draw_target(ax, i, gate, labels, gate_grid, wire_grid, plot_params):
     x = gate_grid[i]
     target_index = _get_flipped_index(target, labels)
     y = wire_grid[target_index]
-    if not symbol:
-        return
     if name in ["CNOT", "TOFFOLI"]:
         _oplus(ax, x, y, plot_params)
-    elif name == "CPHASE":
-        _cdot(ax, x, y, plot_params)
     elif name == "SWAP":
         _swapx(ax, x, y, plot_params)
     else:
@@ -553,12 +535,13 @@ def _make_cluster_gates(gates_items):
     return cluster_gates
 
 
-def _process_gates(array_gates):
+def _process_gates(array_gates, nqubits):
     """
     Transforms the list of gates given by the Qibo circuit into a list of gates with a suitable structre to print on screen with matplotlib.
 
     Args:
         array_gates (list): List of gates provided by the Qibo circuit.
+        nqubits (int): Number of circuit qubits
 
     Returns:
         list: List of suitable gates to plot with matplotlib.
@@ -601,7 +584,7 @@ def _process_gates(array_gates):
                 item += ("q_" + str(qbit),)
                 gates_plot.append(item)
         elif init_label == "ENTANGLEMENTENTROPY":
-            for qbit in list(range(circuit.nqubits)):
+            for qbit in list(range(nqubits)):
                 item = (init_label,)
                 item += ("q_" + str(qbit),)
                 gates_plot.append(item)
@@ -610,16 +593,18 @@ def _process_gates(array_gates):
             item += (init_label,)
 
             for qbit in gate._target_qubits:
-                if qbit is tuple:
+                if type(qbit) is tuple:
                     item += ("q_" + str(qbit[0]),)
                 else:
                     item += ("q_" + str(qbit),)
 
             for qbit in gate._control_qubits:
-                if qbit is tuple:
-                    item += ("q_" + str(qbit[0]),)
-                else:
-                    item += ("q_" + str(qbit),)
+                item_add = (
+                    ("q_" + str(qbit[0]),)
+                    if isinstance(qbit, tuple)
+                    else ("q_" + str(qbit),)
+                )
+                item += item_add
 
             gates_plot.append(item)
 
@@ -637,10 +622,11 @@ def _plot_params(style: Union[dict, str, None]) -> dict:
         dict: Style configuration.
     """
     if not isinstance(style, dict):
-        try:
-            style = STYLE.get(style) if (style is not None) else STYLE["default"]
-        except AttributeError:
-            style = STYLE["default"]
+        style = (
+            STYLE.get(style)
+            if (style is not None and style in STYLE.keys())
+            else STYLE["default"]
+        )
 
     return style
 
@@ -718,9 +704,11 @@ def plot_circuit(circuit, scale=0.6, cluster_gates=True, style=None):
             fgates = None
 
             if cluster_gates:
-                fgates = _make_cluster_gates(_process_gates(gate.gates))
+                fgates = _make_cluster_gates(
+                    _process_gates(gate.gates, circuit.nqubits)
+                )
             else:
-                fgates = _process_gates(gate.gates)
+                fgates = _process_gates(gate.gates, circuit.nqubits)
 
             l_gates = len(gate.gates)
             equal_qbits = False
@@ -736,7 +724,7 @@ def plot_circuit(circuit, scale=0.6, cluster_gates=True, style=None):
         else:
             all_gates.append(gate)
 
-    gates_plot = _process_gates(all_gates)
+    gates_plot = _process_gates(all_gates, circuit.nqubits)
 
     if cluster_gates and len(gates_plot) > 0:
         gates_cluster = _make_cluster_gates(gates_plot)
diff --git a/tests/test_transpiler_router.py b/tests/test_transpiler_router.py
index 806223cbc8..d5901d5f92 100644
--- a/tests/test_transpiler_router.py
+++ b/tests/test_transpiler_router.py
@@ -50,6 +50,15 @@ def grid_connectivity():
     return chip
 
 
+def line_connectivity(n):
+    Q = [i for i in range(n)]
+    chip = nx.Graph()
+    chip.add_nodes_from(Q)
+    graph_list = [(Q[i], (Q[i] + 1) % n) for i in range(n - 1)]
+    chip.add_edges_from(graph_list)
+    return chip
+
+
 def generate_random_circuit(nqubits, ngates, seed=42):
     """Generate a random circuit with RX and CZ gates."""
     np.random.seed(seed)
@@ -136,6 +145,29 @@ def test_random_circuits_5q(gates, placer, connectivity):
     )
 
 
+def test_random_circuits_15q_50g():
+    nqubits, ngates = 15, 50
+    connectivity = line_connectivity(nqubits)
+    placer = Random(connectivity=connectivity)
+    layout_circ = Circuit(nqubits)
+    initial_layout = placer(layout_circ)
+    transpiler = Sabre(connectivity=connectivity)
+    circuit = generate_random_circuit(nqubits=nqubits, ngates=ngates)
+    transpiled_circuit, final_qubit_map = transpiler(circuit, initial_layout)
+    assert transpiler.added_swaps >= 0
+    assert_connectivity(connectivity, transpiled_circuit)
+    assert_placement(transpiled_circuit, final_qubit_map)
+    assert ngates + transpiler.added_swaps == transpiled_circuit.ngates
+    qubit_matcher = Preprocessing(connectivity=connectivity)
+    new_circuit = qubit_matcher(circuit=circuit)
+    assert_circuit_equivalence(
+        original_circuit=new_circuit,
+        transpiled_circuit=transpiled_circuit,
+        final_map=final_qubit_map,
+        initial_map=initial_layout,
+    )
+
+
 def test_star_circuit():
     placer = Subgraph(star_connectivity())
     initial_layout = placer(star_circuit())
@@ -254,8 +286,8 @@ def test_sabre_shortest_path_routing():
     router._preprocessing(circuit=loop_circ, initial_layout=initial_layout)
     router._shortest_path_routing()  # q2 should be moved adjacent to q8
 
-    gate_28 = router.circuit.circuit_blocks.block_list[2]
-    gate_28_qubits = router.circuit.get_physical_qubits(gate_28)
+    gate_28 = router.circuit_map.circuit_blocks.block_list[2]
+    gate_28_qubits = router.circuit_map.get_physical_qubits(gate_28)
 
     # Check if the physical qubits of the gate (2, 8) are adjacent
     assert gate_28_qubits[1] in list(router.connectivity.neighbors(gate_28_qubits[0]))
@@ -275,37 +307,69 @@ def test_circuit_map():
     circuit_map = CircuitMap(initial_layout=initial_layout, circuit=circ)
     block_list = circuit_map.circuit_blocks
     # test blocks_qubits_pairs
-    assert circuit_map.blocks_qubits_pairs() == [(0, 1), (1, 2), (0, 1), (2, 3)]
+    assert circuit_map.blocks_logical_qubits_pairs() == [(0, 1), (1, 2), (0, 1), (2, 3)]
     # test execute_block and routed_circuit
     circuit_map.execute_block(block_list.search_by_index(0))
     routed_circuit = circuit_map.routed_circuit()
     assert isinstance(routed_circuit.queue[0], gates.H)
     assert len(routed_circuit.queue) == 4
-    assert routed_circuit.queue[2].qubits == (1, 2)
-    # test update
+    qubits = routed_circuit.queue[2].qubits
+    assert (
+        routed_circuit.wire_names[qubits[0]] == "q1"
+        and routed_circuit.wire_names[qubits[1]] == "q2"
+    )
+
+    # test update 1
     circuit_map.update((0, 2))
     routed_circuit = circuit_map.routed_circuit()
     assert isinstance(routed_circuit.queue[4], gates.SWAP)
-    assert routed_circuit.queue[4].qubits == (1, 0)
+    qubits = routed_circuit.queue[4].qubits
+    assert (
+        routed_circuit.wire_names[qubits[0]] == "q1"
+        and routed_circuit.wire_names[qubits[1]] == "q0"
+    )
     assert circuit_map._swaps == 1
-    assert circuit_map._circuit_logical == [2, 1, 0, 3]
+    assert circuit_map.physical_to_logical == [0, 2, 1, 3]
+    assert circuit_map.logical_to_physical == [0, 2, 1, 3]
+
+    # test update 2
     circuit_map.update((1, 2))
     routed_circuit = circuit_map.routed_circuit()
-    assert routed_circuit.queue[5].qubits == (2, 0)
-    assert circuit_map._circuit_logical == [1, 2, 0, 3]
-    # test execute_block after multiple swaps
+    assert isinstance(routed_circuit.queue[5], gates.SWAP)
+    qubits = routed_circuit.queue[5].qubits
+    assert (
+        routed_circuit.wire_names[qubits[0]] == "q2"
+        and routed_circuit.wire_names[qubits[1]] == "q1"
+    )
+    assert circuit_map._swaps == 2
+    assert circuit_map.physical_to_logical == [0, 1, 2, 3]
+    assert circuit_map.logical_to_physical == [0, 1, 2, 3]
+
+    # # test execute_block after multiple swaps
     circuit_map.execute_block(block_list.search_by_index(1))
     circuit_map.execute_block(block_list.search_by_index(2))
     circuit_map.execute_block(block_list.search_by_index(3))
     routed_circuit = circuit_map.routed_circuit()
     assert isinstance(routed_circuit.queue[6], gates.CZ)
-    # circuit to logical map: [1,2,0,3]. initial map: {"q0": 2, "q1": 0, "q2": 1, "q3": 3}.
-    assert routed_circuit.queue[6].qubits == (0, 1)  # initial circuit qubits (1,2)
-    assert routed_circuit.queue[7].qubits == (2, 0)  # (0,1)
-    assert routed_circuit.queue[8].qubits == (1, 3)  # (2,3)
+
+    qubits = routed_circuit.queue[6].qubits
+    assert (
+        routed_circuit.wire_names[qubits[0]] == "q1"
+        and routed_circuit.wire_names[qubits[1]] == "q2"
+    )
+    qubits = routed_circuit.queue[7].qubits
+    assert (
+        routed_circuit.wire_names[qubits[0]] == "q0"
+        and routed_circuit.wire_names[qubits[1]] == "q1"
+    )
+    qubits = routed_circuit.queue[8].qubits
+    assert (
+        routed_circuit.wire_names[qubits[0]] == "q2"
+        and routed_circuit.wire_names[qubits[1]] == "q3"
+    )
     assert len(circuit_map.circuit_blocks()) == 0
     # test final layout
-    assert circuit_map.final_layout() == {"q0": 1, "q1": 2, "q2": 0, "q3": 3}
+    assert circuit_map.final_layout() == {"q0": 0, "q1": 1, "q2": 2, "q3": 3}
 
 
 def test_sabre_matched():
@@ -483,17 +547,36 @@ def test_undo():
     # Two SWAP gates are added
     circuit_map.update((1, 2))
     circuit_map.update((2, 3))
-    assert circuit_map._circuit_logical == [0, 3, 1, 2]
+    assert circuit_map.physical_to_logical == [0, 3, 1, 2]
+    assert circuit_map.logical_to_physical == [0, 2, 3, 1]
     assert len(circuit_map._routed_blocks.block_list) == 2
+    assert circuit_map._swaps == 2
 
     # Undo the last SWAP gate
     circuit_map.undo()
-    assert circuit_map._circuit_logical == [0, 2, 1, 3]
+    assert circuit_map.physical_to_logical == [0, 2, 1, 3]
+    assert circuit_map.logical_to_physical == [0, 2, 1, 3]
     assert circuit_map._swaps == 1
     assert len(circuit_map._routed_blocks.block_list) == 1
 
     # Undo the first SWAP gate
     circuit_map.undo()
-    assert circuit_map._circuit_logical == [0, 1, 2, 3]
+    assert circuit_map.physical_to_logical == [0, 1, 2, 3]
+    assert circuit_map.logical_to_physical == [0, 1, 2, 3]
     assert circuit_map._swaps == 0
     assert len(circuit_map._routed_blocks.block_list) == 0
+
+
+def test_circuitmap_no_circuit():
+    # If a `CircuitMap` is not a temporary instance and is created without a circuit, it should raise an error.
+    with pytest.raises(ValueError):
+        circuit_map = CircuitMap()
+
+
+def test_logical_to_physical_setter():
+    circ = Circuit(4)
+    initial_layout = {"q0": 0, "q1": 3, "q2": 2, "q3": 1}
+    circuit_map = CircuitMap(initial_layout=initial_layout, circuit=circ)
+    circuit_map.logical_to_physical = [2, 0, 1, 3]
+    assert circuit_map.logical_to_physical == [2, 0, 1, 3]
+    assert circuit_map.physical_to_logical == [1, 2, 0, 3]