diff --git a/mpqp/core/circuit.py b/mpqp/core/circuit.py
index a681b988..b875192f 100644
--- a/mpqp/core/circuit.py
+++ b/mpqp/core/circuit.py
@@ -383,8 +383,6 @@ def _update_targets_components(self, component: Instruction | NoiseModel):
         if isinstance(component, Barrier):
             component.size = self.nb_qubits
             component.targets = list(range(self.nb_qubits))
-        elif isinstance(component, ExpectationMeasure):
-            component._check_targets_order()  # pyright: ignore[reportPrivateUsage]
         elif isinstance(component, DimensionalNoiseModel):
             component.check_dimension()
         elif isinstance(component, BasisMeasure):
diff --git a/mpqp/core/instruction/measurement/expectation_value.py b/mpqp/core/instruction/measurement/expectation_value.py
index d09a9c0d..39cd5e21 100644
--- a/mpqp/core/instruction/measurement/expectation_value.py
+++ b/mpqp/core/instruction/measurement/expectation_value.py
@@ -8,13 +8,11 @@ class to define your observable, and a :class:`ExpectationMeasure` to perform
 import copy
 from numbers import Real
 from typing import TYPE_CHECKING, Literal, Optional, Union, overload
-from typing_extensions import Never
-from warnings import warn
 
 import numpy as np
 import numpy.typing as npt
+from typing_extensions import Never
 
-from mpqp.core.instruction.gates.native_gates import SWAP
 from mpqp.core.instruction.measurement.measure import Measure
 from mpqp.core.instruction.measurement.pauli_string import (
     CommutingTypes,
@@ -24,7 +22,6 @@ class to define your observable, and a :class:`ExpectationMeasure` to perform
 )
 from mpqp.core.languages import Language
 from mpqp.tools.display import one_lined_repr
-from mpqp.tools.errors import NumberQubitsError
 from mpqp.tools.generics import Matrix
 from mpqp.tools.maths import is_diagonal, is_hermitian, is_power_of_two
 
@@ -38,8 +35,6 @@ class to define your observable, and a :class:`ExpectationMeasure` to perform
     from qiskit.quantum_info import SparsePauliOp
     from sympy import Expr
 
-    from mpqp.core.instruction.gates.custom_controlled_gate import Gate
-
 
 class Observable:
     """Class defining an observable, used for evaluating expectation values.
@@ -337,7 +332,7 @@ def to_other_language(
             return QLMObservable(self.nb_qubits, matrix=self.matrix)
         elif language == Language.BRAKET:
             if self._pauli_string:
-                from braket.circuits.observables import TensorProduct, Sum
+                from braket.circuits.observables import Sum, TensorProduct
 
                 obs = self.pauli_string.to_other_language(Language.BRAKET)
                 if isinstance(obs, TensorProduct):
@@ -469,8 +464,7 @@ def __init__(
             self.observables.append(new_obs)
 
         if targets is None:
-            self.targets = list(range(observable[0].nb_qubits))
-        self._check_targets_order()
+            self.targets = list(range(self.observables[0].nb_qubits))
 
     @property
     def nb_observables(self) -> int:
@@ -480,56 +474,6 @@ def nb_observables(self) -> int:
     def observables_labels(self) -> list[str]:
         return [o.label for o in self.observables if o.label is not None]
 
-    def _check_targets_order(self):
-        """Ensures target qubits are ordered and contiguous, rearranging them if
-        necessary (private)."""
-
-        if len(self.targets) == 0:
-            self._pre_measure: list[Gate] = []
-            return
-
-        if self.nb_qubits != self.observables[0].nb_qubits:
-            raise NumberQubitsError(
-                f"Target size {self.nb_qubits} doesn't match observable size "
-                f"{self.observables[0].nb_qubits}."
-            )
-
-        self._pre_measure: list[Gate] = []
-        """List of Gates added before the expectation measurement to correctly swap
-        target qubits when their are not ordered or contiguous."""
-        targets_is_ordered = all(
-            [self.targets[i] > self.targets[i - 1] for i in range(1, len(self.targets))]
-        )
-        tweaked_tgt = copy.copy(self.targets)
-        if (
-            max(tweaked_tgt) - min(tweaked_tgt) + 1 != len(tweaked_tgt)
-            or not targets_is_ordered
-        ):
-            warn(
-                "Non contiguous or non sorted observable target will introduce "
-                "additional CNOTs."
-            )
-
-            for t_index, target in enumerate(tweaked_tgt):  # sort the targets
-                min_index = tweaked_tgt.index(min(tweaked_tgt[t_index:]))
-                if t_index != min_index:
-                    self._pre_measure.append(SWAP(target, tweaked_tgt[min_index]))
-                    tweaked_tgt[t_index] = tweaked_tgt[min_index]
-                    tweaked_tgt[min_index] = target
-            for t_index, target in enumerate(tweaked_tgt):  # compact the targets
-                if t_index == 0:
-                    continue
-                if target != tweaked_tgt[t_index - 1] + 1:
-                    self._pre_measure.append(SWAP(target, tweaked_tgt[t_index - 1] + 1))
-                    tweaked_tgt[t_index] = tweaked_tgt[t_index - 1] + 1
-        self.rearranged_targets = tweaked_tgt
-        """Adjusted list of target qubits when they are not initially sorted and
-        contiguous."""
-
-    @property
-    def pre_measure(self) -> list[Gate]:
-        return self._pre_measure
-
     def get_pauli_grouping(self) -> list[list[PauliStringMonomial]]:
         """Return the grouped monomials of the Pauli string of the observable.
         The grouping is done according to the grouping method of the expectation
diff --git a/mpqp/core/instruction/measurement/pauli_string.py b/mpqp/core/instruction/measurement/pauli_string.py
index adc0fa4c..8fbd2fea 100644
--- a/mpqp/core/instruction/measurement/pauli_string.py
+++ b/mpqp/core/instruction/measurement/pauli_string.py
@@ -943,6 +943,56 @@ def is_diagonal(self) -> bool:
 
         return all([all([a == pI or a == pZ for a in m.atoms]) for m in self.monomials])
 
+    def rearrange(self, targets: list[int], copy: bool = True) -> "PauliString":
+        """
+        This function aims at reorderring a pauli string's monomials according to unordered targets.
+
+        Note: The targets must be contiguous, otherwise the algorithm won't work.
+        Args:
+            ps: The PauliString to reorder.
+            targets: The list of unordered targets.
+            copy: If set at True will deepcopy the initial string ps.
+
+        Examples:
+            >>> ps = pX @ pI
+            >>> print(ps.rearrange([1,0]))
+            pI@pX
+            >>> ps2 = pX @ pI + pI @ pX
+            >>> print(ps2.rearrange([1,0]))
+            pI@pX + pX@pI
+        """
+        from copy import deepcopy
+
+        if copy:
+            pauli = deepcopy(self)
+        else:
+            pauli = self
+
+        l = len(targets)
+        targets = deepcopy(targets)
+        rearranged = sorted(targets)
+        for index in range(l):
+            if targets[index] == index:
+                continue
+            shuffled_index = rearranged.index(targets[index])
+            for monom in pauli.monomials:
+                atoms = monom.atoms
+                atoms[shuffled_index], atoms[rearranged[index]] = (
+                    atoms[rearranged[index]],
+                    atoms[shuffled_index],
+                )
+            rearranged[index], rearranged[targets[index]] = (
+                rearranged[targets[index]],
+                rearranged[index],
+            )
+
+            i = targets.index(index)
+            targets[i], targets[index] = (
+                targets[index],
+                targets[i],
+            )
+        return pauli
+
 
 class PauliStringMonomial(PauliString):
     """Represents a monomial in a Pauli string, consisting of a coefficient and
diff --git a/mpqp/execution/providers/aws.py b/mpqp/execution/providers/aws.py
index 8c83f6a9..542430a7 100644
--- a/mpqp/execution/providers/aws.py
+++ b/mpqp/execution/providers/aws.py
@@ -187,11 +187,14 @@ def run_braket_observable(job: Job):
 
         if job.measure.pre_transpiled is None:
             grouping = job.measure.get_pauli_grouping()
+            pre_measure = [
+                QCircuit(find_qubitwise_rotations(group)) for group in grouping
+            ]
+            for circuit in pre_measure:
+                for instr in circuit.instructions:
+                    instr.targets = [job.measure.targets[t] for t in instr.targets]
             transpiled_pre_measures = [
-                QCircuit(find_qubitwise_rotations(group)).to_other_language(
-                    Language.BRAKET
-                )
-                for group in grouping
+                pre_m.to_other_language(Language.BRAKET) for pre_m in pre_measure
             ]
             eigenvalues = [
                 {monom.name: pauli_monomial_eigenvalues(monom) for monom in group}
@@ -226,7 +229,7 @@ def run_braket_observable(job: Job):
                 )
                 result = local_result.result()
                 assert isinstance(result, GateModelQuantumTaskResult)
-                length = 2**job.circuit.nb_qubits
+                length = 2**job.measure.nb_qubits
                 sorted_values: list[float] = []
                 for i in range(length):
                     binary_state = f"{bin(i)[2:].zfill(len(bin(length))- 3)}"
diff --git a/mpqp/execution/runner.py b/mpqp/execution/runner.py
index abc2bb3a..7c3a9239 100644
--- a/mpqp/execution/runner.py
+++ b/mpqp/execution/runner.py
@@ -55,48 +55,102 @@
 
 
 def adjust_measure(measure: ExpectationMeasure, circuit: QCircuit):
-    """We allow the measure to not span the entire circuit, but providers
+    """A measure can be incomplete and not span the entire circuit, but providers
     usually do not support this behavior. To make this work, we tweak the measure
     this function to match the expected behavior.
 
-    In order to do this, we add identity measures on the qubits not targeted by
-    the measure. In addition to this, some swaps are automatically added so the
-    the qubits measured are ordered and contiguous (though this is done in
-    :func:`generate_job`)
+    In order to do this, we place identity operators on the qubits not targeted by
+    the measure. If the targets are not ordered, each observable is first reordered so that
+    its local qubit order matches the sorted traget order. Pauli observables are directly embedded
+    on their target qubits, while matrix observables are padded with identity matrices
+    when the targets are ordered and contiguous, and are otherwise embedded through their pauli decomposition.
 
     Args:
         measure: The expectation measure, potentially incomplete.
-        circuit: The circuit to which will be added the potential swaps allowing
-            the user to get the expectation value of the qubits in an arbitrary
-            order (this part is not handled by this function).
+        circuit: The circuit defining the full qubit register.
 
     Returns:
-        The measure padded with identities before and after.
+        A measure targeting all circuit qubits, with observables embedded into the full register.
     """
     # TODO: use this only for specific provider
 
     if measure.targets == list(range(circuit.nb_qubits)):
         return measure
 
-    tweaked_observables = []
-    n_before = measure.rearranged_targets[0]
-    n_after = circuit.nb_qubits - measure.rearranged_targets[-1] - 1
+    nb_qubits = circuit.nb_qubits
+    targets = measure.targets
+
+    targets_is_ordered = all(
+        [targets[i] > targets[i - 1] for i in range(1, len(targets))]
+    )
+    if not targets_is_ordered:
+        ordered_targets = sorted(targets)
+        contiguous_targets = [targets.index(t) for t in ordered_targets]
+        for obs in measure.observables:
+            if (
+                obs._matrix is None  # pyright: ignore[reportPrivateUsage]
+                or measure.optimize_measurement
+            ):  # Order pauli string
+                from mpqp.tools.maths import rearrange_pauli_string
+
+                obs._pauli_string = (  # pyright: ignore[reportPrivateUsage]
+                    rearrange_pauli_string(obs.pauli_string, contiguous_targets)
+                )
+            else:  # Order the matrix
+                from mpqp.tools.maths import rearrange_matrix
+
+                obs.matrix = rearrange_matrix(obs.matrix, contiguous_targets)
+
+    targets_is_contiguous = len(targets) > 0 and (
+        targets[-1] - targets[0] + 1 == len(sorted(targets))
+    )
+
+    tweaked_observables: list[Observable] = []
+
     for obs in measure.observables:
-        if obs._pauli_string is not None:  # pyright: ignore[reportPrivateUsage]
-            from mpqp.measures import pI
+        from mpqp.core.instruction.measurement.pauli_string import (
+            PauliString,
+            PauliStringMonomial,
+        )
+        from mpqp.measures import pI
+
+        if (
+            obs._pauli_string is None  # pyright: ignore[reportPrivateUsage]
+            and targets_is_contiguous
+        ):
+            n_before = targets[0]
+            n_after = nb_qubits - targets[-1] - 1
+
+            full_matrix = obs.matrix
 
-            pauli = pI(n_before - 1) @ obs.pauli_string @ pI(n_after - 1)
-            tweaked_observables.append(Observable(pauli))
-        else:
             Id_before = np.eye(2**n_before)
             Id_after = np.eye(2**n_after)
+
+            if n_before > 0:
+                full_matrix = np.kron(Id_before, full_matrix)
+
+            if n_after > 0:
+                full_matrix = np.kron(full_matrix, Id_after)
+
             tweaked_observables.append(
                 Observable(
-                    np.kron(
-                        np.kron(Id_before, obs.matrix), Id_after
-                    )  # pyright: ignore[reportArgumentType]
+                    full_matrix, label=obs.label  # pyright: ignore[reportArgumentType]
                 )
             )
+            continue
+
+        pauli = obs.pauli_string
+        embedded = PauliString()
+
+        for mono in pauli.monomials:
+            full_register = [pI] * nb_qubits
+
+            for local_idx, target in enumerate(targets):
+                full_register[target] = mono.atoms[local_idx]
+
+            embedded += PauliStringMonomial(mono.coef, full_register)
+
+        tweaked_observables.append(Observable(embedded.simplify(), label=obs.label))
 
     tweaked_measure = ExpectationMeasure(
         tweaked_observables,
@@ -104,7 +158,9 @@ def adjust_measure(measure: ExpectationMeasure, circuit: QCircuit):
         measure.shots,
         measure.commuting_type,
         measure.grouping_method,
+        label=measure.label,
         optimize_measurement=measure.optimize_measurement,
+        optim_diagonal=measure.optim_diagonal,
     )
     return tweaked_measure
 
@@ -145,14 +201,21 @@ def generate_job(
             else:
                 job = Job(JobType.SAMPLE, circuit, device)
         elif isinstance(measurement, ExpectationMeasure):
-            m = adjust_measure(measurement, circuit)
-            c = circuit.without_measurements(deep_copy=False)
-            c.add(m)
-            job = Job(
-                JobType.OBSERVABLE,
-                c,
-                device,
-            )
+            if measurement.optimize_measurement and isinstance(device, AWSDevice):
+                job = Job(
+                    JobType.OBSERVABLE,
+                    circuit,
+                    device,
+                )
+            else:
+                m = adjust_measure(measurement, circuit)
+                c = circuit.without_measurements(deep_copy=False)
+                c.add(m)
+                job = Job(
+                    JobType.OBSERVABLE,
+                    c,
+                    device,
+                )
         else:
             raise NotImplementedError(
                 f"Measurement type {type(measurement)} not handled"
diff --git a/mpqp/tools/circuit.py b/mpqp/tools/circuit.py
index 16b1fbe2..19e12a68 100644
--- a/mpqp/tools/circuit.py
+++ b/mpqp/tools/circuit.py
@@ -13,6 +13,7 @@
     TOF,
     CRk,
     P,
+    OneQubitNoParamGate,
     Rk,
     RotationGate,
     Rx,
@@ -41,6 +42,7 @@ def random_circuit(
     gate_classes: Optional[Sequence[type[Gate]]] = None,
     nb_qubits: int = 5,
     nb_gates: Optional[int] = None,
+    use_all_qubits: bool = False,
     seed: Optional[int] = None,
 ):
     """This function creates a QCircuit with a specified number of qubits and gates.
@@ -86,6 +88,14 @@ def random_circuit(
     qcircuit = QCircuit(nb_qubits)
     for _ in range(nb_gates):
         qcircuit.add(random_gate(gate_classes, nb_qubits, rng))
+
+    if use_all_qubits:  # used in case we want to test braket
+        gates = np.array(
+            [g for g in NATIVE_GATES if issubclass(g, OneQubitNoParamGate)]
+        )
+
+        for i in range(nb_qubits):
+            qcircuit.add(rng.choice(gates)(i))
     return qcircuit
 
 
@@ -107,11 +117,12 @@ def statevector_from_random_circuit(
 
     Examples:
         >>> print(statevector_from_random_circuit(2, seed=123)) # doctest: +NORMALIZE_WHITESPACE
-        [0.70710678+0.j  0. -0.j  0.26893257-0.65396886j  0. -0.j ]
+        [0.70710678+0.j         0.        -0.j         0.26893257-0.65396886j
+         0.        -0.j        ]
     """
     from mpqp.execution import IBMDevice, Result, run
 
-    mpqp_circ = random_circuit(None, nb_qubits, None, seed)
+    mpqp_circ = random_circuit(None, nb_qubits, None, seed=seed)
     res = run(mpqp_circ, IBMDevice.AER_SIMULATOR_STATEVECTOR)
     if TYPE_CHECKING:
         assert isinstance(res, Result)
diff --git a/mpqp/tools/maths.py b/mpqp/tools/maths.py
index dfa5e053..ac341a12 100644
--- a/mpqp/tools/maths.py
+++ b/mpqp/tools/maths.py
@@ -11,10 +11,11 @@
 import numpy as np
 import numpy.typing as npt
 from scipy.linalg import inv, sqrtm
+from typeguard import typechecked
 
 if TYPE_CHECKING:
     from sympy import Expr
-
+    from mpqp.core.instruction.measurement.pauli_string import PauliString
     from mpqp.tools.generics import Matrix
 
 
@@ -478,3 +479,136 @@ def is_power_of_two(n: int) -> bool:
 
     """
     return n >= 1 and (n & (n - 1)) == 0
+
+
+@typechecked
+def rearrange_matrix(m: Matrix, targets: list[int], copy: bool = True) -> Matrix:
+    """Function to reorder the rows and columns of a matrix in order to change the targets of a gate.
+    The intended order for a gate is having continuous targets in growing order.
+
+    For example the targets for a 3 qubit gate should be [1,2,3], changing it for [3,2,1] would
+    reverse the effects on the qubits 3 and 1 (akin to a SWAP gate on those qubits).
+
+    Note: This function's goal is not to move around a gate in a circuit but to shuffle the targets in a sense.
+
+    Args:
+        m: The matrix for which we want to reorder the targets.
+        targets: The targets
+        copy: If True performs the copy of the matrix, to prevent overwriting the original matrix.
+
+    Returns:
+        The shuffled matrix according to the given targets.
+
+    Example:
+    >>> I = np.eye(2)
+    >>> Z = np.array([[1,0], [0,-1]])
+    >>> matrix = np.kron(I, Z)
+    >>> pprint(matrix)
+    [[1, 0 , 0, 0 ],
+     [0, -1, 0, 0 ],
+     [0, 0 , 1, 0 ],
+     [0, 0 , 0, -1]]
+    >>> pprint(rearrange_matrix(matrix, [1,0]))
+    [[1, 0, 0 , 0 ],
+     [0, 1, 0 , 0 ],
+     [0, 0, -1, 0 ],
+     [0, 0, 0 , -1]]
+    """
+    from copy import deepcopy
+
+    if copy:
+        matrix = deepcopy(m)
+    else:
+        matrix = m
+    l = len(targets)
+    targets = deepcopy(targets)
+    shuffled = sorted(targets)
+    for index in range(l - 1):
+        if targets[index] == index:
+            continue
+        # If no swaps happened of the target then shuffled_index = targets[index]
+        shuffled_index = shuffled.index(targets[index])
+
+        i = 1 << (l - 1 - shuffled_index)
+        j = 1 << (l - 1 - index)
+        for change in range(1 << l):
+            current = bin(change)[2:].zfill(l)
+            if current[shuffled_index - l] == "0" and current[index - l] == "1":
+                current = int(current, 2)
+                conjugate = current + i - j
+                for k in range(len(matrix)):
+                    hold = matrix[k][current]
+                    matrix[k][current] = matrix[k][conjugate]
+                    matrix[k][conjugate] = hold
+
+                for k in range(len(matrix)):
+                    hold = matrix[current][k]
+                    matrix[current][k] = matrix[conjugate][k]
+                    matrix[conjugate][k] = hold
+
+        # keeps tracks of the position of the targets in the matrix
+
+        shuffled[index], shuffled[targets[index]] = (
+            shuffled[targets[index]],
+            shuffled[index],
+        )
+
+        i = targets.index(index)
+        targets[i], targets[index] = (
+            targets[index],
+            targets[i],
+        )
+    return matrix
+
+
+def rearrange_pauli_string(
+    ps: PauliString, targets: list[int], copy: bool = True
+) -> PauliString:
+    """
+    This function aims at reorderring a pauli string's monomials according to unordered targets.
+
+    Note: The targets must be contiguous, otherwise the algorithm won't work.
+    Args:
+        ps: The PauliString to reorder.
+        targets: The list of unordered targets.
+        copy: If set at True will deepcopy the initial string ps.
+
+    Examples:
+        >>> ps = pX @ pI
+        >>> print(rearrange_pauli_string(ps, [1,0]))
+        pI@pX
+        >>> ps2 = pX @ pI + pI @ pX
+        >>> print(rearrange_pauli_string(ps2, [1,0]))
+        pI@pX + pX@pI
+    """
+    from copy import deepcopy
+
+    if copy:
+        pauli = deepcopy(ps)
+    else:
+        pauli = ps
+
+    l = len(targets)
+    targets = deepcopy(targets)
+    shuffled = sorted(targets)
+    for index in range(l):
+        if targets[index] == index:
+            continue
+        shuffled_index = shuffled.index(targets[index])
+        for monom in pauli.monomials:
+            atoms = monom.atoms
+            atoms[shuffled_index], atoms[shuffled[index]] = (
+                atoms[shuffled[index]],
+                atoms[shuffled_index],
+            )
+        shuffled[index], shuffled[targets[index]] = (
+            shuffled[targets[index]],
+            shuffled[index],
+        )
+
+        i = targets.index(index)
+        targets[i], targets[index] = (
+            targets[index],
+            targets[i],
+        )
+    return pauli
diff --git a/tests/core/instruction/gates/test_custom_gate.py b/tests/core/instruction/gates/test_custom_gate.py
index 37c6c3aa..70d91b45 100644
--- a/tests/core/instruction/gates/test_custom_gate.py
+++ b/tests/core/instruction/gates/test_custom_gate.py
@@ -64,7 +64,6 @@ def exec_random_orthogonal_matrix(circ_size: int, device: AvailableDevice):
     result = run(c, device)
 
     # we reduce the precision because of approximation errors coming from CustomGate usage
-    assert isinstance(result, Result)
     assert matrix_eq(result.amplitudes, exp_state_vector, 1e-5, 1e-5)
 
 
@@ -124,7 +123,9 @@ def test_custom_gate_with_random_circuit_qiskit(circ_size: int):
 @pytest.mark.provider("braket")
 @pytest.mark.parametrize("circ_size", range(1, 6))
 def test_custom_gate_with_random_circuit_braket(circ_size: int):
-    exec_custom_gate_with_random_circuit(circ_size, AWSDevice.BRAKET_LOCAL_SIMULATOR)
+    exec_custom_gate_with_random_circuit(
+        circ_size, AWSDevice.BRAKET_LOCAL_SIMULATOR, True
+    )
 
 
 @pytest.mark.provider("cirq")
@@ -139,8 +140,10 @@ def test_custom_gate_with_random_circuit_myqlm(circ_size: int):
     exec_custom_gate_with_random_circuit(circ_size, ATOSDevice.MYQLM_PYLINALG)
 
 
-def exec_custom_gate_with_random_circuit(circ_size: int, device: AvailableDevice):
-    random_circ = random_circuit(nb_qubits=circ_size)
+def exec_custom_gate_with_random_circuit(
+    circ_size: int, device: AvailableDevice, use_all_qubits: bool = False
+):
+    random_circ = random_circuit(nb_qubits=circ_size, use_all_qubits=use_all_qubits)
     matrix = random_circ.to_matrix()
     custom_gate_circ = QCircuit([CustomGate(matrix, list(range(circ_size)))])
 
diff --git a/tests/core/instruction/measurement/test_expectation_value.py b/tests/core/instruction/measurement/test_expectation_value.py
index 82ae65e3..9ba7a3e5 100644
--- a/tests/core/instruction/measurement/test_expectation_value.py
+++ b/tests/core/instruction/measurement/test_expectation_value.py
@@ -24,23 +24,6 @@ def test_expectation_measure_right_targets(targets: list[int]):
         ExpectationMeasure(obs, targets)
 
 
-@pytest.mark.parametrize(
-    "targets, expected_swaps",
-    [
-        ([1, 3, 4], [{2, 3}, {3, 4}]),
-        ([1, 0, 2], [{1, 0}]),
-        ([2, 0, 3], [{0, 2}, {1, 2}, {2, 3}]),
-    ],
-)
-def test_expectation_measure_wrong_targets(
-    targets: list[int], expected_swaps: list[tuple[int, int]]
-):
-    obs = Observable(np.diag([1] * 2 ** len(targets)))
-    with pytest.warns(UserWarning):
-        measure = ExpectationMeasure(obs, targets)
-    assert [set(swap.targets) for swap in measure.pre_measure] == expected_swaps
-
-
 # TODO: complete this
 @pytest.fixture
 def list_to_cirq_pauli() -> (
diff --git a/tests/core/instruction/measurement/test_pauli_string.py b/tests/core/instruction/measurement/test_pauli_string.py
index 4afc9c76..5ec6b295 100644
--- a/tests/core/instruction/measurement/test_pauli_string.py
+++ b/tests/core/instruction/measurement/test_pauli_string.py
@@ -826,3 +826,18 @@ def test_pauli_monomial_from_atom(
 ):
     result = atom(prefix) if postfix is None else atom(prefix, postfix)
     assert result == expected_ps
+
+
+@pytest.mark.parametrize(
+    ("ps", "targets", "expected"),
+    [
+        (pX @ pY, [1, 0], pY @ pX),
+        (pX @ pI @ pZ, [1, 0, 2], pI @ pX @ pZ),
+        (pX @ pY @ pZ, [2, 0, 1], pY @ pZ @ pX),
+    ],
+)
+def test_rearrange_pauli_string(
+    ps: PauliString, targets: list[int], expected: PauliString
+):
+
+    assert ps.rearrange(targets) == expected
diff --git a/tests/test_doc.py b/tests/test_doc.py
index cc1179ea..dd561b69 100644
--- a/tests/test_doc.py
+++ b/tests/test_doc.py
@@ -133,6 +133,8 @@
     rand_product_local_unitaries,
     rand_unitary_2x2_matrix,
     rand_unitary_matrix,
+    rearrange_matrix,
+    rearrange_pauli_string,
 )
 from mpqp.tools.operators import *
 from mpqp.tools.pauli_grouping import CommutingTypes, pauli_grouping_greedy
diff --git a/tests/tools/test_math.py b/tests/tools/test_math.py
index f1b3df9e..18bb0f46 100644
--- a/tests/tools/test_math.py
+++ b/tests/tools/test_math.py
@@ -3,7 +3,7 @@
 from sympy import symbols
 
 from mpqp.tools.generics import Matrix
-from mpqp.tools.maths import is_hermitian, rand_hermitian_matrix
+from mpqp.tools.maths import is_hermitian, rand_hermitian_matrix, rand_unitary_matrix
 
 x = symbols("x", real=True)
 
@@ -25,3 +25,22 @@ def test_is_hermitian(matrix: Matrix, isHermitian: bool):
 
 def test_rand_hermitian():
     assert is_hermitian(rand_hermitian_matrix(3))
+
+
+@pytest.mark.parametrize(
+    ("matrix", "targets"),
+    [
+        (rand_unitary_matrix(4), [1, 0]),
+        (rand_unitary_matrix(8), [1, 0, 2]),
+        (rand_unitary_matrix(8), [2, 0, 1]),
+    ],
+)
+def test_rearrange_matrix(matrix: Matrix, targets: list[int]):
+    from mpqp.gates import CustomGate
+    from mpqp.tools.maths import rearrange_matrix, matrix_eq
+    from mpqp import QCircuit
+
+    g = CustomGate(matrix, targets)
+    m = rearrange_matrix(matrix, targets)
+    g2 = CustomGate(m, sorted(targets))
+    assert matrix_eq(QCircuit([g]).to_matrix(), g2.to_matrix())