From 1ff915bc27bbcb778cf77757c047e831f8657f6c Mon Sep 17 00:00:00 2001
From: Alejandro Sopena <alexsopena15@gmail.com>
Date: Mon, 16 Oct 2023 17:26:35 +0200
Subject: [PATCH 01/47] clifford sampling

---
 src/qibo/models/error_mitigation.py | 177 ++++++++++++++++++++++++++++
 1 file changed, 177 insertions(+)

diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index f0d333203d..c17a845655 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -1,5 +1,7 @@
 """Error Mitigation Methods."""
 
+from itertools import product
+
 import numpy as np
 from scipy.optimize import curve_fit
 
@@ -595,3 +597,178 @@ def apply_randomized_readout_mitigation(
         results[j]._frequencies = freq_sum
 
     return results
+
+
+def sample_clifford_training_circuit(
+    circuit,
+    backend=None,
+):
+    """Samples a training circuit for CDR by susbtituting some of the non-Clifford gates.
+
+    Args:
+        circuit (:class:`qibo.models.Circuit`): circuit to sample from,
+            decomposed in ``RX(pi/2)``, ``X``, ``CNOT`` and ``RZ`` gates.
+        replacement_gates (list, optional): candidates for the substitution of the
+            non-Clifford gates. The ``list`` should be composed by ``tuples`` of the
+            form (``gates.XYZ``, ``kwargs``). For example, phase gates are used by default:
+            ``list((RZ, {'theta':0}), (RZ, {'theta':pi/2}), (RZ, {'theta':pi}), (RZ, {'theta':3*pi/2}))``.
+        sigma (float, optional): standard devation of the Gaussian distribution used for sampling.
+        backend (:class:`qibo.backends.abstract.Backend`, optional): Calculation engine.
+
+    Returns:
+        :class:`qibo.models.Circuit`: The sampled circuit.
+    """
+    from qibo.quantum_info import random_clifford
+
+    if backend is None:  # pragma: no cover
+        backend = GlobalBackend()
+
+    # Find all the non-Clifford gates
+    gates_to_replace = []
+    for i, gate in enumerate(circuit.queue):
+        if gate.clifford is False and isinstance(gate, gates.M) is False:
+            gates_to_replace.append([i, gate])
+    gates_to_replace = np.array(gates_to_replace, dtype=object)
+
+    if len(gates_to_replace) == 0:
+        raise_error(ValueError, "No non-Clifford gate found, no circuit sampled.")
+
+    # Build the training circuit by substituting the sampled gates
+    sampled_circuit = circuit.__class__(**circuit.init_kwargs)
+
+    for i, gate in enumerate(circuit.queue):
+        if isinstance(gate, gates.M):
+            gate_rand = gates.Unitary(
+                random_clifford(1, backend=backend, return_circuit=False),
+                gate.qubits[0],
+            )
+            gate_rand.clifford = True
+
+            sampled_circuit.add(gate_rand)
+            sampled_circuit.add(gate)
+
+        else:
+            if i in gates_to_replace[:, 0]:
+                gate = gates.Unitary(
+                    random_clifford(1, backend=backend, return_circuit=False),
+                    gate.qubits[0],
+                )
+                gate.clifford = True
+            sampled_circuit.add(gate)
+
+    return sampled_circuit
+
+
+def escircuit(circuit, obs, backend=None):
+    from qibo.quantum_info import comp_basis_to_pauli, random_clifford, vectorization
+
+    if backend is None:  # pragma: no cover
+        backend = GlobalBackend()
+
+    circ_cliff = sample_clifford_training_circuit(circuit, backend=backend)
+
+    c_unitary = circ_cliff.unitary(backend=backend)
+    nqubits = circ_cliff.nqubits
+    U_c2p = comp_basis_to_pauli(nqubits, backend=backend)
+    obs_liouville = vectorization(
+        np.transpose(np.conjugate(c_unitary)) @ obs.matrix @ c_unitary, order="row"
+    )
+    obs_pauli_liouville = U_c2p @ obs_liouville
+    index = np.where(abs(obs_pauli_liouville) >= 1e-5)[0][0]
+    obs1 = list(product(["I", "X", "Y", "Z"], repeat=nqubits))[index]
+
+    paulis = {
+        "I": gates.I(0).matrix(),
+        "X": gates.X(0).matrix(),
+        "Y": gates.Y(0).matrix(),
+        "Z": gates.Z(0).matrix(),
+    }
+
+    adjust_gates = []
+    for i in range(nqubits):
+        obs_i = paulis[obs1[i]]
+        R = paulis["I"]
+        while np.any(abs(obs_i - paulis["Z"]) > 1e-5) and np.any(
+            abs(obs_i - paulis["I"]) > 1e-5
+        ):
+            R = random_clifford(1, backend=backend, return_circuit=False)
+            obs_i = np.conjugate(np.transpose(R)) @ paulis[obs1[i]] @ R
+
+        adjust_gate = gates.Unitary(R, i)
+        adjust_gate.clifford = True
+        adjust_gates.append(adjust_gate)
+
+    circ_cliff1 = circ_cliff.__class__(**circ_cliff.init_kwargs)
+
+    for gate in adjust_gates:
+        circ_cliff1.add(gate)
+
+    for gate in circ_cliff.queue:
+        circ_cliff1.add(gate)
+
+    return circ_cliff1.fuse(max_qubits=1), circ_cliff, adjust_gates
+
+
+def mit_obs(
+    circuit,
+    observable,
+    noise_model,
+    nshots=int(1e4),
+    n_training_samples=10,
+    backend=None,
+):
+    if backend is None:  # pragma: no cover
+        backend = GlobalBackend()
+
+    training_circs = [
+        escircuit(circuit, observable, backend=backend)[0]
+        for _ in range(n_training_samples)
+    ]
+
+    data = {"exact": {"-1": [], "1": []}, "noisy": {"-1": [], "1": []}}
+
+    a_list = {"-1": [], "1": []}
+    for c in training_circs:
+        exp = c(nshots=nshots).expectation_from_samples(observable)
+        # state = c().state()
+        # exp = obs.expectation(state)
+        if noise_model is not None and backend.name != "qibolab":
+            c = noise_model.apply(c)
+            c.density_matrix = True
+        circuit_result = backend.execute_circuit(c, nshots=nshots)
+        exp_noisy = circuit_result.expectation_from_samples(observable)
+        # state = c_noisy().state()
+        # exp_noisy = obs.expectation(state)
+
+        if exp > 0:
+            data["exact"]["1"].append(exp)
+            data["noisy"]["1"].append(exp_noisy)
+            a_list["1"].append(1 - exp_noisy / exp)
+        else:
+            data["exact"]["-1"].append(exp)
+            data["noisy"]["-1"].append(exp_noisy)
+            a_list["-1"].append(1 - exp_noisy / exp)
+
+    a_std_1 = np.std(a_list["1"])
+    a_std_m1 = np.std(a_list["-1"])
+    a = (np.sum(a_list["1"]) + np.sum(a_list["-1"])) / n_training_samples
+    a_std = np.sqrt(
+        (a_std_1**2 * len(a_list["1"]) + a_std_m1**2 * len(a_list["-1"]))
+        / n_training_samples
+    )
+
+    if noise_model is not None and backend.name != "qibolab":
+        circuit = noise_model.apply(circuit)
+        circuit.density_matrix = True
+    circuit_result = backend.execute_circuit(circuit, nshots=nshots)
+    exp_noisy = circuit_result.expectation_from_samples(observable)
+    # exp_noisy = obs.expectation(c_noisy().state())
+    exp_mit = (1 - a) * exp_noisy / ((1 - a) ** 2 + a_std**2)
+    exp_mit_std = (
+        a_std
+        * abs(exp_noisy)
+        * abs((1 - a) ** 2 - a_std**2)
+        / ((1 - a) ** 2 + a_std**2) ** 2
+    )
+
+    return exp_mit, exp_mit_std, a, a_std, a_list, data

From 65630298e69ce2871d4daef76ff927ddcde91f6b Mon Sep 17 00:00:00 2001
From: Alejandro Sopena <alexsopena15@gmail.com>
Date: Mon, 16 Oct 2023 22:29:02 +0200
Subject: [PATCH 02/47] output compatible with cdr

---
 src/qibo/models/error_mitigation.py | 20 +++++++++++++-------
 1 file changed, 13 insertions(+), 7 deletions(-)

diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index c17a845655..6af4c09334 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -655,11 +655,10 @@ def sample_clifford_training_circuit(
                 )
                 gate.clifford = True
             sampled_circuit.add(gate)
-
     return sampled_circuit
 
 
-def escircuit(circuit, obs, backend=None):
+def escircuit(circuit, obs, fuse=True, backend=None):
     from qibo.quantum_info import comp_basis_to_pauli, random_clifford, vectorization
 
     if backend is None:  # pragma: no cover
@@ -706,7 +705,10 @@ def escircuit(circuit, obs, backend=None):
     for gate in circ_cliff.queue:
         circ_cliff1.add(gate)
 
-    return circ_cliff1.fuse(max_qubits=1), circ_cliff, adjust_gates
+    if fuse:
+        circ_cliff1 = circ_cliff1.fuse(max_qubits=1)
+
+    return circ_cliff1, circ_cliff, adjust_gates
 
 
 def mit_obs(
@@ -715,6 +717,7 @@ def mit_obs(
     noise_model,
     nshots=int(1e4),
     n_training_samples=10,
+    full_output=False,
     backend=None,
 ):
     if backend is None:  # pragma: no cover
@@ -725,7 +728,7 @@ def mit_obs(
         for _ in range(n_training_samples)
     ]
 
-    data = {"exact": {"-1": [], "1": []}, "noisy": {"-1": [], "1": []}}
+    data = {"noise-free": {"-1": [], "1": []}, "noisy": {"-1": [], "1": []}}
 
     a_list = {"-1": [], "1": []}
     for c in training_circs:
@@ -741,11 +744,11 @@ def mit_obs(
         # exp_noisy = obs.expectation(state)
 
         if exp > 0:
-            data["exact"]["1"].append(exp)
+            data["noise-free"]["1"].append(exp)
             data["noisy"]["1"].append(exp_noisy)
             a_list["1"].append(1 - exp_noisy / exp)
         else:
-            data["exact"]["-1"].append(exp)
+            data["noise-free"]["-1"].append(exp)
             data["noisy"]["-1"].append(exp_noisy)
             a_list["-1"].append(1 - exp_noisy / exp)
 
@@ -771,4 +774,7 @@ def mit_obs(
         / ((1 - a) ** 2 + a_std**2) ** 2
     )
 
-    return exp_mit, exp_mit_std, a, a_std, a_list, data
+    if full_output:
+        return exp_mit, exp_mit_std, a, a_std, a_list, data
+
+    return exp_mit

From 40251cab3e87af887af08473c2cead4605b4fb84 Mon Sep 17 00:00:00 2001
From: Alejandro Sopena <alexsopena15@gmail.com>
Date: Tue, 17 Oct 2023 13:03:24 +0200
Subject: [PATCH 03/47] do not fuse I gates

---
 src/qibo/gates/special.py           | 3 ++-
 src/qibo/models/circuit.py          | 2 +-
 src/qibo/models/error_mitigation.py | 3 +--
 3 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/qibo/gates/special.py b/src/qibo/gates/special.py
index 236f19d72d..8968e82ba2 100644
--- a/src/qibo/gates/special.py
+++ b/src/qibo/gates/special.py
@@ -1,5 +1,6 @@
 from qibo.backends import GlobalBackend
 from qibo.gates.abstract import SpecialGate
+from qibo.gates.gates import I
 from qibo.gates.measurements import M
 
 
@@ -54,7 +55,7 @@ def __init__(self, *q):
     def from_gate(cls, gate):
         fgate = cls(*gate.qubits)
         fgate.append(gate)
-        if isinstance(gate, (M, SpecialGate)):
+        if isinstance(gate, (M, SpecialGate, I)):
             # special gates do not participate in fusion
             fgate.marked = True
         return fgate
diff --git a/src/qibo/models/circuit.py b/src/qibo/models/circuit.py
index 3e5d132b1f..f333736bd4 100644
--- a/src/qibo/models/circuit.py
+++ b/src/qibo/models/circuit.py
@@ -74,7 +74,7 @@ def from_fused(self):
                     queue.append(gate.gates[0])
                 else:
                     queue.append(gate)
-            elif isinstance(gate.gates[0], (gates.SpecialGate, gates.M)):
+            elif isinstance(gate.gates[0], (gates.SpecialGate, gates.M, gates.I)):
                 # special gates are marked by default so we need
                 # to add them manually
                 queue.append(gate.gates[0])
diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index 6af4c09334..b2b06f8925 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -737,7 +737,6 @@ def mit_obs(
         # exp = obs.expectation(state)
         if noise_model is not None and backend.name != "qibolab":
             c = noise_model.apply(c)
-            c.density_matrix = True
         circuit_result = backend.execute_circuit(c, nshots=nshots)
         exp_noisy = circuit_result.expectation_from_samples(observable)
         # state = c_noisy().state()
@@ -760,9 +759,9 @@ def mit_obs(
         / n_training_samples
     )
 
+    circuit.fuse(max_qubits=1)
     if noise_model is not None and backend.name != "qibolab":
         circuit = noise_model.apply(circuit)
-        circuit.density_matrix = True
     circuit_result = backend.execute_circuit(circuit, nshots=nshots)
     exp_noisy = circuit_result.expectation_from_samples(observable)
     # exp_noisy = obs.expectation(c_noisy().state())

From 5a261e487b6698f4d7d383ff2b6bd5b268b637de Mon Sep 17 00:00:00 2001
From: Alejandro Sopena <alexsopena15@gmail.com>
Date: Tue, 17 Oct 2023 13:39:22 +0200
Subject: [PATCH 04/47] fix fuse

---
 src/qibo/models/error_mitigation.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index b2b06f8925..558e3f3a62 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -759,7 +759,7 @@ def mit_obs(
         / n_training_samples
     )
 
-    circuit.fuse(max_qubits=1)
+    circuit = circuit.fuse(max_qubits=1)
     if noise_model is not None and backend.name != "qibolab":
         circuit = noise_model.apply(circuit)
     circuit_result = backend.execute_circuit(circuit, nshots=nshots)

From 204920796e579be4fcc569f7aa46d4a8b96c72bb Mon Sep 17 00:00:00 2001
From: Alejandro Sopena <alexsopena15@gmail.com>
Date: Fri, 20 Oct 2023 15:06:45 +0200
Subject: [PATCH 05/47] transpiling circuit

---
 src/qibo/models/error_mitigation.py | 35 +++++++++++++++++++++++++----
 1 file changed, 31 insertions(+), 4 deletions(-)

diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index 558e3f3a62..a20480ea03 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -658,6 +658,25 @@ def sample_clifford_training_circuit(
     return sampled_circuit
 
 
+def transpile_circ(circuit, qubit_map, backend):
+    if backend.name == "qibolab":
+        new_c = circuit.__class__(5)
+        for gate in circuit.queue:
+            qubits = [qubit_map[j] for j in gate.qubits]
+            if isinstance(gate, gates.M):
+                new_gate = gates.M(*tuple(qubits), **gate.init_kwargs)
+                new_gate.result = gate.result
+                new_c.add(new_gate)
+            elif isinstance(gate, gates.I):
+                new_c.add(gate.__class__(*tuple(qubits), **gate.init_kwargs))
+            else:
+                matrix = gate.matrix()
+                new_c.add(gates.Unitary(matrix, *tuple(qubits), **gate.init_kwargs))
+        return new_c
+    else:
+        return circuit
+
+
 def escircuit(circuit, obs, fuse=True, backend=None):
     from qibo.quantum_info import comp_basis_to_pauli, random_clifford, vectorization
 
@@ -714,7 +733,8 @@ def escircuit(circuit, obs, fuse=True, backend=None):
 def mit_obs(
     circuit,
     observable,
-    noise_model,
+    qubit_map=None,
+    noise_model=None,
     nshots=int(1e4),
     n_training_samples=10,
     full_output=False,
@@ -723,6 +743,9 @@ def mit_obs(
     if backend is None:  # pragma: no cover
         backend = GlobalBackend()
 
+    if qubit_map is None:
+        qubit_map = list(range(circuit.nqubits))
+
     training_circs = [
         escircuit(circuit, observable, backend=backend)[0]
         for _ in range(n_training_samples)
@@ -732,13 +755,16 @@ def mit_obs(
 
     a_list = {"-1": [], "1": []}
     for c in training_circs:
-        exp = c(nshots=nshots).expectation_from_samples(observable)
+        circuit_result = c(nshots=nshots)
+        exp = observable.expectation_from_samples(circuit_result.frequencies())
         # state = c().state()
         # exp = obs.expectation(state)
         if noise_model is not None and backend.name != "qibolab":
             c = noise_model.apply(c)
+
+        c = transpile_circ(c, qubit_map, backend)
         circuit_result = backend.execute_circuit(c, nshots=nshots)
-        exp_noisy = circuit_result.expectation_from_samples(observable)
+        exp_noisy = observable.expectation_from_samples(circuit_result.frequencies())
         # state = c_noisy().state()
         # exp_noisy = obs.expectation(state)
 
@@ -762,8 +788,9 @@ def mit_obs(
     circuit = circuit.fuse(max_qubits=1)
     if noise_model is not None and backend.name != "qibolab":
         circuit = noise_model.apply(circuit)
+    circuit = transpile_circ(circuit, qubit_map, backend)
     circuit_result = backend.execute_circuit(circuit, nshots=nshots)
-    exp_noisy = circuit_result.expectation_from_samples(observable)
+    exp_noisy = observable.expectation_from_samples(circuit_result.frequencies())
     # exp_noisy = obs.expectation(c_noisy().state())
     exp_mit = (1 - a) * exp_noisy / ((1 - a) ** 2 + a_std**2)
     exp_mit_std = (

From cdf56471454df17e939ae786117f4755debf4e29 Mon Sep 17 00:00:00 2001
From: Alejandro Sopena <alexsopena15@gmail.com>
Date: Fri, 20 Oct 2023 15:07:45 +0200
Subject: [PATCH 06/47] fix

---
 src/qibo/models/error_mitigation.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index a20480ea03..738aedb22a 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -660,7 +660,7 @@ def sample_clifford_training_circuit(
 
 def transpile_circ(circuit, qubit_map, backend):
     if backend.name == "qibolab":
-        new_c = circuit.__class__(5)
+        new_c = circuit.__class__(backend.platform.nqubits)
         for gate in circuit.queue:
             qubits = [qubit_map[j] for j in gate.qubits]
             if isinstance(gate, gates.M):

From c1fc36bc59b971531265123d5850999196093774 Mon Sep 17 00:00:00 2001
From: Alejandro <alexsopena15@gmail.com>
Date: Fri, 20 Oct 2023 21:26:11 +0400
Subject: [PATCH 07/47] fix transpiling

---
 src/qibo/models/error_mitigation.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index 738aedb22a..a77903c20e 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -659,6 +659,7 @@ def sample_clifford_training_circuit(
 
 
 def transpile_circ(circuit, qubit_map, backend):
+    from qibolab.transpilers.unitary_decompositions import u3_decomposition
     if backend.name == "qibolab":
         new_c = circuit.__class__(backend.platform.nqubits)
         for gate in circuit.queue:
@@ -671,7 +672,7 @@ def transpile_circ(circuit, qubit_map, backend):
                 new_c.add(gate.__class__(*tuple(qubits), **gate.init_kwargs))
             else:
                 matrix = gate.matrix()
-                new_c.add(gates.Unitary(matrix, *tuple(qubits), **gate.init_kwargs))
+                new_c.add(gates.U3(qubits[0], *u3_decomposition(matrix)))
         return new_c
     else:
         return circuit

From 6c49e7123137ce137df7c58dfae8c9d9cd29e8cb Mon Sep 17 00:00:00 2001
From: Alejandro Sopena <alexsopena15@gmail.com>
Date: Fri, 20 Oct 2023 20:56:58 +0200
Subject: [PATCH 08/47] fix identities

---
 src/qibo/models/error_mitigation.py | 19 +++++++++++++++----
 1 file changed, 15 insertions(+), 4 deletions(-)

diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index a77903c20e..dad2870686 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -637,7 +637,7 @@ def sample_clifford_training_circuit(
     sampled_circuit = circuit.__class__(**circuit.init_kwargs)
 
     for i, gate in enumerate(circuit.queue):
-        if isinstance(gate, gates.M):
+        if gate.name == "id_end":  # isinstance(gate, gates.M):
             gate_rand = gates.Unitary(
                 random_clifford(1, backend=backend, return_circuit=False),
                 gate.qubits[0],
@@ -660,6 +660,7 @@ def sample_clifford_training_circuit(
 
 def transpile_circ(circuit, qubit_map, backend):
     from qibolab.transpilers.unitary_decompositions import u3_decomposition
+
     if backend.name == "qibolab":
         new_c = circuit.__class__(backend.platform.nqubits)
         for gate in circuit.queue:
@@ -719,11 +720,19 @@ def escircuit(circuit, obs, fuse=True, backend=None):
 
     circ_cliff1 = circ_cliff.__class__(**circ_cliff.init_kwargs)
 
-    for gate in adjust_gates:
-        circ_cliff1.add(gate)
+    # for gate in adjust_gates:
+    #     circ_cliff1.add(gate)
 
+    # for gate in circ_cliff.queue:
+    #     circ_cliff1.add(gate)
+    j = 0
     for gate in circ_cliff.queue:
-        circ_cliff1.add(gate)
+        if gate.name == "id_init":
+            circ_cliff1.add(gate)
+            circ_cliff1.add(adjust_gates[j])
+            j += 1
+        else:
+            circ_cliff1.add(gate)
 
     if fuse:
         circ_cliff1 = circ_cliff1.fuse(max_qubits=1)
@@ -764,6 +773,7 @@ def mit_obs(
             c = noise_model.apply(c)
 
         c = transpile_circ(c, qubit_map, backend)
+        print(c.draw())
         circuit_result = backend.execute_circuit(c, nshots=nshots)
         exp_noisy = observable.expectation_from_samples(circuit_result.frequencies())
         # state = c_noisy().state()
@@ -787,6 +797,7 @@ def mit_obs(
     )
 
     circuit = circuit.fuse(max_qubits=1)
+    print(circuit.draw())
     if noise_model is not None and backend.name != "qibolab":
         circuit = noise_model.apply(circuit)
     circuit = transpile_circ(circuit, qubit_map, backend)

From b1b3971841bddf84aa77b106bba006adefb20227 Mon Sep 17 00:00:00 2001
From: Alejandro <alexsopena15@gmail.com>
Date: Fri, 20 Oct 2023 23:47:32 +0400
Subject: [PATCH 09/47] remove prints

---
 src/qibo/models/error_mitigation.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index dad2870686..db208c6837 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -773,7 +773,6 @@ def mit_obs(
             c = noise_model.apply(c)
 
         c = transpile_circ(c, qubit_map, backend)
-        print(c.draw())
         circuit_result = backend.execute_circuit(c, nshots=nshots)
         exp_noisy = observable.expectation_from_samples(circuit_result.frequencies())
         # state = c_noisy().state()
@@ -797,7 +796,6 @@ def mit_obs(
     )
 
     circuit = circuit.fuse(max_qubits=1)
-    print(circuit.draw())
     if noise_model is not None and backend.name != "qibolab":
         circuit = noise_model.apply(circuit)
     circuit = transpile_circ(circuit, qubit_map, backend)

From 57db2ec77919bf3ed2b2e78f8d23e296b9a033a7 Mon Sep 17 00:00:00 2001
From: Alejandro <alexsopena15@gmail.com>
Date: Sun, 22 Oct 2023 21:27:17 +0400
Subject: [PATCH 10/47] add minus for iqm5q

---
 src/qibo/models/error_mitigation.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index db208c6837..6385381a1a 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -774,7 +774,7 @@ def mit_obs(
 
         c = transpile_circ(c, qubit_map, backend)
         circuit_result = backend.execute_circuit(c, nshots=nshots)
-        exp_noisy = observable.expectation_from_samples(circuit_result.frequencies())
+        exp_noisy = observable.expectation_from_samples(circuit_result.frequencies()) #Add minus for iqm
         # state = c_noisy().state()
         # exp_noisy = obs.expectation(state)
 
@@ -800,7 +800,7 @@ def mit_obs(
         circuit = noise_model.apply(circuit)
     circuit = transpile_circ(circuit, qubit_map, backend)
     circuit_result = backend.execute_circuit(circuit, nshots=nshots)
-    exp_noisy = observable.expectation_from_samples(circuit_result.frequencies())
+    exp_noisy = observable.expectation_from_samples(circuit_result.frequencies()) # Add - for iqm
     # exp_noisy = obs.expectation(c_noisy().state())
     exp_mit = (1 - a) * exp_noisy / ((1 - a) ** 2 + a_std**2)
     exp_mit_std = (

From 332d88b78169259908bd90877ed58033a3b77833 Mon Sep 17 00:00:00 2001
From: Alejandro <alexsopena15@gmail.com>
Date: Wed, 25 Oct 2023 20:31:02 +0400
Subject: [PATCH 11/47] readout

---
 src/qibo/models/error_mitigation.py | 45 ++++++++++++++++++++++-------
 1 file changed, 34 insertions(+), 11 deletions(-)

diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index 6385381a1a..84d5bf502c 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -463,7 +463,7 @@ def vnCDR(
     return mit_val
 
 
-def calibration_matrix(nqubits, noise_model=None, nshots: int = 1000, backend=None):
+def calibration_matrix(nqubits, qubit_map, noise_model=None, nshots: int = 1000, backend=None):
     """Computes the calibration matrix for readout mitigation.
 
     Args:
@@ -485,7 +485,7 @@ def calibration_matrix(nqubits, noise_model=None, nshots: int = 1000, backend=No
         backend = GlobalBackend()
 
     matrix = np.zeros((2**nqubits, 2**nqubits))
-
+    from qibo.config import log
     for i in range(2**nqubits):
         state = format(i, f"0{nqubits}b")
 
@@ -495,17 +495,21 @@ def calibration_matrix(nqubits, noise_model=None, nshots: int = 1000, backend=No
                 circuit.add(gates.X(q))
         circuit.add(gates.M(*range(nqubits)))
 
+        
         if noise_model is not None and backend.name != "qibolab":
             circuit = noise_model.apply(circuit)
-
+        circuit = transpile_circ(circuit, qubit_map, backend)
         freq = backend.execute_circuit(circuit, nshots=nshots).frequencies()
-
         column = np.zeros(2**nqubits)
         for key in freq.keys():
             f = freq[key] / nshots
             column[int(key, 2)] = f
         matrix[:, i] = column
 
+    # iqm 
+
+    # matrix[:, [0, 1]] = matrix[:, [1, 0]]
+
     return np.linalg.inv(matrix)
 
 
@@ -520,12 +524,12 @@ def apply_readout_mitigation(state, calibration_matrix):
     Returns:
         :class:`qibo.states.CircuitResult`: the input state with the updated frequencies.
     """
-    freq = np.zeros(2**state.nqubits)
+    from qibo.config import log
+    freq = np.zeros(2**len(state.measurements))
     for k, v in state.frequencies().items():
         freq[int(k, 2)] = v
 
     freq = freq.reshape(-1, 1)
-
     for i, val in enumerate(calibration_matrix @ freq):
         state._frequencies[i] = float(val)
 
@@ -685,8 +689,9 @@ def escircuit(circuit, obs, fuse=True, backend=None):
     if backend is None:  # pragma: no cover
         backend = GlobalBackend()
 
+    sign = -1
+    #while sign == -1:
     circ_cliff = sample_clifford_training_circuit(circuit, backend=backend)
-
     c_unitary = circ_cliff.unitary(backend=backend)
     nqubits = circ_cliff.nqubits
     U_c2p = comp_basis_to_pauli(nqubits, backend=backend)
@@ -695,6 +700,8 @@ def escircuit(circuit, obs, fuse=True, backend=None):
     )
     obs_pauli_liouville = U_c2p @ obs_liouville
     index = np.where(abs(obs_pauli_liouville) >= 1e-5)[0][0]
+    sign = np.sign(obs_pauli_liouville[index]) 
+
     obs1 = list(product(["I", "X", "Y", "Z"], repeat=nqubits))[index]
 
     paulis = {
@@ -743,6 +750,7 @@ def escircuit(circuit, obs, fuse=True, backend=None):
 def mit_obs(
     circuit,
     observable,
+    readout = {},
     qubit_map=None,
     noise_model=None,
     nshots=int(1e4),
@@ -762,7 +770,7 @@ def mit_obs(
     ]
 
     data = {"noise-free": {"-1": [], "1": []}, "noisy": {"-1": [], "1": []}}
-
+    from qibo.config import log
     a_list = {"-1": [], "1": []}
     for c in training_circs:
         circuit_result = c(nshots=nshots)
@@ -774,10 +782,17 @@ def mit_obs(
 
         c = transpile_circ(c, qubit_map, backend)
         circuit_result = backend.execute_circuit(c, nshots=nshots)
-        exp_noisy = observable.expectation_from_samples(circuit_result.frequencies()) #Add minus for iqm
+
+        if "calibration_matrix" in readout.keys():
+            circuit_result = apply_readout_mitigation(
+                circuit_result, readout["calibration_matrix"]
+            )
+
+        log.info(c.draw())
+        exp_noisy =   observable.expectation_from_samples(circuit_result.frequencies()) #Add minus for iqm
         # state = c_noisy().state()
         # exp_noisy = obs.expectation(state)
-
+        log.info(str(exp)+' '+str(exp_noisy))
         if exp > 0:
             data["noise-free"]["1"].append(exp)
             data["noisy"]["1"].append(exp_noisy)
@@ -795,12 +810,20 @@ def mit_obs(
         / n_training_samples
     )
 
+    # a = np.mean(a_list["1"])
+    # a_std = a_std_1
+
     circuit = circuit.fuse(max_qubits=1)
     if noise_model is not None and backend.name != "qibolab":
         circuit = noise_model.apply(circuit)
     circuit = transpile_circ(circuit, qubit_map, backend)
+    log.info(circuit.draw())
     circuit_result = backend.execute_circuit(circuit, nshots=nshots)
-    exp_noisy = observable.expectation_from_samples(circuit_result.frequencies()) # Add - for iqm
+    if "calibration_matrix" in readout.keys():
+        circuit_result = apply_readout_mitigation(
+            circuit_result, readout["calibration_matrix"]
+        )
+    exp_noisy =  observable.expectation_from_samples(circuit_result.frequencies()) # Add - for iqm
     # exp_noisy = obs.expectation(c_noisy().state())
     exp_mit = (1 - a) * exp_noisy / ((1 - a) ** 2 + a_std**2)
     exp_mit_std = (

From e24014652929d374c0c166f95f8f8a55371007d7 Mon Sep 17 00:00:00 2001
From: Alejandro <alexsopena15@gmail.com>
Date: Mon, 30 Oct 2023 17:00:16 +0400
Subject: [PATCH 12/47] bayesian unfolding readout

---
 src/qibo/models/error_mitigation.py | 40 +++++++++++++++++++++--------
 1 file changed, 30 insertions(+), 10 deletions(-)

diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index 84d5bf502c..39c3caf721 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -463,7 +463,14 @@ def vnCDR(
     return mit_val
 
 
-def calibration_matrix(nqubits, qubit_map, noise_model=None, nshots: int = 1000, backend=None):
+def ibu(prob, mat, iters=10):
+    t = np.ones((len(prob),1))/len(prob)
+    for k in range(iters):
+        t = t*(np.transpose(mat)@(prob/(mat@t)))
+
+    return t
+
+def calibration_matrix(nqubits, qubit_map, inv=False, noise_model=None, nshots: int = 1000, backend=None):
     """Computes the calibration matrix for readout mitigation.
 
     Args:
@@ -508,12 +515,16 @@ def calibration_matrix(nqubits, qubit_map, noise_model=None, nshots: int = 1000,
 
     # iqm 
 
-    # matrix[:, [0, 1]] = matrix[:, [1, 0]]
+    #matrix[[0, 1], :] = matrix[[1, 0], :]
+
 
-    return np.linalg.inv(matrix)
+    if inv == True:
+        matrix = np.linalg.inv(matrix)
 
+    return matrix
 
-def apply_readout_mitigation(state, calibration_matrix):
+
+def apply_readout_mitigation(state, calibration_matrix, inv=False):
     """Updates the frequencies of the input state with the mitigated ones obtained with
     ``calibration_matrix * state.frequencies()``.
 
@@ -530,8 +541,17 @@ def apply_readout_mitigation(state, calibration_matrix):
         freq[int(k, 2)] = v
 
     freq = freq.reshape(-1, 1)
-    for i, val in enumerate(calibration_matrix @ freq):
-        state._frequencies[i] = float(val)
+
+    if inv:
+        for i, val in enumerate(calibration_matrix @ freq):
+            state._frequencies[i] = float(val)
+    else:
+        prob_mit = ibu(freq/np.sum(freq),calibration_matrix)
+        freq_mit = np.round(prob_mit*np.sum(freq),0)
+        freq_mit = (freq_mit/np.sum(freq_mit))*np.sum(freq)
+
+        for i, val in enumerate(freq_mit):
+            state._frequencies[i] = float(freq_mit[i])
 
     return state
 
@@ -785,7 +805,7 @@ def mit_obs(
 
         if "calibration_matrix" in readout.keys():
             circuit_result = apply_readout_mitigation(
-                circuit_result, readout["calibration_matrix"]
+                circuit_result, readout["calibration_matrix"], readout["inv"]
             )
 
         log.info(c.draw())
@@ -821,9 +841,9 @@ def mit_obs(
     circuit_result = backend.execute_circuit(circuit, nshots=nshots)
     if "calibration_matrix" in readout.keys():
         circuit_result = apply_readout_mitigation(
-            circuit_result, readout["calibration_matrix"]
+            circuit_result, readout["calibration_matrix"], readout["inv"]
         )
-    exp_noisy =  observable.expectation_from_samples(circuit_result.frequencies()) # Add - for iqm
+    exp_noisy =   observable.expectation_from_samples(circuit_result.frequencies()) # Add - for iqm
     # exp_noisy = obs.expectation(c_noisy().state())
     exp_mit = (1 - a) * exp_noisy / ((1 - a) ** 2 + a_std**2)
     exp_mit_std = (
@@ -832,7 +852,7 @@ def mit_obs(
         * abs((1 - a) ** 2 - a_std**2)
         / ((1 - a) ** 2 + a_std**2) ** 2
     )
-
+    log.info(str(exp_mit)+' '+str(exp_noisy))
     if full_output:
         return exp_mit, exp_mit_std, a, a_std, a_list, data
 

From 81f44fe3dbf20ae53b85e77897ddff2f5a8d8c20 Mon Sep 17 00:00:00 2001
From: MatteoRobbiati <robbiati.matteo@gmail.com>
Date: Mon, 20 Nov 2023 16:38:06 +0100
Subject: [PATCH 13/47] remove logs

---
 src/qibo/models/error_mitigation.py | 45 ++++++++--------
 src/qibo/models/variational.py      | 79 +++++++----------------------
 2 files changed, 40 insertions(+), 84 deletions(-)

diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index 39c3caf721..b5fb1657b9 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -464,13 +464,16 @@ def vnCDR(
 
 
 def ibu(prob, mat, iters=10):
-    t = np.ones((len(prob),1))/len(prob)
+    t = np.ones((len(prob), 1)) / len(prob)
     for k in range(iters):
-        t = t*(np.transpose(mat)@(prob/(mat@t)))
+        t = t * (np.transpose(mat) @ (prob / (mat @ t)))
 
     return t
 
-def calibration_matrix(nqubits, qubit_map, inv=False, noise_model=None, nshots: int = 1000, backend=None):
+
+def calibration_matrix(
+    nqubits, qubit_map, inv=False, noise_model=None, nshots: int = 1000, backend=None
+):
     """Computes the calibration matrix for readout mitigation.
 
     Args:
@@ -493,6 +496,7 @@ def calibration_matrix(nqubits, qubit_map, inv=False, noise_model=None, nshots:
 
     matrix = np.zeros((2**nqubits, 2**nqubits))
     from qibo.config import log
+
     for i in range(2**nqubits):
         state = format(i, f"0{nqubits}b")
 
@@ -502,7 +506,6 @@ def calibration_matrix(nqubits, qubit_map, inv=False, noise_model=None, nshots:
                 circuit.add(gates.X(q))
         circuit.add(gates.M(*range(nqubits)))
 
-        
         if noise_model is not None and backend.name != "qibolab":
             circuit = noise_model.apply(circuit)
         circuit = transpile_circ(circuit, qubit_map, backend)
@@ -513,10 +516,9 @@ def calibration_matrix(nqubits, qubit_map, inv=False, noise_model=None, nshots:
             column[int(key, 2)] = f
         matrix[:, i] = column
 
-    # iqm 
-
-    #matrix[[0, 1], :] = matrix[[1, 0], :]
+    # iqm
 
+    # matrix[[0, 1], :] = matrix[[1, 0], :]
 
     if inv == True:
         matrix = np.linalg.inv(matrix)
@@ -536,7 +538,8 @@ def apply_readout_mitigation(state, calibration_matrix, inv=False):
         :class:`qibo.states.CircuitResult`: the input state with the updated frequencies.
     """
     from qibo.config import log
-    freq = np.zeros(2**len(state.measurements))
+
+    freq = np.zeros(2 ** len(state.measurements))
     for k, v in state.frequencies().items():
         freq[int(k, 2)] = v
 
@@ -546,9 +549,9 @@ def apply_readout_mitigation(state, calibration_matrix, inv=False):
         for i, val in enumerate(calibration_matrix @ freq):
             state._frequencies[i] = float(val)
     else:
-        prob_mit = ibu(freq/np.sum(freq),calibration_matrix)
-        freq_mit = np.round(prob_mit*np.sum(freq),0)
-        freq_mit = (freq_mit/np.sum(freq_mit))*np.sum(freq)
+        prob_mit = ibu(freq / np.sum(freq), calibration_matrix)
+        freq_mit = np.round(prob_mit * np.sum(freq), 0)
+        freq_mit = (freq_mit / np.sum(freq_mit)) * np.sum(freq)
 
         for i, val in enumerate(freq_mit):
             state._frequencies[i] = float(freq_mit[i])
@@ -710,7 +713,7 @@ def escircuit(circuit, obs, fuse=True, backend=None):
         backend = GlobalBackend()
 
     sign = -1
-    #while sign == -1:
+    # while sign == -1:
     circ_cliff = sample_clifford_training_circuit(circuit, backend=backend)
     c_unitary = circ_cliff.unitary(backend=backend)
     nqubits = circ_cliff.nqubits
@@ -720,7 +723,7 @@ def escircuit(circuit, obs, fuse=True, backend=None):
     )
     obs_pauli_liouville = U_c2p @ obs_liouville
     index = np.where(abs(obs_pauli_liouville) >= 1e-5)[0][0]
-    sign = np.sign(obs_pauli_liouville[index]) 
+    sign = np.sign(obs_pauli_liouville[index])
 
     obs1 = list(product(["I", "X", "Y", "Z"], repeat=nqubits))[index]
 
@@ -770,7 +773,7 @@ def escircuit(circuit, obs, fuse=True, backend=None):
 def mit_obs(
     circuit,
     observable,
-    readout = {},
+    readout={},
     qubit_map=None,
     noise_model=None,
     nshots=int(1e4),
@@ -791,6 +794,7 @@ def mit_obs(
 
     data = {"noise-free": {"-1": [], "1": []}, "noisy": {"-1": [], "1": []}}
     from qibo.config import log
+
     a_list = {"-1": [], "1": []}
     for c in training_circs:
         circuit_result = c(nshots=nshots)
@@ -808,11 +812,8 @@ def mit_obs(
                 circuit_result, readout["calibration_matrix"], readout["inv"]
             )
 
-        log.info(c.draw())
-        exp_noisy =   observable.expectation_from_samples(circuit_result.frequencies()) #Add minus for iqm
-        # state = c_noisy().state()
-        # exp_noisy = obs.expectation(state)
-        log.info(str(exp)+' '+str(exp_noisy))
+        exp_noisy = observable.expectation_from_samples(circuit_result.frequencies())
+
         if exp > 0:
             data["noise-free"]["1"].append(exp)
             data["noisy"]["1"].append(exp_noisy)
@@ -837,13 +838,14 @@ def mit_obs(
     if noise_model is not None and backend.name != "qibolab":
         circuit = noise_model.apply(circuit)
     circuit = transpile_circ(circuit, qubit_map, backend)
-    log.info(circuit.draw())
     circuit_result = backend.execute_circuit(circuit, nshots=nshots)
     if "calibration_matrix" in readout.keys():
         circuit_result = apply_readout_mitigation(
             circuit_result, readout["calibration_matrix"], readout["inv"]
         )
-    exp_noisy =   observable.expectation_from_samples(circuit_result.frequencies()) # Add - for iqm
+    exp_noisy = observable.expectation_from_samples(
+        circuit_result.frequencies()
+    )  # Add - for iqm
     # exp_noisy = obs.expectation(c_noisy().state())
     exp_mit = (1 - a) * exp_noisy / ((1 - a) ** 2 + a_std**2)
     exp_mit_std = (
@@ -852,7 +854,6 @@ def mit_obs(
         * abs((1 - a) ** 2 - a_std**2)
         / ((1 - a) ** 2 + a_std**2) ** 2
     )
-    log.info(str(exp_mit)+' '+str(exp_noisy))
     if full_output:
         return exp_mit, exp_mit_std, a, a_std, a_list, data
 
diff --git a/src/qibo/models/variational.py b/src/qibo/models/variational.py
index 7172bf07d8..c9cfefd938 100644
--- a/src/qibo/models/variational.py
+++ b/src/qibo/models/variational.py
@@ -36,45 +36,23 @@ def __init__(self, circuit, hamiltonian):
 
     def minimize(
         self,
-        initial_state,
-        method="Powell",
-        jac=None,
-        hess=None,
-        hessp=None,
-        bounds=None,
-        constraints=(),
-        tol=None,
-        callback=None,
-        options=None,
-        compile=False,
-        processes=None,
+        initial_parameters,
+        optimizer=None,
+        optimizer_options={},
     ):
         """Search for parameters which minimizes the hamiltonian expectation.
 
         Args:
-            initial_state (array): a initial guess for the parameters of the
+            initial_parameters (array or list): a initial guess for the parameters of the
                 variational circuit.
-            method (str): the desired minimization method.
-                See :meth:`qibo.optimizers.optimize` for available optimization
-                methods.
-            jac (dict): Method for computing the gradient vector for scipy optimizers.
-            hess (dict): Method for computing the hessian matrix for scipy optimizers.
-            hessp (callable): Hessian of objective function times an arbitrary
-                vector for scipy optimizers.
-            bounds (sequence or Bounds): Bounds on variables for scipy optimizers.
-            constraints (dict): Constraints definition for scipy optimizers.
-            tol (float): Tolerance of termination for scipy optimizers.
-            callback (callable): Called after each iteration for scipy optimizers.
-            options (dict): a dictionary with options for the different optimizers.
-            compile (bool): whether the TensorFlow graph should be compiled.
-            processes (int): number of processes when using the paralle BFGS method.
+            optimizer (qibo.optimizers.Optimizer): Qibo optimizer which is used to
+                minimize the energy [default is `ScipyMinimizer(method="Powell")`].
+            optimizer_options (dict): extra options for the optimizers.
 
         Return:
-            The final expectation value.
-            The corresponding best parameters.
-            The optimization result object. For scipy methods it returns
-            the ``OptimizeResult``, for ``'cma'`` the ``CMAEvolutionStrategy.result``,
-            and for ``'sgd'`` the options used during the optimization.
+            output of `qibo.optimizers.Optimizer` classes: [best function value,
+                best set of trainable parameters, complete result object according
+                to specific optimizer configuration].
         """
 
         def _loss(params, circuit, hamiltonian):
@@ -83,38 +61,15 @@ def _loss(params, circuit, hamiltonian):
             final_state = result.state()
             return hamiltonian.expectation(final_state)
 
-        if compile:
-            loss = self.hamiltonian.backend.compile(_loss)
-        else:
-            loss = _loss
-
-        if method == "cma":
-            # TODO: check if we can use this shortcut
-            # dtype = getattr(self.hamiltonian.backend.np, self.hamiltonian.backend._dtypes.get('DTYPE'))
-            dtype = self.hamiltonian.backend.np.float64
-            loss = lambda p, c, h: dtype(_loss(p, c, h))
-        elif method != "sgd":
-            loss = lambda p, c, h: self.hamiltonian.backend.to_numpy(_loss(p, c, h))
-
-        result, parameters, extra = self.optimizers.optimize(
-            loss,
-            initial_state,
+        opt = optimizer(
+            initial_parameters=initial_parameters,
+            loss=_loss,
             args=(self.circuit, self.hamiltonian),
-            method=method,
-            jac=jac,
-            hess=hess,
-            hessp=hessp,
-            bounds=bounds,
-            constraints=constraints,
-            tol=tol,
-            callback=callback,
-            options=options,
-            compile=compile,
-            processes=processes,
-            backend=self.hamiltonian.backend,
+            **optimizer_options,
         )
-        self.circuit.set_parameters(parameters)
-        return result, parameters, extra
+
+        best_f, best_p, results = opt.fit()
+        return best_f, best_p, results
 
 
 class AAVQE:

From 1aed5d4a1468ad2258fe5f07cef0025b77bb0e30 Mon Sep 17 00:00:00 2001
From: MatteoRobbiati <robbiati.matteo@gmail.com>
Date: Mon, 20 Nov 2023 17:14:36 +0100
Subject: [PATCH 14/47] rm comments and add preliminary docstrings

---
 src/qibo/models/error_mitigation.py | 28 +++++++++++-----------------
 1 file changed, 11 insertions(+), 17 deletions(-)

diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index b5fb1657b9..ceff367882 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -464,6 +464,7 @@ def vnCDR(
 
 
 def ibu(prob, mat, iters=10):
+    """Iterative Bayesian Unfolding."""
     t = np.ones((len(prob), 1)) / len(prob)
     for k in range(iters):
         t = t * (np.transpose(mat) @ (prob / (mat @ t)))
@@ -516,10 +517,6 @@ def calibration_matrix(
             column[int(key, 2)] = f
         matrix[:, i] = column
 
-    # iqm
-
-    # matrix[[0, 1], :] = matrix[[1, 0], :]
-
     if inv == True:
         matrix = np.linalg.inv(matrix)
 
@@ -707,6 +704,10 @@ def transpile_circ(circuit, qubit_map, backend):
 
 
 def escircuit(circuit, obs, fuse=True, backend=None):
+    """
+    Error sensitive circuit algorithm.
+    Implement what proposed in https://arxiv.org/abs/2112.06255.
+    """
     from qibo.quantum_info import comp_basis_to_pauli, random_clifford, vectorization
 
     if backend is None:  # pragma: no cover
@@ -750,11 +751,6 @@ def escircuit(circuit, obs, fuse=True, backend=None):
 
     circ_cliff1 = circ_cliff.__class__(**circ_cliff.init_kwargs)
 
-    # for gate in adjust_gates:
-    #     circ_cliff1.add(gate)
-
-    # for gate in circ_cliff.queue:
-    #     circ_cliff1.add(gate)
     j = 0
     for gate in circ_cliff.queue:
         if gate.name == "id_init":
@@ -770,7 +766,7 @@ def escircuit(circuit, obs, fuse=True, backend=None):
     return circ_cliff1, circ_cliff, adjust_gates
 
 
-def mit_obs(
+def ICS(
     circuit,
     observable,
     readout={},
@@ -781,6 +777,10 @@ def mit_obs(
     full_output=False,
     backend=None,
 ):
+    """
+    Compute the important Clifford Sampling methos proposed in
+    https://arxiv.org/abs/2112.06255.
+    """
     if backend is None:  # pragma: no cover
         backend = GlobalBackend()
 
@@ -831,9 +831,6 @@ def mit_obs(
         / n_training_samples
     )
 
-    # a = np.mean(a_list["1"])
-    # a_std = a_std_1
-
     circuit = circuit.fuse(max_qubits=1)
     if noise_model is not None and backend.name != "qibolab":
         circuit = noise_model.apply(circuit)
@@ -843,10 +840,7 @@ def mit_obs(
         circuit_result = apply_readout_mitigation(
             circuit_result, readout["calibration_matrix"], readout["inv"]
         )
-    exp_noisy = observable.expectation_from_samples(
-        circuit_result.frequencies()
-    )  # Add - for iqm
-    # exp_noisy = obs.expectation(c_noisy().state())
+    exp_noisy = observable.expectation_from_samples(circuit_result.frequencies())
     exp_mit = (1 - a) * exp_noisy / ((1 - a) ** 2 + a_std**2)
     exp_mit_std = (
         a_std

From c238ad9ea44ffcf3eefef2bbb82238bc358092d8 Mon Sep 17 00:00:00 2001
From: Alejandro Sopena <alexsopena15@gmail.com>
Date: Mon, 11 Dec 2023 01:46:31 +0100
Subject: [PATCH 15/47] get_gammas

---
 src/qibo/models/__init__.py         |  2 +-
 src/qibo/models/error_mitigation.py | 37 ++++++++++++++++-------------
 2 files changed, 21 insertions(+), 18 deletions(-)

diff --git a/src/qibo/models/__init__.py b/src/qibo/models/__init__.py
index 0adff068c3..f0d07da334 100644
--- a/src/qibo/models/__init__.py
+++ b/src/qibo/models/__init__.py
@@ -5,7 +5,7 @@
     ZNE,
     get_gammas,
     get_noisy_circuit,
-    sample_training_circuit,
+    sample_training_circuit_cdr,
     vnCDR,
 )
 from qibo.models.evolution import AdiabaticEvolution, StateEvolution
diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index ceff367882..832dca5a25 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -10,42 +10,45 @@
 from qibo.config import raise_error
 
 
-def get_gammas(c, solve: bool = True):
+def get_gammas(noise_levels, analytical: bool = True):
     """Standalone function to compute the ZNE coefficients given the noise levels.
 
     Args:
-        c (numpy.ndarray): array containing the different noise levels.
+        noise_levels (numpy.ndarray): array containing the different noise levels.
             Note that in the CNOT insertion paradigm this corresponds to
             the number of CNOT pairs to be inserted. The canonical ZNE
             noise levels are obtained as ``2 * c + 1``.
-        solve (bool, optional): If ``True``, computes the coeffients by solving the
+        analytical (bool, optional): if ``True``, computes the coeffients by solving the
             linear system. If ``False``, use the analytical solution valid
             for the CNOT insertion method. Default is ``True``.
 
     Returns:
-        numpy.ndarray: The computed coefficients.
+        numpy.ndarray: the computed coefficients.
     """
-    if solve:
-        c = 2 * c + 1
-        a = np.array([c**i for i in range(len(c))])
-        b = np.zeros(len(c))
-        b[0] = 1
-        gammas = np.linalg.solve(a, b)
+    if analytical:
+        noise_levels = 2 * noise_levels + 1
+        a_matrix = np.array([noise_levels**i for i in range(len(noise_levels))])
+        b_vector = np.zeros(len(noise_levels))
+        b_vector[0] = 1
+        zne_coefficients = np.linalg.solve(a_matrix, b_vector)
     else:
-        cmax = c[-1]
-        gammas = np.array(
+        max_noise_level = noise_levels[-1]
+        zne_coefficients = np.array(
             [
                 1
-                / (2 ** (2 * cmax) * np.math.factorial(i))
+                / (2 ** (2 * max_noise_level) * np.math.factorial(i))
                 * (-1) ** i
                 / (1 + 2 * i)
-                * np.math.factorial(1 + 2 * cmax)
-                / (np.math.factorial(cmax) * np.math.factorial(cmax - i))
-                for i in c
+                * np.math.factorial(1 + 2 * max_noise_level)
+                / (
+                    np.math.factorial(max_noise_level)
+                    * np.math.factorial(max_noise_level - i)
+                )
+                for i in noise_levels
             ]
         )
 
-    return gammas
+    return zne_coefficients
 
 
 def get_noisy_circuit(circuit, num_insertions: int, insertion_gate: str = "CNOT"):

From b57184ef32ede094ccd45a7b907671dbc0329596 Mon Sep 17 00:00:00 2001
From: Alejandro Sopena <alexsopena15@gmail.com>
Date: Mon, 11 Dec 2023 01:47:46 +0100
Subject: [PATCH 16/47] get_noisy_circuit

---
 src/qibo/models/error_mitigation.py | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index 832dca5a25..7f5ca9f1e2 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -62,7 +62,7 @@ def get_noisy_circuit(circuit, num_insertions: int, insertion_gate: str = "CNOT"
             Default is ``"CNOT"``.
 
     Returns:
-        :class:`qibo.models.Circuit`: The circuit with the inserted CNOT pairs.
+        :class:`qibo.models.Circuit`: circuit with the inserted gate pairs.
     """
     if insertion_gate not in ("CNOT", "RX"):  # pragma: no cover
         raise_error(
@@ -83,16 +83,17 @@ def get_noisy_circuit(circuit, num_insertions: int, insertion_gate: str = "CNOT"
 
     for gate in circuit.queue:
         noisy_circuit.add(gate)
+
         if isinstance(gate, i_gate):
             if insertion_gate == "CNOT":
                 control = gate.control_qubits[0]
                 target = gate.target_qubits[0]
-                for i in range(num_insertions):
+                for _ in range(num_insertions):
                     noisy_circuit.add(gates.CNOT(control, target))
                     noisy_circuit.add(gates.CNOT(control, target))
             elif gate.init_kwargs["theta"] == theta:
                 qubit = gate.qubits[0]
-                for i in range(num_insertions):
+                for _ in range(num_insertions):
                     noisy_circuit.add(gates.RX(qubit, theta=theta))
                     noisy_circuit.add(gates.RX(qubit, theta=-theta))
 

From 7fc65603e864eb1c7c25e08936936b5fbb892c09 Mon Sep 17 00:00:00 2001
From: Alejandro Sopena <alexsopena15@gmail.com>
Date: Mon, 11 Dec 2023 01:48:30 +0100
Subject: [PATCH 17/47] zero noise extrapolation

---
 src/qibo/models/error_mitigation.py | 36 +++++++++++++++++------------
 1 file changed, 21 insertions(+), 15 deletions(-)

diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index 7f5ca9f1e2..c899ac0745 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -123,26 +123,32 @@ def ZNE(
         noise_levels (numpy.ndarray): Sequence of noise levels.
         noise_model (:class:`qibo.noise.NoiseModel`, optional): Noise model applied
             to simulate noisy computation.
-        nshots (int, optional): Number of shots.
+        nshots (int, optional): Number of shots. Defauylts to 10000.
         solve_for_gammas (bool, optional): If ``True``, explicitly solve the
-            equations to obtain the ``gamma`` coefficients.
+            equations to obtain the ``gamma`` coefficients. Default is ``False``.
         insertion_gate (str, optional): gate to be used in the insertion.
             If ``"RX"``, the gate used is :math:``RX(\\pi / 2)``.
-            Default is ``"CNOT"``.
-        readout (dict, optional): It has the structure
-            {'calibration_matrix': `numpy.ndarray`, 'ncircuits': `int`}.
-            If passed, the calibration matrix or the randomized method is
-            used to mitigate readout errors.
-        backend (:class:`qibo.backends.abstract.Backend`, optional): Calculation engine.
+            Defaults to ``"CNOT"``.
+        readout (dict, optional): A dictionary that may contain the following keys:
+            - 'calibration_matrix': numpy.ndarray, used for applying a pre-computed calibration matrix for readout error mitigation.
+            - 'random_ncircuits': int, specifies the number of random circuits to use for the randomized method of readout error mitigation.
+            - 'ibu_iters': int, specifies the number of iterations for the iterative Bayesian update method of readout error mitigation.
+            If provided, the corresponding readout error mitigation method is used.
+        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``.
 
     Returns:
         numpy.ndarray: Estimate of the expected value of ``observable`` in the noise free condition.
-    """
 
+    Reference:
+        1. K. Temme, S. Bravyi et al, *Error mitigation for short-depth quantum circuits*.
+           `arXiv:1612.02058 [quant-ph] <https://arxiv.org/abs/1612.02058>`_.
+    """
     if backend is None:  # pragma: no cover
         backend = GlobalBackend()
 
-    expected_val = []
+    expected_values = []
     for num_insertions in noise_levels:
         noisy_circuit = get_noisy_circuit(
             circuit, num_insertions, insertion_gate=insertion_gate
@@ -156,17 +162,17 @@ def ZNE(
                 noisy_circuit = noise_model.apply(noisy_circuit)
             circuit_result = backend.execute_circuit(noisy_circuit, nshots=nshots)
         if "calibration_matrix" in readout.keys() is not None:
-            circuit_result = apply_readout_mitigation(
-                circuit_result, readout["calibration_matrix"]
+            circuit_result = apply_cal_mat_readout_mitigation(
+                circuit_result, readout["calibration_matrix"], readout["ibu_iters"]
             )
         val = circuit_result.expectation_from_samples(observable)
         if "ncircuits" in readout.keys():
             val /= circuit_result_cal.expectation_from_samples(observable)
-        expected_val.append(val)
+        expected_values.append(val)
 
-    gamma = get_gammas(noise_levels, solve=solve_for_gammas)
+    gamma = get_gammas(noise_levels, analytical=solve_for_gammas)
 
-    return np.sum(gamma * expected_val)
+    return np.sum(gamma * expected_values)
 
 
 def sample_training_circuit(

From cf792bffd172d625adb3d878db850f42866a7986 Mon Sep 17 00:00:00 2001
From: Alejandro Sopena <alexsopena15@gmail.com>
Date: Mon, 11 Dec 2023 01:49:46 +0100
Subject: [PATCH 18/47] clifford data regression

---
 src/qibo/models/error_mitigation.py | 84 ++++++++++++++---------------
 1 file changed, 42 insertions(+), 42 deletions(-)

diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index c899ac0745..65ab5bf1dc 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -175,7 +175,7 @@ def ZNE(
     return np.sum(gamma * expected_values)
 
 
-def sample_training_circuit(
+def sample_training_circuit_cdr(
     circuit,
     replacement_gates: list = None,
     sigma: float = 0.5,
@@ -191,7 +191,9 @@ def sample_training_circuit(
             form (``gates.XYZ``, ``kwargs``). For example, phase gates are used by default:
             ``list((RZ, {'theta':0}), (RZ, {'theta':pi/2}), (RZ, {'theta':pi}), (RZ, {'theta':3*pi/2}))``.
         sigma (float, optional): standard devation of the Gaussian distribution used for sampling.
-        backend (:class:`qibo.backends.abstract.Backend`, optional): Calculation engine.
+        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``.
 
     Returns:
         :class:`qibo.models.Circuit`: The sampled circuit.
@@ -202,18 +204,15 @@ def sample_training_circuit(
     if replacement_gates is None:
         replacement_gates = [(gates.RZ, {"theta": n * np.pi / 2}) for n in range(4)]
 
-    # Find all the non-Clifford RZ gates
     gates_to_replace = []
     for i, gate in enumerate(circuit.queue):
         if isinstance(gate, gates.RZ):
             if gate.init_kwargs["theta"] % (np.pi / 2) != 0.0:
                 gates_to_replace.append((i, gate))
 
-    if len(gates_to_replace) == 0:
+    if not gates_to_replace:
         raise_error(ValueError, "No non-Clifford RZ gate found, no circuit sampled.")
 
-    # For each RZ gate build the possible candidates and
-    # compute the frobenius distance to the candidates
     replacement, distance = [], []
     for _, gate in gates_to_replace:
         rep_gates = np.array(
@@ -231,31 +230,28 @@ def sample_training_circuit(
         )
 
     distance = np.vstack(distance)
-    # Compute the scores
     prob = np.exp(-(distance**2) / sigma**2)
-    # Sample which of the RZ found to substitute
+
     index = np.random.choice(
         range(len(gates_to_replace)),
         size=min(int(len(gates_to_replace) / 2), 50),
         replace=False,
         p=prob.sum(-1) / prob.sum(),
     )
+
     gates_to_replace = np.array([gates_to_replace[i] for i in index])
     prob = [prob[i] for i in index]
-    # Sample which replacement gate to substitute with
+
     replacement = np.array([replacement[i] for i in index])
     replacement = [
         replacement[i][np.random.choice(range(len(p)), size=1, p=p / p.sum())[0]]
         for i, p in enumerate(prob)
     ]
     replacement = {i[0]: g for i, g in zip(gates_to_replace, replacement)}
-    # Build the training circuit by substituting the sampled gates
+
     sampled_circuit = circuit.__class__(**circuit.init_kwargs)
     for i, gate in enumerate(circuit.queue):
-        if i in replacement.keys():
-            sampled_circuit.add(replacement[i])
-        else:
-            sampled_circuit.add(gate)
+        sampled_circuit.add(replacement.get(i, gate))
 
     return sampled_circuit
 
@@ -279,57 +275,62 @@ def CDR(
         observable (numpy.ndarray): observable to be measured.
         noise_model (:class:`qibo.noise.NoiseModel`): noise model used for simulating
             noisy computation.
-        nshots (int, optional): number of shots.
+        nshots (int, optional): number of shots. Defaults 10000.
         model (callable, optional): model used for fitting. This should be a callable
             function object ``f(x, *params)``, taking as input the predictor variable
             and the parameters. Default is a simple linear model ``f(x,a,b) := a*x + b``.
-        n_training_samples (int, optional): number of training circuits to sample.
+        n_training_samples (int, optional): number of training circuits to sample. Defaults to 100.
         full_output (bool, optional): if ``True``, this function returns additional
-            information: ``val``, ``optimal_params``, ``train_val``.
-        readout (dict, optional): It has the structure
-            {'calibration_matrix': `numpy.ndarray`, 'ncircuits': `int`}.
-            If passed, the calibration matrix or the randomized method is
-            used to mitigate readout errors.
-        backend (:class:`qibo.backends.abstract.Backend`, optional): calculation engine.
+            information: ``val``, ``optimal_params``, ``train_val``. Defaults to ``False``.
+        readout (dict, optional): A dictionary that may contain the following keys:
+            - 'calibration_matrix': numpy.ndarray, used for applying a pre-computed calibration matrix for readout error mitigation.
+            - 'random_ncircuits': int, specifies the number of random circuits to use for the randomized method of readout error mitigation.
+            - 'ibu_iters': int, specifies the number of iterations for the iterative Bayesian update method of readout error mitigation.
+            If provided, the corresponding readout error mitigation method is used.
+        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``.
 
     Returns:
         mit_val (float): Mitigated expectation value of `observable`.
         val (float): Noisy expectation value of `observable`.
         optimal_params (list): Optimal values for `params`.
         train_val (dict): Contains the noise-free and noisy expectation values obtained with the training circuits.
-    """
 
-    # Set backend
+    Reference:
+        1. P. Czarnik, A. Arrasmith et al, *Error mitigation with Clifford quantum-circuit data*.
+           `arXiv:2005.10189 [quant-ph] <https://arxiv.org/abs/2005.10189>`_.
+    """
     if backend is None:  # pragma: no cover
         backend = GlobalBackend()
-    # Sample the training set
+
     training_circuits = [
-        sample_training_circuit(circuit) for n in range(n_training_samples)
+        sample_training_circuit_cdr(circuit) for _ in range(n_training_samples)
     ]
-    # Run the sampled circuits
+
     train_val = {"noise-free": [], "noisy": []}
-    for c in training_circuits:
-        val = c(nshots=nshots).expectation_from_samples(observable)
+    for circ in training_circuits:
+        val = circ(nshots=nshots).expectation_from_samples(observable)
         train_val["noise-free"].append(val)
         if "ncircuits" in readout.keys():
             circuit_result, circuit_result_cal = apply_randomized_readout_mitigation(
-                c, noise_model, nshots, readout["ncircuits"], backend
+                circ, noise_model, nshots, readout["ncircuits"], backend
             )
         else:
             if noise_model is not None and backend.name != "qibolab":
-                c = noise_model.apply(c)
-            circuit_result = backend.execute_circuit(c, nshots=nshots)
+                circ = noise_model.apply(circ)
+            circuit_result = backend.execute_circuit(circ, nshots=nshots)
         if "calibration_matrix" in readout.keys() is not None:
-            circuit_result = apply_readout_mitigation(
-                circuit_result, readout["calibration_matrix"]
+            circuit_result = apply_cal_mat_readout_mitigation(
+                circuit_result, readout["calibration_matrix"], readout["ibu_iters"]
             )
         val = circuit_result.expectation_from_samples(observable)
         if "ncircuits" in readout.keys():
             val /= circuit_result_cal.expectation_from_samples(observable)
         train_val["noisy"].append(val)
-    # Fit the model
+
     optimal_params = curve_fit(model, train_val["noisy"], train_val["noise-free"])[0]
-    # Run the input circuit
+
     if "ncircuits" in readout.keys():
         circuit_result, circuit_result_cal = apply_randomized_readout_mitigation(
             circuit, noise_model, nshots, readout["ncircuits"], backend
@@ -339,19 +340,18 @@ def CDR(
             circuit = noise_model.apply(circuit)
         circuit_result = backend.execute_circuit(circuit, nshots=nshots)
     if "calibration_matrix" in readout.keys() is not None:
-        circuit_result = apply_readout_mitigation(
-            circuit_result, readout["calibration_matrix"]
+        circuit_result = apply_cal_mat_readout_mitigation(
+            circuit_result, readout["calibration_matrix"], readout["ibu_iters"]
         )
     val = circuit_result.expectation_from_samples(observable)
     if "ncircuits" in readout.keys():
         val /= circuit_result_cal.expectation_from_samples(observable)
     mit_val = model(val, *optimal_params)
 
-    # Return data
-    if full_output == True:
+    if full_output is True:
         return mit_val, val, optimal_params, train_val
-    else:
-        return mit_val
+
+    return mit_val
 
 
 def vnCDR(

From d4414acf1e136ce6bbc905777c6a476f7649f9b4 Mon Sep 17 00:00:00 2001
From: Alejandro Sopena <alexsopena15@gmail.com>
Date: Mon, 11 Dec 2023 01:50:12 +0100
Subject: [PATCH 19/47] vnCDR

---
 src/qibo/models/error_mitigation.py | 45 ++++++++++++++---------------
 1 file changed, 22 insertions(+), 23 deletions(-)

diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index 65ab5bf1dc..eb533a00aa 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -385,12 +385,15 @@ def vnCDR(
             If ``"RX"``, the gate used is :math:``RX(\\pi / 2)``.
             Default is ``"CNOT"``.
         full_output (bool, optional): if ``True``, this function returns additional
-            information: ``val``, ``optimal_params``, ``train_val``.
-        readout (dict, optional): It has the structure
-            {'calibration_matrix': `numpy.ndarray`, 'ncircuits': `int`}.
-            If passed, the calibration matrix or the randomized method is
-            used to mitigate readout errors.
-        backend (:class:`qibo.backends.abstract.Backend`, optional): calculation engine.
+            information: ``val``, ``optimal_params``, ``train_val``. Defaults to ``False``.
+        readout (dict, optional): A dictionary that may contain the following keys:
+            - 'calibration_matrix': numpy.ndarray, used for applying a pre-computed calibration matrix for readout error mitigation.
+            - 'random_ncircuits': int, specifies the number of random circuits to use for the randomized method of readout error mitigation.
+            - 'ibu_iters': int, specifies the number of iterations for the iterative Bayesian update method of readout error mitigation.
+            If provided, the corresponding readout error mitigation method is used.
+        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``.
 
     Returns:
         mit_val (float): Mitigated expectation value of `observable`.
@@ -398,24 +401,24 @@ def vnCDR(
         optimal_params (list): Optimal values for `params`.
         train_val (dict): Contains the noise-free and noisy expectation values obtained
         with the training circuits.
-    """
 
-    # Set backend
+    Reference:
+        1. A. Lowe, MH. Gordon et al, *Unified approach to data-driven quantum error mitigation*.
+           `arXiv:2011.01157 [quant-ph] <https://arxiv.org/abs/2011.01157>`_.
+    """
     if backend is None:  # pragma: no cover
         backend = GlobalBackend()
 
-    # Sample the training circuits
     training_circuits = [
-        sample_training_circuit(circuit) for n in range(n_training_samples)
+        sample_training_circuit_cdr(circuit) for _ in range(n_training_samples)
     ]
     train_val = {"noise-free": [], "noisy": []}
 
-    # Add the different noise levels and run the circuits
-    for c in training_circuits:
-        val = c(nshots=nshots).expectation_from_samples(observable)
+    for circ in training_circuits:
+        val = circ(nshots=nshots).expectation_from_samples(observable)
         train_val["noise-free"].append(val)
         for level in noise_levels:
-            noisy_c = get_noisy_circuit(c, level, insertion_gate=insertion_gate)
+            noisy_c = get_noisy_circuit(circ, level, insertion_gate=insertion_gate)
             if "ncircuits" in readout.keys():
                 (
                     circuit_result,
@@ -428,22 +431,19 @@ def vnCDR(
                     noisy_c = noise_model.apply(noisy_c)
                 circuit_result = backend.execute_circuit(noisy_c, nshots=nshots)
             if "calibration_matrix" in readout.keys():
-                circuit_result = apply_readout_mitigation(
-                    circuit_result, readout["calibration_matrix"]
+                circuit_result = apply_cal_mat_readout_mitigation(
+                    circuit_result, readout["calibration_matrix"], readout["ibu_iters"]
                 )
             val = circuit_result.expectation_from_samples(observable)
             if "ncircuits" in readout.keys():
                 val /= circuit_result_cal.expectation_from_samples(observable)
             train_val["noisy"].append(val)
 
-    # Repeat noise-free values for each noise level
     noisy_array = np.array(train_val["noisy"]).reshape(-1, len(noise_levels))
 
-    # Fit the model
     params = np.random.rand(len(noise_levels))
     optimal_params = curve_fit(model, noisy_array.T, train_val["noise-free"], p0=params)
 
-    # Run the input circuit
     val = []
     for level in noise_levels:
         noisy_c = get_noisy_circuit(circuit, level, insertion_gate=insertion_gate)
@@ -456,8 +456,8 @@ def vnCDR(
                 noisy_c = noise_model.apply(noisy_c)
             circuit_result = backend.execute_circuit(noisy_c, nshots=nshots)
         if "calibration_matrix" in readout.keys():
-            circuit_result = apply_readout_mitigation(
-                circuit_result, readout["calibration_matrix"]
+            circuit_result = apply_cal_mat_readout_mitigation(
+                circuit_result, readout["calibration_matrix"], readout["ibu_iters"]
             )
         expval = circuit_result.expectation_from_samples(observable)
         if "ncircuits" in readout.keys():
@@ -466,8 +466,7 @@ def vnCDR(
 
     mit_val = model(np.array(val).reshape(-1, 1), *optimal_params[0])[0]
 
-    # Return data
-    if full_output == True:
+    if full_output is True:
         return mit_val, val, optimal_params, train_val
 
     return mit_val

From 72e5ef8e4aaef95837e1a742a7fbd3403774e028 Mon Sep 17 00:00:00 2001
From: Alejandro Sopena <alexsopena15@gmail.com>
Date: Mon, 11 Dec 2023 01:51:57 +0100
Subject: [PATCH 20/47] readout mitigation

---
 src/qibo/models/error_mitigation.py | 148 +++++++++++++++++-----------
 1 file changed, 90 insertions(+), 58 deletions(-)

diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index eb533a00aa..a202d8ef95 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -472,17 +472,35 @@ def vnCDR(
     return mit_val
 
 
-def ibu(prob, mat, iters=10):
-    """Iterative Bayesian Unfolding."""
-    t = np.ones((len(prob), 1)) / len(prob)
-    for k in range(iters):
-        t = t * (np.transpose(mat) @ (prob / (mat @ t)))
+def iterative_bayesian_unfolding(probabilities, response_matrix, iterations=10):
+    """
+    Iterative Bayesian Unfolding (IBU) method for readout mitigation.
+
+    Args:
+        probabilities (numpy.ndarray): the input probabilities to be unfolded.
+        response_matrix (numpy.ndarray): the response matrix.
+        iterations (int, optional): the number of iterations to perform. Defaults to 10.
+
+    Returns:
+        numpy.ndarray: the unfolded probabilities.
 
-    return t
+    Reference:
+        1. S. Srinivasan, B. Pokharel et al, *Scalable Measurement Error Mitigation via Iterative Bayesian Unfolding*.
+           `arXiv:2210.12284 [quant-ph] <https://arxiv.org/abs/2210.12284>`_.
+    """
+    unfolded_probabilities = np.ones((len(probabilities), 1)) / len(probabilities)
 
+    for _ in range(iterations):
+        unfolded_probabilities = unfolded_probabilities * (
+            np.transpose(response_matrix)
+            @ (probabilities / (response_matrix @ unfolded_probabilities))
+        )
 
-def calibration_matrix(
-    nqubits, qubit_map, inv=False, noise_model=None, nshots: int = 1000, backend=None
+    return unfolded_probabilities
+
+
+def get_calibration_matrix(
+    nqubits, qubit_map, noise_model=None, nshots: int = 10000, backend=None
 ):
     """Computes the calibration matrix for readout mitigation.
 
@@ -491,76 +509,82 @@ def calibration_matrix(
         noise_model (:class:`qibo.noise.NoiseModel`, optional): noise model used for simulating
             noisy computation. This matrix can be used to mitigate the effect of
             `qibo.noise.ReadoutError`.
-        nshots (int, optional): number of shots.
-        backend (:class:`qibo.backends.abstract.Backend`, optional): calculation engine.
+        nshots (int, optional): number of shots. Defaults to 10000.
+        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``.
 
     Returns:
         numpy.ndarray : The computed (`nqubits`, `nqubits`) calibration matrix for
             readout mitigation.
     """
-
     from qibo import Circuit  # pylint: disable=import-outside-toplevel
 
     if backend is None:  # pragma: no cover
         backend = GlobalBackend()
 
-    matrix = np.zeros((2**nqubits, 2**nqubits))
-    from qibo.config import log
+    calibration_matrix = np.zeros((2**nqubits, 2**nqubits))
 
     for i in range(2**nqubits):
-        state = format(i, f"0{nqubits}b")
+        binary_state = format(i, f"0{nqubits}b")
 
         circuit = Circuit(nqubits, density_matrix=True)
-        for q, bit in enumerate(state):
+        for qubit, bit in enumerate(binary_state):
             if bit == "1":
-                circuit.add(gates.X(q))
+                circuit.add(gates.X(qubit))
         circuit.add(gates.M(*range(nqubits)))
 
         if noise_model is not None and backend.name != "qibolab":
             circuit = noise_model.apply(circuit)
         circuit = transpile_circ(circuit, qubit_map, backend)
-        freq = backend.execute_circuit(circuit, nshots=nshots).frequencies()
-        column = np.zeros(2**nqubits)
-        for key in freq.keys():
-            f = freq[key] / nshots
-            column[int(key, 2)] = f
-        matrix[:, i] = column
 
-    if inv == True:
-        matrix = np.linalg.inv(matrix)
+        frequencies = backend.execute_circuit(circuit, nshots=nshots).frequencies()
+
+        column = np.zeros(2**nqubits)
+        for key, value in frequencies.items():
+            column[int(key, 2)] = value / nshots
+        calibration_matrix[:, i] = column
 
-    return matrix
+    return calibration_matrix
 
 
-def apply_readout_mitigation(state, calibration_matrix, inv=False):
-    """Updates the frequencies of the input state with the mitigated ones obtained with
-    ``calibration_matrix * state.frequencies()``.
+def apply_cal_mat_readout_mitigation(state, calibration_matrix, iterations=None):
+    """
+    Applies readout error mitigation to the given state using the provided calibration matrix.
 
     Args:
-        state (:class:`qibo.states.CircuitResult`): input state to be updated.
-        calibration_matrix (numpy.ndarray): calibration matrix for readout mitigation.
+        state (:class:`qibo.states.CircuitResult`): the input state to be updated. This state should contain the
+            frequencies that need to be mitigated.
+        calibration_matrix (numpy.ndarray, optional): the calibration matrix for readout mitigation.
+        iterations (int, optional): the number of iterations to use for the Iterative Bayesian Unfolding method.
+            If ``None`` the 'inverse' method is used. Defaults to ``None``.
 
     Returns:
-        :class:`qibo.states.CircuitResult`: the input state with the updated frequencies.
+        :class:`qibo.states.CircuitResult`: The input state with the updated (mitigated) frequencies.
     """
-    from qibo.config import log
+    frequencies = np.zeros(2 ** len(state.measurements[0].qubits))
+    for key, value in state.frequencies().items():
+        frequencies[int(key, 2)] = value
 
-    freq = np.zeros(2 ** len(state.measurements))
-    for k, v in state.frequencies().items():
-        freq[int(k, 2)] = v
+    frequencies = frequencies.reshape(-1, 1)
 
-    freq = freq.reshape(-1, 1)
-
-    if inv:
-        for i, val in enumerate(calibration_matrix @ freq):
-            state._frequencies[i] = float(val)
+    if iterations is None:
+        calibration_matrix = np.linalg.inv(calibration_matrix)
+        for i, value in enumerate(calibration_matrix @ frequencies):
+            state._frequencies[i] = float(value)
     else:
-        prob_mit = ibu(freq / np.sum(freq), calibration_matrix)
-        freq_mit = np.round(prob_mit * np.sum(freq), 0)
-        freq_mit = (freq_mit / np.sum(freq_mit)) * np.sum(freq)
+        mitigated_probabilities = iterative_bayesian_unfolding(
+            frequencies / np.sum(frequencies), calibration_matrix, iterations
+        )
+        mitigated_frequencies = np.round(
+            mitigated_probabilities * np.sum(frequencies), 0
+        )
+        mitigated_frequencies = (
+            mitigated_frequencies / np.sum(mitigated_frequencies)
+        ) * np.sum(frequencies)
 
-        for i, val in enumerate(freq_mit):
-            state._frequencies[i] = float(freq_mit[i])
+        for i, value in enumerate(mitigated_frequencies):
+            state._frequencies[i] = float(value)
 
     return state
 
@@ -568,22 +592,27 @@ def apply_readout_mitigation(state, calibration_matrix, inv=False):
 def apply_randomized_readout_mitigation(
     circuit, noise_model=None, nshots: int = int(1e3), ncircuits: int = 10, backend=None
 ):
-    """Implements the readout mitigation method proposed in https://arxiv.org/abs/2012.09738.
+    """Readout mitigation method that transforms the bias in an expectation value into a measurable multiplicative factor. This factor can be eliminated at the expense of increased sampling complexity for the observable.
 
     Args:
         circuit (:class:`qibo.models.Circuit`): input circuit.
         noise_model(:class:`qibo.noise.NoiseModel`, optional): noise model used for
-            simulating noisy computation. This matrix can be used to mitigate the
-            effects of :class:`qibo.noise.ReadoutError`.
-        nshots (int, optional): number of shots.
+            simulating noisy computation. Defaults to ``None``.
+        nshots (int, optional): number of shots. Defaults to 10000.
         ncircuits (int, optional): number of randomized circuits. Each of them uses
-            ``int(nshots / ncircuits)`` shots.
-        backend (:class:`qibo.backends.abstract.Backend`): calculation engine.
+            ``int(nshots / ncircuits)`` shots. Defaults to 10.
+        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``.
 
     Return:
         :class:`qibo.states.CircuitResult`: the state of the input circuit with
             mitigated frequencies.
 
+
+    Reference:
+        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>`_.
     """
     from qibo import Circuit  # pylint: disable=import-outside-toplevel
     from qibo.quantum_info import (  # pylint: disable=import-outside-toplevel
@@ -593,7 +622,7 @@ def apply_randomized_readout_mitigation(
     if backend is None:  # pragma: no cover
         backend = GlobalBackend()
 
-    qubits = circuit.queue[-1].qubits
+    meas_qubits = circuit.measurements[0].qubits
     nshots_r = int(nshots / ncircuits)
     freq = np.zeros((ncircuits, 2), object)
     for k in range(ncircuits):
@@ -601,19 +630,22 @@ def apply_randomized_readout_mitigation(
         circuit_c.queue.pop()
         cal_circuit = Circuit(circuit.nqubits, density_matrix=True)
 
-        x_gate = random_pauli(len(qubits), 1, subset=["I", "X"]).queue
+        x_gate = random_pauli(circuit.nqubits, 1, subset=["I", "X"]).queue
 
         error_map = {}
-        for gate in x_gate:
+        for j, gate in enumerate(x_gate):
             if gate.name == "x":
-                error_map[gate.qubits[0]] = 1
+                if gate.qubits[0] in meas_qubits:
+                    error_map[gate.qubits[0]] = 1
+                else:
+                    x_gate.queue[j] = gates.I(gate.qubits[0])
 
         circuits = [circuit_c, cal_circuit]
         results = []
         freqs = []
         for circ in circuits:
             circ.add(x_gate)
-            circ.add(gates.M(*qubits))
+            circ.add(gates.M(*meas_qubits))
             if noise_model is not None and backend.name != "qibolab":
                 circ = noise_model.apply(circ)
             result = backend.execute_circuit(circ, nshots=nshots_r)
@@ -625,8 +657,8 @@ def apply_randomized_readout_mitigation(
     for j in range(2):
         results[j].nshots = nshots
         freq_sum = freq[0, j]
-        for f in freq[1::, j]:
-            freq_sum += f
+        for frs in freq[1::, j]:
+            freq_sum += frs
         results[j]._frequencies = freq_sum
 
     return results

From aa7efb1b6814f80e25b9fdc8b890d03e91535d03 Mon Sep 17 00:00:00 2001
From: Alejandro Sopena <alexsopena15@gmail.com>
Date: Mon, 11 Dec 2023 01:52:38 +0100
Subject: [PATCH 21/47] importance clifford sampling

---
 src/qibo/models/error_mitigation.py | 280 ++++++++++++++++------------
 1 file changed, 161 insertions(+), 119 deletions(-)

diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index a202d8ef95..60f6a67433 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -678,31 +678,33 @@ def sample_clifford_training_circuit(
             form (``gates.XYZ``, ``kwargs``). For example, phase gates are used by default:
             ``list((RZ, {'theta':0}), (RZ, {'theta':pi/2}), (RZ, {'theta':pi}), (RZ, {'theta':3*pi/2}))``.
         sigma (float, optional): standard devation of the Gaussian distribution used for sampling.
-        backend (:class:`qibo.backends.abstract.Backend`, optional): Calculation engine.
+        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``.
 
     Returns:
         :class:`qibo.models.Circuit`: The sampled circuit.
     """
-    from qibo.quantum_info import random_clifford
+    from qibo.quantum_info import (  # pylint: disable=import-outside-toplevel
+        random_clifford,
+    )
 
     if backend is None:  # pragma: no cover
         backend = GlobalBackend()
 
-    # Find all the non-Clifford gates
-    gates_to_replace = []
-    for i, gate in enumerate(circuit.queue):
-        if gate.clifford is False and isinstance(gate, gates.M) is False:
-            gates_to_replace.append([i, gate])
-    gates_to_replace = np.array(gates_to_replace, dtype=object)
+    non_clifford_gates = [
+        (i, gate)
+        for i, gate in enumerate(circuit.queue)
+        if not gate.clifford and not isinstance(gate, gates.M)
+    ]
 
-    if len(gates_to_replace) == 0:
+    if not non_clifford_gates:
         raise_error(ValueError, "No non-Clifford gate found, no circuit sampled.")
 
-    # Build the training circuit by substituting the sampled gates
     sampled_circuit = circuit.__class__(**circuit.init_kwargs)
 
     for i, gate in enumerate(circuit.queue):
-        if gate.name == "id_end":  # isinstance(gate, gates.M):
+        if isinstance(gate, gates.M):
             gate_rand = gates.Unitary(
                 random_clifford(1, backend=backend, return_circuit=False),
                 gate.qubits[0],
@@ -713,98 +715,99 @@ def sample_clifford_training_circuit(
             sampled_circuit.add(gate)
 
         else:
-            if i in gates_to_replace[:, 0]:
+            if i in [index for index, _ in non_clifford_gates]:
                 gate = gates.Unitary(
                     random_clifford(1, backend=backend, return_circuit=False),
                     gate.qubits[0],
                 )
                 gate.clifford = True
             sampled_circuit.add(gate)
+
     return sampled_circuit
 
 
-def transpile_circ(circuit, qubit_map, backend):
-    from qibolab.transpilers.unitary_decompositions import u3_decomposition
+def error_sensitive_circuit(circuit, observable, fuse=True, backend=None):
+    """
+    Generates a Clifford circuit that preserves the same circuit frame as the input circuit, and stabilizes the specified Pauli observable.
 
-    if backend.name == "qibolab":
-        new_c = circuit.__class__(backend.platform.nqubits)
-        for gate in circuit.queue:
-            qubits = [qubit_map[j] for j in gate.qubits]
-            if isinstance(gate, gates.M):
-                new_gate = gates.M(*tuple(qubits), **gate.init_kwargs)
-                new_gate.result = gate.result
-                new_c.add(new_gate)
-            elif isinstance(gate, gates.I):
-                new_c.add(gate.__class__(*tuple(qubits), **gate.init_kwargs))
-            else:
-                matrix = gate.matrix()
-                new_c.add(gates.U3(qubits[0], *u3_decomposition(matrix)))
-        return new_c
-    else:
-        return circuit
+    Args:
+        circuit (:class:`qibo.models.Circuit`): input circuit.
+        observable (numpy.ndarray): Pauli observable to be measured.
+        fuse (bool, optional): if True, fuse the single qubits gates in the circuit.
+        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``.
 
+    Returns:
+        :class:`qibo.models.Circuit`: the error sensitive circuit.
+        :class:`qibo.models.Circuit`: the sampled Clifford circuit.
+        list: the list of adjustment gates.
 
-def escircuit(circuit, obs, fuse=True, backend=None):
-    """
-    Error sensitive circuit algorithm.
-    Implement what proposed in https://arxiv.org/abs/2112.06255.
+    Reference:
+        1. Dayue Qin, Yanzhu Chen et al, *Error statistics and scalability of quantum error mitigation formulas*.
+           `arXiv:2112.06255 [quant-ph] <https://arxiv.org/abs/2112.06255>`_.
     """
-    from qibo.quantum_info import comp_basis_to_pauli, random_clifford, vectorization
+    from qibo.quantum_info import (  # pylint: disable=import-outside-toplevel
+        comp_basis_to_pauli,
+        random_clifford,
+        vectorization,
+    )
 
     if backend is None:  # pragma: no cover
         backend = GlobalBackend()
 
-    sign = -1
-    # while sign == -1:
-    circ_cliff = sample_clifford_training_circuit(circuit, backend=backend)
-    c_unitary = circ_cliff.unitary(backend=backend)
-    nqubits = circ_cliff.nqubits
-    U_c2p = comp_basis_to_pauli(nqubits, backend=backend)
-    obs_liouville = vectorization(
-        np.transpose(np.conjugate(c_unitary)) @ obs.matrix @ c_unitary, order="row"
+    sampled_circuit = sample_clifford_training_circuit(circuit, backend=backend)
+    unitary_matrix = sampled_circuit.unitary(backend=backend)
+    num_qubits = sampled_circuit.nqubits
+
+    comp_to_pauli = comp_basis_to_pauli(num_qubits, backend=backend)
+    observable_liouville = vectorization(
+        np.transpose(np.conjugate(unitary_matrix)) @ observable.matrix @ unitary_matrix,
+        order="row",
+        backend=backend,
     )
-    obs_pauli_liouville = U_c2p @ obs_liouville
-    index = np.where(abs(obs_pauli_liouville) >= 1e-5)[0][0]
-    sign = np.sign(obs_pauli_liouville[index])
+    observable_pauli_liouville = comp_to_pauli @ observable_liouville
+
+    index = int(np.where(abs(observable_pauli_liouville) >= 1e-5)[0][0])
 
-    obs1 = list(product(["I", "X", "Y", "Z"], repeat=nqubits))[index]
+    observable_pauli = list(product(["I", "X", "Y", "Z"], repeat=num_qubits))[index]
 
-    paulis = {
-        "I": gates.I(0).matrix(),
-        "X": gates.X(0).matrix(),
-        "Y": gates.Y(0).matrix(),
-        "Z": gates.Z(0).matrix(),
+    pauli_gates = {
+        "I": backend.cast(gates.I(0).matrix(backend=backend)),
+        "X": backend.cast(gates.X(0).matrix(backend=backend)),
+        "Y": backend.cast(gates.Y(0).matrix(backend=backend)),
+        "Z": backend.cast(gates.Z(0).matrix(backend=backend)),
     }
 
-    adjust_gates = []
-    for i in range(nqubits):
-        obs_i = paulis[obs1[i]]
-        R = paulis["I"]
-        while np.any(abs(obs_i - paulis["Z"]) > 1e-5) and np.any(
-            abs(obs_i - paulis["I"]) > 1e-5
+    adjustment_gates = []
+    for i in range(num_qubits):
+        observable_i = pauli_gates[observable_pauli[i]]
+        random_init = pauli_gates["I"]
+        while np.any(abs(observable_i - pauli_gates["Z"]) > 1e-5) and np.any(
+            abs(observable_i - pauli_gates["I"]) > 1e-5
         ):
-            R = random_clifford(1, backend=backend, return_circuit=False)
-            obs_i = np.conjugate(np.transpose(R)) @ paulis[obs1[i]] @ R
+            random_init = random_clifford(1, backend=backend, return_circuit=False)
+            observable_i = (
+                np.conjugate(np.transpose(random_init))
+                @ pauli_gates[observable_pauli[i]]
+                @ random_init
+            )
 
-        adjust_gate = gates.Unitary(R, i)
-        adjust_gate.clifford = True
-        adjust_gates.append(adjust_gate)
+        adjustment_gate = gates.Unitary(random_init, i)
+        adjustment_gate.clifford = True
+        adjustment_gates.append(adjustment_gate)
 
-    circ_cliff1 = circ_cliff.__class__(**circ_cliff.init_kwargs)
+    sensitive_circuit = sampled_circuit.__class__(**sampled_circuit.init_kwargs)
 
-    j = 0
-    for gate in circ_cliff.queue:
-        if gate.name == "id_init":
-            circ_cliff1.add(gate)
-            circ_cliff1.add(adjust_gates[j])
-            j += 1
-        else:
-            circ_cliff1.add(gate)
+    for gate in adjustment_gates:
+        sensitive_circuit.add(gate)
+    for gate in sampled_circuit.queue:
+        sensitive_circuit.add(gate)
 
     if fuse:
-        circ_cliff1 = circ_cliff1.fuse(max_qubits=1)
+        sensitive_circuit = sensitive_circuit.fuse(max_qubits=1)
 
-    return circ_cliff1, circ_cliff, adjust_gates
+    return sensitive_circuit, sampled_circuit, adjustment_gates
 
 
 def ICS(
@@ -819,8 +822,37 @@ def ICS(
     backend=None,
 ):
     """
-    Compute the important Clifford Sampling methos proposed in
-    https://arxiv.org/abs/2112.06255.
+    Computes the Important Clifford Sampling method.
+
+    Args:
+        circuit (:class:`qibo.models.Circuit`): input circuit.
+        observable (numpy.ndarray): the observable to be measured.
+        readout (dict, optional): A dictionary that may contain the following keys:
+            - 'calibration_matrix': numpy.ndarray, used for applying a pre-computed calibration matrix for readout error mitigation.
+            - 'random_ncircuits': int, specifies the number of random circuits to use for the randomized method of readout error mitigation.
+            - 'ibu_iters': int, specifies the number of iterations for the iterative Bayesian update method of readout error mitigation.
+                            If ``None`` the 'inverse' method is used.
+            If provided, the corresponding readout error mitigation method is used.
+        qubit_map (list, optional): the qubit map. If None, a list of range of circuit's qubits is used. Defaults to ``None``.
+        noise_model (qibo.models.noise.Noise, optional): the noise model to be applied. Defaults to ``None``.
+        nshots (int, optional): the number of shots for the circuit execution. Defaults to 10000.
+        n_training_samples (int, optional): the number of training samples. Defaults to 10.
+        full_output (bool, optional): if ``True``, this function returns additional
+            information: ``val``, ``optimal_params``, ``train_val``. Defaults to ``False``.
+        backend (qibo.backends.abstract.Backend, optional): the backend to be used in the execution.
+            If None, it uses the global backend. Defaults to ``None``.
+
+    Returns:
+        mitigated_expectation (float): the mitigated expectated value.
+        mitigated_expectation_std (float): the standard deviation of the mitigated expectated value.
+        dep_param (float): the depolarizing parameter.
+        dep_param_std (float): the standard deviation of the depolarizing parameter.
+        lambda_list (list): the list of the depolarizing parameters.
+        data (dict): the data dictionary containing the noise-free and noisy expectation values obtained with the training circuits.
+
+    Reference:
+        1. Dayue Qin, Yanzhu Chen et al, *Error statistics and scalability of quantum error mitigation formulas*.
+           `arXiv:2112.06255 [quant-ph] <https://arxiv.org/abs/2112.06255>`_.
     """
     if backend is None:  # pragma: no cover
         backend = GlobalBackend()
@@ -828,68 +860,78 @@ def ICS(
     if qubit_map is None:
         qubit_map = list(range(circuit.nqubits))
 
-    training_circs = [
-        escircuit(circuit, observable, backend=backend)[0]
+    training_circuits = [
+        error_sensitive_circuit(circuit, observable, backend=backend)[0]
         for _ in range(n_training_samples)
     ]
 
-    data = {"noise-free": {"-1": [], "1": []}, "noisy": {"-1": [], "1": []}}
-    from qibo.config import log
+    data = {"noise-free": [], "noisy": []}
+    lambda_list = []
+
+    for training_circuit in training_circuits:
+        circuit_result = training_circuit(nshots=nshots)
+        expectation = observable.expectation_from_samples(circuit_result.frequencies())
 
-    a_list = {"-1": [], "1": []}
-    for c in training_circs:
-        circuit_result = c(nshots=nshots)
-        exp = observable.expectation_from_samples(circuit_result.frequencies())
-        # state = c().state()
-        # exp = obs.expectation(state)
         if noise_model is not None and backend.name != "qibolab":
-            c = noise_model.apply(c)
+            training_circuit = noise_model.apply(training_circuit)
 
-        c = transpile_circ(c, qubit_map, backend)
-        circuit_result = backend.execute_circuit(c, nshots=nshots)
+        training_circuit = transpile_circ(training_circuit, qubit_map, backend)
+        circuit_result = backend.execute_circuit(training_circuit, nshots=nshots)
 
         if "calibration_matrix" in readout.keys():
-            circuit_result = apply_readout_mitigation(
-                circuit_result, readout["calibration_matrix"], readout["inv"]
+            circuit_result = apply_cal_mat_readout_mitigation(
+                circuit_result,
+                readout["calibration_matrix"],
+                readout["ibu_iters"],
             )
 
-        exp_noisy = observable.expectation_from_samples(circuit_result.frequencies())
+        noisy_expectation = observable.expectation_from_samples(
+            circuit_result.frequencies()
+        )
 
-        if exp > 0:
-            data["noise-free"]["1"].append(exp)
-            data["noisy"]["1"].append(exp_noisy)
-            a_list["1"].append(1 - exp_noisy / exp)
-        else:
-            data["noise-free"]["-1"].append(exp)
-            data["noisy"]["-1"].append(exp_noisy)
-            a_list["-1"].append(1 - exp_noisy / exp)
-
-    a_std_1 = np.std(a_list["1"])
-    a_std_m1 = np.std(a_list["-1"])
-    a = (np.sum(a_list["1"]) + np.sum(a_list["-1"])) / n_training_samples
-    a_std = np.sqrt(
-        (a_std_1**2 * len(a_list["1"]) + a_std_m1**2 * len(a_list["-1"]))
-        / n_training_samples
-    )
+        data["noise-free"].append(expectation)
+        data["noisy"].append(noisy_expectation)
+        lambda_list.append(1 - noisy_expectation / expectation)
+
+    dep_param = np.mean(lambda_list)
+    dep_param_std = np.std(lambda_list)
 
     circuit = circuit.fuse(max_qubits=1)
+
     if noise_model is not None and backend.name != "qibolab":
         circuit = noise_model.apply(circuit)
+
     circuit = transpile_circ(circuit, qubit_map, backend)
     circuit_result = backend.execute_circuit(circuit, nshots=nshots)
+
     if "calibration_matrix" in readout.keys():
-        circuit_result = apply_readout_mitigation(
-            circuit_result, readout["calibration_matrix"], readout["inv"]
+        circuit_result = apply_cal_mat_readout_mitigation(
+            circuit_result, readout["calibration_matrix"], readout["ibu_iters"]
         )
-    exp_noisy = observable.expectation_from_samples(circuit_result.frequencies())
-    exp_mit = (1 - a) * exp_noisy / ((1 - a) ** 2 + a_std**2)
-    exp_mit_std = (
-        a_std
-        * abs(exp_noisy)
-        * abs((1 - a) ** 2 - a_std**2)
-        / ((1 - a) ** 2 + a_std**2) ** 2
+
+    noisy_expectation = observable.expectation_from_samples(
+        circuit_result.frequencies()
+    )
+    mitigated_expectation = (
+        (1 - dep_param)
+        * noisy_expectation
+        / ((1 - dep_param) ** 2 + dep_param_std**2)
     )
+    mitigated_expectation_std = (
+        dep_param_std
+        * abs(noisy_expectation)
+        * abs((1 - dep_param) ** 2 - dep_param_std**2)
+        / ((1 - dep_param) ** 2 + dep_param_std**2) ** 2
+    )
+
     if full_output:
-        return exp_mit, exp_mit_std, a, a_std, a_list, data
+        return (
+            mitigated_expectation,
+            mitigated_expectation_std,
+            dep_param,
+            dep_param_std,
+            lambda_list,
+            data,
+        )
 
-    return exp_mit
+    return mitigated_expectation

From 626540b92a83e9009501e9061c0175a4336210ae Mon Sep 17 00:00:00 2001
From: Alejandro Sopena <alexsopena15@gmail.com>
Date: Mon, 11 Dec 2023 01:53:47 +0100
Subject: [PATCH 22/47] transpile

---
 src/qibo/models/error_mitigation.py | 23 +++++++++++++++++++++++
 1 file changed, 23 insertions(+)

diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index 60f6a67433..1996bb69f7 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -935,3 +935,26 @@ def ICS(
         )
 
     return mitigated_expectation
+
+
+def transpile_circ(circuit, qubit_map, backend):
+    from qibolab.transpilers.unitary_decompositions import u3_decomposition
+
+    if backend.name == "qibolab":
+        if qubit_map is None:
+            qubit_map = list(range(circuit.nqubits))
+        new_c = circuit.__class__(backend.platform.nqubits)
+        for gate in circuit.queue:
+            qubits = [qubit_map[j] for j in gate.qubits]
+            if isinstance(gate, gates.M):
+                new_gate = gates.M(*tuple(qubits), **gate.init_kwargs)
+                new_gate.result = gate.result
+                new_c.add(new_gate)
+            elif isinstance(gate, gates.I):
+                new_c.add(gate.__class__(*tuple(qubits), **gate.init_kwargs))
+            else:
+                matrix = gate.matrix()
+                new_c.add(gates.U3(qubits[0], *u3_decomposition(matrix)))
+        return new_c
+
+    return circuit

From 6b13d0f5892b522f7ddd03fc625bc972edab048d Mon Sep 17 00:00:00 2001
From: Alejandro Sopena <alexsopena15@gmail.com>
Date: Mon, 11 Dec 2023 01:54:22 +0100
Subject: [PATCH 23/47] updating tests

---
 tests/test_models_error_mitigation.py | 130 +++++++++++++++++++-------
 1 file changed, 95 insertions(+), 35 deletions(-)

diff --git a/tests/test_models_error_mitigation.py b/tests/test_models_error_mitigation.py
index d602d0ae4c..0d19979e60 100644
--- a/tests/test_models_error_mitigation.py
+++ b/tests/test_models_error_mitigation.py
@@ -2,14 +2,17 @@
 import pytest
 
 from qibo import Circuit, gates
+from qibo.backends import construct_backend
 from qibo.hamiltonians import SymbolicHamiltonian
 from qibo.models.error_mitigation import (
     CDR,
+    ICS,
     ZNE,
+    apply_cal_mat_readout_mitigation,
     apply_randomized_readout_mitigation,
-    apply_readout_mitigation,
-    calibration_matrix,
-    sample_training_circuit,
+    get_calibration_matrix,
+    sample_clifford_training_circuit,
+    sample_training_circuit_cdr,
     vnCDR,
 )
 from qibo.noise import DepolarizingError, NoiseModel, ReadoutError
@@ -22,6 +25,7 @@ def get_noise_model(error, gate, cal_matrix=[False, None]):
     noise.add(error, gate)
     if cal_matrix[0]:
         noise.add(ReadoutError(probabilities=cal_matrix[1]), gate=gates.M)
+
     return noise
 
 
@@ -43,10 +47,9 @@ def get_circuit(nqubits):
     c.add(gates.RZ(q + 1, theta=-2 * dt) for q in range(1, nqubits, 2))
     c.add(gates.CNOT(q, q + 1) for q in range(1, nqubits, 2))
     c.add(gates.M(*range(nqubits)))
-    return c
 
+    return c
 
-from qibo.backends import construct_backend
 
 backend = construct_backend("numpy")
 # # Generate random calibration matrices
@@ -66,7 +69,7 @@ def get_circuit(nqubits):
             3,
             get_noise_model(DepolarizingError(0.1), gates.CNOT),
             "CNOT",
-            {"calibration_matrix": cal_matrix_3q},
+            {"calibration_matrix": cal_matrix_3q, "ibu_iters": None},
         ),
         (
             3,
@@ -79,7 +82,13 @@ def get_circuit(nqubits):
             1,
             get_noise_model(DepolarizingError(0.3), gates.RX),
             "RX",
-            {"calibration_matrix": cal_matrix_1q},
+            {"calibration_matrix": cal_matrix_1q, "ibu_iters": None},
+        ),
+        (
+            1,
+            get_noise_model(DepolarizingError(0.3), gates.RX),
+            "RX",
+            {"calibration_matrix": cal_matrix_1q, "ibu_iters": 10},
         ),
         (1, get_noise_model(DepolarizingError(0.1), gates.RX), "RX", {"ncircuits": 2}),
     ],
@@ -116,6 +125,7 @@ def test_zne(backend, nqubits, noise, solve, insertion_gate, readout):
         readout=readout,
         backend=backend,
     )
+
     assert np.abs(exact - estimate) <= np.abs(exact - noisy)
 
 
@@ -127,7 +137,11 @@ def test_zne(backend, nqubits, noise, solve, insertion_gate, readout):
         (get_noise_model(DepolarizingError(0.1), gates.CNOT), {}),
         (
             get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, cal_matrix_3q]),
-            {"calibration_matrix": cal_matrix_3q},
+            {"calibration_matrix": cal_matrix_3q, "ibu_iters": None},
+        ),
+        (
+            get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, cal_matrix_3q]),
+            {"calibration_matrix": cal_matrix_3q, "ibu_iters": 10},
         ),
         (
             get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, cal_matrix_3q]),
@@ -167,6 +181,7 @@ def test_cdr(backend, nqubits, noise, full_output, readout):
     )
     if full_output:
         estimate = estimate[0]
+
     assert np.abs(exact - estimate) <= np.abs(exact - noisy)
 
 
@@ -189,32 +204,22 @@ def test_sample_training_circuit(nqubits):
     c.add(gates.RZ(q + 1, theta=-2 * dt) for q in range(1, nqubits, 2))
     c.add(gates.CNOT(q, q + 1) for q in range(1, nqubits, 2))
     c.add(gates.M(q) for q in range(nqubits))
+
+    for j in range(len(c.queue)):
+        c.queue[j].clifford = True
     with pytest.raises(ValueError):
-        sample_training_circuit(c)
+        sample_training_circuit_cdr(c)
+        sample_clifford_training_circuit(c)
 
 
 @pytest.mark.parametrize(
     "nqubits,noise,insertion_gate,readout",
     [
-        (3, get_noise_model(DepolarizingError(0.1), gates.CNOT), "CNOT", {}),
-        (
-            3,
-            get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, cal_matrix_3q]),
-            "CNOT",
-            {"calibration_matrix": cal_matrix_3q},
-        ),
-        (
-            3,
-            get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, cal_matrix_3q]),
-            "CNOT",
-            {"ncircuits": 2},
-        ),
-        (1, get_noise_model(DepolarizingError(0.1), gates.RX), "RX", {}),
         (
             1,
             get_noise_model(DepolarizingError(0.1), gates.RX, [True, cal_matrix_1q]),
             "RX",
-            {"calibration_matrix": cal_matrix_1q},
+            {"calibration_matrix": cal_matrix_1q, "ibu_iters": 10},
         ),
         (
             1,
@@ -242,12 +247,6 @@ def test_vncdr(backend, nqubits, noise, full_output, insertion_gate, readout):
     # Noise-free expected value
     exact = obs.expectation(backend.execute_circuit(c).state())
     # Noisy expected value without mitigation
-    if "calibration_matrix" in readout.keys() or "ncircuits" in readout.keys():
-        if nqubits == 1:
-            p = cal_matrix_1q
-        elif nqubits == 3:
-            p = cal_matrix_3q
-        # noise.add(ReadoutError(probabilities=p),gate=gates.M)
     state = backend.execute_circuit(noise.apply(c), nshots=10000)
     noisy = state.expectation_from_samples(obs)
     # Mitigated expected value
@@ -265,12 +264,14 @@ def test_vncdr(backend, nqubits, noise, full_output, insertion_gate, readout):
     )
     if full_output:
         estimate = estimate[0]
+
     assert np.abs(exact - estimate) <= np.abs(exact - noisy)
 
 
 @pytest.mark.parametrize("nqubits", [3])
-@pytest.mark.parametrize("method", ["cal_matrix", "temme"])
-def test_readout_mitigation(backend, nqubits, method):
+@pytest.mark.parametrize("method", ["cal_matrix", "randomized"])
+@pytest.mark.parametrize("ibu_iters", [None, 10])
+def test_readout_mitigation(backend, nqubits, method, ibu_iters):
     if backend.name == "tensorflow":
         import tensorflow as tf
 
@@ -283,7 +284,9 @@ def test_readout_mitigation(backend, nqubits, method):
     noise = NoiseModel()
     noise.add(ReadoutError(probabilities=p), gate=gates.M)
     if method == "cal_matrix":
-        calibration = calibration_matrix(nqubits, noise, nshots=nshots, backend=backend)
+        calibration = get_calibration_matrix(
+            nqubits, None, noise, nshots=nshots, backend=backend
+        )
     # Define the observable
     obs = np.prod([Z(i) for i in range(nqubits)])
     obs = SymbolicHamiltonian(obs, backend=backend)
@@ -295,9 +298,9 @@ def test_readout_mitigation(backend, nqubits, method):
     state = backend.execute_circuit(noise.apply(c), nshots=nshots)
     noisy_val = state.expectation_from_samples(obs)
     if method == "cal_matrix":
-        mit_state = apply_readout_mitigation(state, calibration)
+        mit_state = apply_cal_mat_readout_mitigation(state, calibration, ibu_iters)
         mit_val = mit_state.expectation_from_samples(obs)
-    elif method == "temme":
+    elif method == "randomized":
         ncircuits = 10
         result, result_cal = apply_randomized_readout_mitigation(
             c, noise, nshots, ncircuits, backend
@@ -305,4 +308,61 @@ def test_readout_mitigation(backend, nqubits, method):
         mit_val = result.expectation_from_samples(
             obs
         ) / result_cal.expectation_from_samples(obs)
+
     assert np.abs(true_val - mit_val) <= np.abs(true_val - noisy_val)
+
+
+@pytest.mark.parametrize("nqubits", [3])
+@pytest.mark.parametrize("full_output", [False, True])
+@pytest.mark.parametrize(
+    "noise,readout",
+    [
+        (get_noise_model(DepolarizingError(0.1), gates.CNOT), {}),
+        (
+            get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, cal_matrix_3q]),
+            {"calibration_matrix": cal_matrix_3q, "ibu_iters": None},
+        ),
+        (
+            get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, cal_matrix_3q]),
+            {"calibration_matrix": cal_matrix_3q, "ibu_iters": 10},
+        ),
+        (
+            get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, cal_matrix_3q]),
+            {"ncircuits": 2},
+        ),
+    ],
+)
+def test_ics(backend, nqubits, noise, full_output, readout):
+    if backend.name == "tensorflow":
+        import tensorflow as tf
+
+        tf.config.threading.set_inter_op_parallelism_threads = 1
+        tf.config.threading.set_intra_op_parallelism_threads = 1
+    else:
+        backend.set_threads(1)
+    """Test that ICS reduces the noise."""
+    # Define the circuit
+    c = get_circuit(nqubits)
+    # Define the observable
+    obs = np.prod([Z(i) for i in range(nqubits)])
+    obs = SymbolicHamiltonian(obs, backend=backend)
+    # Noise-free expected value
+    exact = obs.expectation(backend.execute_circuit(c).state())
+    # Noisy expected value without mitigation
+    state = backend.execute_circuit(noise.apply(c), nshots=10000)
+    noisy = state.expectation_from_samples(obs)
+    # Mitigated expected value
+    estimate = ICS(
+        circuit=c,
+        observable=obs,
+        noise_model=noise,
+        nshots=10000,
+        n_training_samples=20,
+        full_output=full_output,
+        readout=readout,
+        backend=backend,
+    )
+    if full_output:
+        estimate = estimate[0]
+
+    assert np.abs(exact - estimate) <= np.abs(exact - noisy)

From ea56f406bcab9c234eaf462a679c1195234502d1 Mon Sep 17 00:00:00 2001
From: Alejandro Sopena <alexsopena15@gmail.com>
Date: Mon, 11 Dec 2023 02:22:48 +0100
Subject: [PATCH 24/47] fix docstring

---
 src/qibo/models/error_mitigation.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index c7e3b0c654..0226e1914b 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -553,14 +553,14 @@ def apply_cal_mat_readout_mitigation(state, calibration_matrix, iterations=None)
     Applies readout error mitigation to the given state using the provided calibration matrix.
 
     Args:
-        state (:class:`qibo.states.CircuitResult`): the input state to be updated. This state should contain the
+        state (:class:`qibo.measurements.CircuitResult`): the input state to be updated. This state should contain the
             frequencies that need to be mitigated.
         calibration_matrix (numpy.ndarray, optional): the calibration matrix for readout mitigation.
         iterations (int, optional): the number of iterations to use for the Iterative Bayesian Unfolding method.
             If ``None`` the 'inverse' method is used. Defaults to ``None``.
 
     Returns:
-        :class:`qibo.states.CircuitResult`: The input state with the updated (mitigated) frequencies.
+        :class:`qibo.measurements.CircuitResult`: The input state with the updated (mitigated) frequencies.
     """
     frequencies = np.zeros(2 ** len(state.measurements[0].qubits))
     for key, value in state.frequencies().items():
@@ -938,7 +938,7 @@ def ICS(
 
 
 def transpile_circ(circuit, qubit_map, backend):
-    from qibolab.transpilers.unitary_decompositions import u3_decomposition
+    from qibo.transpiler.unitary_decompositions import u3_decomposition
 
     if backend.name == "qibolab":
         if qubit_map is None:

From 1dd96c78f6acd93ba7e2e163816d092da3f46a07 Mon Sep 17 00:00:00 2001
From: Alejandro Sopena <alexsopena15@gmail.com>
Date: Mon, 11 Dec 2023 02:29:07 +0100
Subject: [PATCH 25/47] fix variational

---
 src/qibo/models/variational.py | 79 ++++++++++++++++++++++++++--------
 1 file changed, 62 insertions(+), 17 deletions(-)

diff --git a/src/qibo/models/variational.py b/src/qibo/models/variational.py
index 46bc417536..d031e58ae0 100644
--- a/src/qibo/models/variational.py
+++ b/src/qibo/models/variational.py
@@ -40,23 +40,45 @@ def __init__(self, circuit, hamiltonian):
 
     def minimize(
         self,
-        initial_parameters,
-        optimizer=None,
-        optimizer_options={},
+        initial_state,
+        method="Powell",
+        jac=None,
+        hess=None,
+        hessp=None,
+        bounds=None,
+        constraints=(),
+        tol=None,
+        callback=None,
+        options=None,
+        compile=False,
+        processes=None,
     ):
         """Search for parameters which minimizes the hamiltonian expectation.
 
         Args:
-            initial_parameters (array or list): a initial guess for the parameters of the
+            initial_state (array): a initial guess for the parameters of the
                 variational circuit.
-            optimizer (qibo.optimizers.Optimizer): Qibo optimizer which is used to
-                minimize the energy [default is `ScipyMinimizer(method="Powell")`].
-            optimizer_options (dict): extra options for the optimizers.
+            method (str): the desired minimization method.
+                See :meth:`qibo.optimizers.optimize` for available optimization
+                methods.
+            jac (dict): Method for computing the gradient vector for scipy optimizers.
+            hess (dict): Method for computing the hessian matrix for scipy optimizers.
+            hessp (callable): Hessian of objective function times an arbitrary
+                vector for scipy optimizers.
+            bounds (sequence or Bounds): Bounds on variables for scipy optimizers.
+            constraints (dict): Constraints definition for scipy optimizers.
+            tol (float): Tolerance of termination for scipy optimizers.
+            callback (callable): Called after each iteration for scipy optimizers.
+            options (dict): a dictionary with options for the different optimizers.
+            compile (bool): whether the TensorFlow graph should be compiled.
+            processes (int): number of processes when using the paralle BFGS method.
 
         Return:
-            output of `qibo.optimizers.Optimizer` classes: [best function value,
-                best set of trainable parameters, complete result object according
-                to specific optimizer configuration].
+            The final expectation value.
+            The corresponding best parameters.
+            The optimization result object. For scipy methods it returns
+            the ``OptimizeResult``, for ``'cma'`` the ``CMAEvolutionStrategy.result``,
+            and for ``'sgd'`` the options used during the optimization.
         """
 
         def _loss(params, circuit, hamiltonian):
@@ -65,15 +87,38 @@ def _loss(params, circuit, hamiltonian):
             final_state = result.state()
             return hamiltonian.expectation(final_state)
 
-        opt = optimizer(
-            initial_parameters=initial_parameters,
-            loss=_loss,
+        if compile:
+            loss = self.hamiltonian.backend.compile(_loss)
+        else:
+            loss = _loss
+
+        if method == "cma":
+            # TODO: check if we can use this shortcut
+            # dtype = getattr(self.hamiltonian.backend.np, self.hamiltonian.backend._dtypes.get('DTYPE'))
+            dtype = self.hamiltonian.backend.np.float64
+            loss = lambda p, c, h: dtype(_loss(p, c, h))
+        elif method != "sgd":
+            loss = lambda p, c, h: self.hamiltonian.backend.to_numpy(_loss(p, c, h))
+
+        result, parameters, extra = self.optimizers.optimize(
+            loss,
+            initial_state,
             args=(self.circuit, self.hamiltonian),
-            **optimizer_options,
+            method=method,
+            jac=jac,
+            hess=hess,
+            hessp=hessp,
+            bounds=bounds,
+            constraints=constraints,
+            tol=tol,
+            callback=callback,
+            options=options,
+            compile=compile,
+            processes=processes,
+            backend=self.hamiltonian.backend,
         )
-
-        best_f, best_p, results = opt.fit()
-        return best_f, best_p, results
+        self.circuit.set_parameters(parameters)
+        return result, parameters, extra
 
     def energy_fluctuation(self, state):
         """

From 60c7299f92c93021478cab0cc6069d8d8c575b3e Mon Sep 17 00:00:00 2001
From: Alejandro Sopena <alexsopena15@gmail.com>
Date: Mon, 11 Dec 2023 02:35:06 +0100
Subject: [PATCH 26/47] update init

---
 src/qibo/models/__init__.py | 9 +--------
 1 file changed, 1 insertion(+), 8 deletions(-)

diff --git a/src/qibo/models/__init__.py b/src/qibo/models/__init__.py
index 7093858901..cf55bf7d52 100644
--- a/src/qibo/models/__init__.py
+++ b/src/qibo/models/__init__.py
@@ -1,14 +1,7 @@
 from qibo.models import hep, tsp
 from qibo.models.circuit import Circuit
 from qibo.models.encodings import unary_encoder
-from qibo.models.error_mitigation import (
-    CDR,
-    ZNE,
-    get_gammas,
-    get_noisy_circuit,
-    sample_training_circuit_cdr,
-    vnCDR,
-)
+from qibo.models.error_mitigation import CDR, ICS, ZNE, vnCDR
 from qibo.models.evolution import AdiabaticEvolution, StateEvolution
 from qibo.models.grover import Grover
 from qibo.models.qft import QFT

From 452661e387844c1b78cc2fa243a03da011b7d9ed Mon Sep 17 00:00:00 2001
From: Alejandro Sopena <alexsopena15@gmail.com>
Date: Mon, 11 Dec 2023 02:44:30 +0100
Subject: [PATCH 27/47] update pre-commit

---
 .pre-commit-config.yaml | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index a8aa2dea52..2d641a95aa 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -9,7 +9,7 @@ repos:
       - id: check-yaml
       - id: debug-statements
   - repo: https://github.com/psf/black
-    rev: 23.9.1
+    rev: 23.11.0
     hooks:
       - id: black
   - repo: https://github.com/pycqa/isort
@@ -22,7 +22,7 @@ repos:
     hooks:
       - id: pyupgrade
   - repo: https://github.com/hadialqattan/pycln
-    rev: v2.2.2
+    rev: v2.4.0
     hooks:
       - id: pycln
         args: [--config=pyproject.toml]

From a4b691a7e2030b7105a39897776b0553a07528a9 Mon Sep 17 00:00:00 2001
From: Alejandro Sopena <alexsopena15@gmail.com>
Date: Thu, 14 Dec 2023 20:34:53 +0100
Subject: [PATCH 28/47] fix documentation tests

---
 doc/source/code-examples/advancedexamples.rst | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/doc/source/code-examples/advancedexamples.rst b/doc/source/code-examples/advancedexamples.rst
index 821b6e1141..4a3fd766f8 100644
--- a/doc/source/code-examples/advancedexamples.rst
+++ b/doc/source/code-examples/advancedexamples.rst
@@ -1293,17 +1293,17 @@ calibration matrix and use it modify the final state after the circuit execution
 
 .. testcode::
 
-   from qibo.models.error_mitigation import apply_readout_mitigation, calibration_matrix
+   from qibo.models.error_mitigation import apply_cal_mat_readout_mitigation, get_calibration_matrix
 
    nshots = 10000
    # compute the calibration matrix
-   calibration = calibration_matrix(
+   calibration = get_calibration_matrix(
        nqubits, backend=backend, noise_model=noise, nshots=nshots
    )
    # execute the circuit
    state = backend.execute_circuit(noise.apply(c), nshots=nshots)
    # mitigate the readout errors
-   mit_state = apply_readout_mitigation(state, calibration)
+   mit_state = apply_cal_mat_readout_mitigation(state, calibration)
    mit_val = mit_state.expectation_from_samples(obs)
    print(mit_val)
    # 0.5945794816381054

From 8fcb3a208b038323f944ea63edfc970eec4a1752 Mon Sep 17 00:00:00 2001
From: Alejandro Sopena <alexsopena15@gmail.com>
Date: Thu, 14 Dec 2023 23:17:42 +0100
Subject: [PATCH 29/47] unify readout mitigation

---
 src/qibo/models/error_mitigation.py | 232 +++++++++++++---------------
 1 file changed, 111 insertions(+), 121 deletions(-)

diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index 0226e1914b..075c1111ab 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -109,6 +109,7 @@ def ZNE(
     solve_for_gammas=False,
     insertion_gate="CNOT",
     readout: dict = {},
+    qubit_map=None,
     backend=None,
 ):
     """Runs the Zero Noise Extrapolation method for error mitigation.
@@ -153,21 +154,15 @@ def ZNE(
         noisy_circuit = get_noisy_circuit(
             circuit, num_insertions, insertion_gate=insertion_gate
         )
-        if "ncircuits" in readout.keys():
-            circuit_result, circuit_result_cal = apply_randomized_readout_mitigation(
-                noisy_circuit, noise_model, nshots, readout["ncircuits"], backend
-            )
-        else:
-            if noise_model is not None and backend.name != "qibolab":
-                noisy_circuit = noise_model.apply(noisy_circuit)
-            circuit_result = backend.execute_circuit(noisy_circuit, nshots=nshots)
-        if "calibration_matrix" in readout.keys() is not None:
-            circuit_result = apply_cal_mat_readout_mitigation(
-                circuit_result, readout["calibration_matrix"], readout["ibu_iters"]
-            )
-        val = circuit_result.expectation_from_samples(observable)
-        if "ncircuits" in readout.keys():
-            val /= circuit_result_cal.expectation_from_samples(observable)
+        val = apply_readout_mitigation(
+            noisy_circuit,
+            observable,
+            noise_model,
+            nshots,
+            readout,
+            qubit_map,
+            backend=backend,
+        )
         expected_values.append(val)
 
     gamma = get_gammas(noise_levels, analytical=solve_for_gammas)
@@ -265,6 +260,7 @@ def CDR(
     n_training_samples: int = 100,
     full_output: bool = False,
     readout: dict = {},
+    qubit_map=None,
     backend=None,
 ):
     """Runs the Clifford Data Regression error mitigation method.
@@ -312,40 +308,16 @@ def CDR(
     for circ in training_circuits:
         val = circ(nshots=nshots).expectation_from_samples(observable)
         train_val["noise-free"].append(val)
-        if "ncircuits" in readout.keys():
-            circuit_result, circuit_result_cal = apply_randomized_readout_mitigation(
-                circ, noise_model, nshots, readout["ncircuits"], backend
-            )
-        else:
-            if noise_model is not None and backend.name != "qibolab":
-                circ = noise_model.apply(circ)
-            circuit_result = backend.execute_circuit(circ, nshots=nshots)
-        if "calibration_matrix" in readout.keys() is not None:
-            circuit_result = apply_cal_mat_readout_mitigation(
-                circuit_result, readout["calibration_matrix"], readout["ibu_iters"]
-            )
-        val = circuit_result.expectation_from_samples(observable)
-        if "ncircuits" in readout.keys():
-            val /= circuit_result_cal.expectation_from_samples(observable)
+        val = apply_readout_mitigation(
+            circ, observable, noise_model, nshots, readout, qubit_map, backend=backend
+        )
         train_val["noisy"].append(val)
 
     optimal_params = curve_fit(model, train_val["noisy"], train_val["noise-free"])[0]
 
-    if "ncircuits" in readout.keys():
-        circuit_result, circuit_result_cal = apply_randomized_readout_mitigation(
-            circuit, noise_model, nshots, readout["ncircuits"], backend
-        )
-    else:
-        if noise_model is not None and backend.name != "qibolab":
-            circuit = noise_model.apply(circuit)
-        circuit_result = backend.execute_circuit(circuit, nshots=nshots)
-    if "calibration_matrix" in readout.keys() is not None:
-        circuit_result = apply_cal_mat_readout_mitigation(
-            circuit_result, readout["calibration_matrix"], readout["ibu_iters"]
-        )
-    val = circuit_result.expectation_from_samples(observable)
-    if "ncircuits" in readout.keys():
-        val /= circuit_result_cal.expectation_from_samples(observable)
+    val = apply_readout_mitigation(
+        circuit, observable, noise_model, nshots, readout, qubit_map, backend=backend
+    )
     mit_val = model(val, *optimal_params)
 
     if full_output is True:
@@ -365,6 +337,7 @@ def vnCDR(
     insertion_gate: str = "CNOT",
     full_output: bool = False,
     readout: dict = {},
+    qubit_map=None,
     backend=None,
 ):
     """Runs the variable-noise Clifford Data Regression error mitigation method.
@@ -419,24 +392,15 @@ def vnCDR(
         train_val["noise-free"].append(val)
         for level in noise_levels:
             noisy_c = get_noisy_circuit(circ, level, insertion_gate=insertion_gate)
-            if "ncircuits" in readout.keys():
-                (
-                    circuit_result,
-                    circuit_result_cal,
-                ) = apply_randomized_readout_mitigation(
-                    noisy_c, noise_model, nshots, readout["ncircuits"], backend
-                )
-            else:
-                if noise_model is not None and backend.name != "qibolab":
-                    noisy_c = noise_model.apply(noisy_c)
-                circuit_result = backend.execute_circuit(noisy_c, nshots=nshots)
-            if "calibration_matrix" in readout.keys():
-                circuit_result = apply_cal_mat_readout_mitigation(
-                    circuit_result, readout["calibration_matrix"], readout["ibu_iters"]
-                )
-            val = circuit_result.expectation_from_samples(observable)
-            if "ncircuits" in readout.keys():
-                val /= circuit_result_cal.expectation_from_samples(observable)
+            val = apply_readout_mitigation(
+                noisy_c,
+                observable,
+                noise_model,
+                nshots,
+                readout,
+                qubit_map,
+                backend=backend,
+            )
             train_val["noisy"].append(val)
 
     noisy_array = np.array(train_val["noisy"]).reshape(-1, len(noise_levels))
@@ -447,21 +411,15 @@ def vnCDR(
     val = []
     for level in noise_levels:
         noisy_c = get_noisy_circuit(circuit, level, insertion_gate=insertion_gate)
-        if "ncircuits" in readout.keys():
-            circuit_result, circuit_result_cal = apply_randomized_readout_mitigation(
-                noisy_c, noise_model, nshots, readout["ncircuits"], backend
-            )
-        else:
-            if noise_model is not None and backend.name != "qibolab":
-                noisy_c = noise_model.apply(noisy_c)
-            circuit_result = backend.execute_circuit(noisy_c, nshots=nshots)
-        if "calibration_matrix" in readout.keys():
-            circuit_result = apply_cal_mat_readout_mitigation(
-                circuit_result, readout["calibration_matrix"], readout["ibu_iters"]
-            )
-        expval = circuit_result.expectation_from_samples(observable)
-        if "ncircuits" in readout.keys():
-            expval /= circuit_result_cal.expectation_from_samples(observable)
+        expval = apply_readout_mitigation(
+            noisy_c,
+            observable,
+            noise_model,
+            nshots,
+            readout,
+            qubit_map,
+            backend=backend,
+        )
         val.append(expval)
 
     mit_val = model(np.array(val).reshape(-1, 1), *optimal_params[0])[0]
@@ -500,7 +458,7 @@ def iterative_bayesian_unfolding(probabilities, response_matrix, iterations=10):
 
 
 def get_calibration_matrix(
-    nqubits, qubit_map, noise_model=None, nshots: int = 10000, backend=None
+    nqubits, qubit_map=None, noise_model=None, nshots: int = 10000, backend=None
 ):
     """Computes the calibration matrix for readout mitigation.
 
@@ -534,11 +492,11 @@ def get_calibration_matrix(
                 circuit.add(gates.X(qubit))
         circuit.add(gates.M(*range(nqubits)))
 
-        if noise_model is not None and backend.name != "qibolab":
-            circuit = noise_model.apply(circuit)
-        circuit = transpile_circ(circuit, qubit_map, backend)
+        circuit_result = _circuit_conf(
+            circuit, qubit_map, noise_model, nshots, backend=backend
+        )
 
-        frequencies = backend.execute_circuit(circuit, nshots=nshots).frequencies()
+        frequencies = circuit_result.frequencies()
 
         column = np.zeros(2**nqubits)
         for key, value in frequencies.items():
@@ -590,7 +548,12 @@ def apply_cal_mat_readout_mitigation(state, calibration_matrix, iterations=None)
 
 
 def apply_randomized_readout_mitigation(
-    circuit, noise_model=None, nshots: int = int(1e3), ncircuits: int = 10, backend=None
+    circuit,
+    noise_model=None,
+    nshots: int = int(1e3),
+    ncircuits: int = 10,
+    qubit_map=None,
+    backend=None,
 ):
     """Readout mitigation method that transforms the bias in an expectation value into a measurable multiplicative factor. This factor can be eliminated at the expense of increased sampling complexity for the observable.
 
@@ -646,9 +609,10 @@ def apply_randomized_readout_mitigation(
         for circ in circuits:
             circ.add(x_gate)
             circ.add(gates.M(*meas_qubits))
-            if noise_model is not None and backend.name != "qibolab":
-                circ = noise_model.apply(circ)
-            result = backend.execute_circuit(circ, nshots=nshots_r)
+
+            result = _circuit_conf(
+                circ, qubit_map, noise_model, nshots_r, backend=backend
+            )
             result._samples = result.apply_bitflips(error_map)
             results.append(result)
             freqs.append(result.frequencies(binary=False))
@@ -664,6 +628,39 @@ def apply_randomized_readout_mitigation(
     return results
 
 
+def apply_readout_mitigation(
+    circuit,
+    observable,
+    noise_model=None,
+    nshots: int = int(1e3),
+    readout={},
+    qubit_map=None,
+    backend=None,
+):
+    if "ncircuits" in readout:
+        circuit_result, circuit_result_cal = apply_randomized_readout_mitigation(
+            circuit, noise_model, nshots, readout["ncircuits"], backend
+        )
+    else:
+        circuit_result = _circuit_conf(
+            circuit, qubit_map, noise_model, nshots, backend=backend
+        )
+
+    if "calibration_matrix" in readout:
+        circuit_result = apply_cal_mat_readout_mitigation(
+            circuit_result,
+            readout["calibration_matrix"],
+            readout.get("ibu_iters", None),
+        )
+
+    exp_val = circuit_result.expectation_from_samples(observable)
+
+    if "ncircuits" in readout:
+        exp_val /= circuit_result_cal.expectation_from_samples(observable)
+
+    return exp_val
+
+
 def sample_clifford_training_circuit(
     circuit,
     backend=None,
@@ -726,7 +723,7 @@ def sample_clifford_training_circuit(
     return sampled_circuit
 
 
-def error_sensitive_circuit(circuit, observable, fuse=True, backend=None):
+def error_sensitive_circuit(circuit, observable, backend=None):
     """
     Generates a Clifford circuit that preserves the same circuit frame as the input circuit, and stabilizes the specified Pauli observable.
 
@@ -804,9 +801,6 @@ def error_sensitive_circuit(circuit, observable, fuse=True, backend=None):
     for gate in sampled_circuit.queue:
         sensitive_circuit.add(gate)
 
-    if fuse:
-        sensitive_circuit = sensitive_circuit.fuse(max_qubits=1)
-
     return sensitive_circuit, sampled_circuit, adjustment_gates
 
 
@@ -872,21 +866,14 @@ def ICS(
         circuit_result = training_circuit(nshots=nshots)
         expectation = observable.expectation_from_samples(circuit_result.frequencies())
 
-        if noise_model is not None and backend.name != "qibolab":
-            training_circuit = noise_model.apply(training_circuit)
-
-        training_circuit = transpile_circ(training_circuit, qubit_map, backend)
-        circuit_result = backend.execute_circuit(training_circuit, nshots=nshots)
-
-        if "calibration_matrix" in readout.keys():
-            circuit_result = apply_cal_mat_readout_mitigation(
-                circuit_result,
-                readout["calibration_matrix"],
-                readout["ibu_iters"],
-            )
-
-        noisy_expectation = observable.expectation_from_samples(
-            circuit_result.frequencies()
+        noisy_expectation = apply_readout_mitigation(
+            training_circuit,
+            observable,
+            noise_model,
+            nshots,
+            readout,
+            qubit_map,
+            backend=backend,
         )
 
         data["noise-free"].append(expectation)
@@ -896,22 +883,10 @@ def ICS(
     dep_param = np.mean(lambda_list)
     dep_param_std = np.std(lambda_list)
 
-    circuit = circuit.fuse(max_qubits=1)
-
-    if noise_model is not None and backend.name != "qibolab":
-        circuit = noise_model.apply(circuit)
-
-    circuit = transpile_circ(circuit, qubit_map, backend)
-    circuit_result = backend.execute_circuit(circuit, nshots=nshots)
-
-    if "calibration_matrix" in readout.keys():
-        circuit_result = apply_cal_mat_readout_mitigation(
-            circuit_result, readout["calibration_matrix"], readout["ibu_iters"]
-        )
-
-    noisy_expectation = observable.expectation_from_samples(
-        circuit_result.frequencies()
+    noisy_expectation = apply_readout_mitigation(
+        circuit, observable, noise_model, nshots, readout, qubit_map, backend=backend
     )
+
     mitigated_expectation = (
         (1 - dep_param)
         * noisy_expectation
@@ -958,3 +933,18 @@ def transpile_circ(circuit, qubit_map, backend):
         return new_c
 
     return circuit
+
+
+def _circuit_conf(circuit, qubit_map, noise_model, nshots, backend=None):
+    from qibo.transpiler.placer import Custom
+
+    if backend.name == "qibolab":
+        backend.transpiler.passes[1] = Custom(
+            map=qubit_map, connectivity=backend.platform.topology
+        )
+    elif noise_model is not None:
+        circuit = noise_model.apply(circuit)
+
+    circuit_result = backend.execute_circuit(circuit, nshots=nshots)
+
+    return circuit_result

From e0e8c716bc0a8efe3620f36de51b27c8d6eb715e Mon Sep 17 00:00:00 2001
From: Alejandro Sopena <alexsopena15@gmail.com>
Date: Fri, 15 Dec 2023 00:22:31 +0100
Subject: [PATCH 30/47] fix coverage

---
 src/qibo/models/error_mitigation.py   | 56 +++++++++++++++------------
 tests/test_models_error_mitigation.py | 19 +++------
 2 files changed, 37 insertions(+), 38 deletions(-)

diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index 075c1111ab..1f66cf2a12 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -597,11 +597,12 @@ def apply_randomized_readout_mitigation(
 
         error_map = {}
         for j, gate in enumerate(x_gate):
-            if gate.name == "x":
-                if gate.qubits[0] in meas_qubits:
-                    error_map[gate.qubits[0]] = 1
-                else:
-                    x_gate.queue[j] = gates.I(gate.qubits[0])
+            if isinstance(gate, gates.X) and gate.qubits[0] in meas_qubits:
+                error_map[gate.qubits[0]] = 1
+                # if gate.qubits[0] in meas_qubits:
+                #     error_map[gate.qubits[0]] = 1
+                # else:
+                #     x_gate.queue[j] = gates.I(gate.qubits[0])
 
         circuits = [circuit_c, cal_circuit]
         results = []
@@ -637,6 +638,9 @@ def apply_readout_mitigation(
     qubit_map=None,
     backend=None,
 ):
+    if backend is None:  # pragma: no cover
+        backend = GlobalBackend()
+
     if "ncircuits" in readout:
         circuit_result, circuit_result_cal = apply_randomized_readout_mitigation(
             circuit, noise_model, nshots, readout["ncircuits"], backend
@@ -912,33 +916,35 @@ def ICS(
     return mitigated_expectation
 
 
-def transpile_circ(circuit, qubit_map, backend):
-    from qibo.transpiler.unitary_decompositions import u3_decomposition
+# def transpile_circ(circuit, qubit_map, backend):
+#     from qibo.transpiler.unitary_decompositions import u3_decomposition
 
-    if backend.name == "qibolab":
-        if qubit_map is None:
-            qubit_map = list(range(circuit.nqubits))
-        new_c = circuit.__class__(backend.platform.nqubits)
-        for gate in circuit.queue:
-            qubits = [qubit_map[j] for j in gate.qubits]
-            if isinstance(gate, gates.M):
-                new_gate = gates.M(*tuple(qubits), **gate.init_kwargs)
-                new_gate.result = gate.result
-                new_c.add(new_gate)
-            elif isinstance(gate, gates.I):
-                new_c.add(gate.__class__(*tuple(qubits), **gate.init_kwargs))
-            else:
-                matrix = gate.matrix()
-                new_c.add(gates.U3(qubits[0], *u3_decomposition(matrix)))
-        return new_c
+#     if backend.name == "qibolab":
+#         if qubit_map is None:
+#             qubit_map = list(range(circuit.nqubits))
+#         new_c = circuit.__class__(backend.platform.nqubits)
+#         for gate in circuit.queue:
+#             qubits = [qubit_map[j] for j in gate.qubits]
+#             if isinstance(gate, gates.M):
+#                 new_gate = gates.M(*tuple(qubits), **gate.init_kwargs)
+#                 new_gate.result = gate.result
+#                 new_c.add(new_gate)
+#             elif isinstance(gate, gates.I):
+#                 new_c.add(gate.__class__(*tuple(qubits), **gate.init_kwargs))
+#             else:
+#                 matrix = gate.matrix()
+#                 new_c.add(gates.U3(qubits[0], *u3_decomposition(matrix)))
+#         return new_c
 
-    return circuit
+#     return circuit
 
 
 def _circuit_conf(circuit, qubit_map, noise_model, nshots, backend=None):
     from qibo.transpiler.placer import Custom
 
-    if backend.name == "qibolab":
+    if backend is None:  # pragma: no cover
+        backend = GlobalBackend()
+    elif backend.name == "qibolab":  # pragma: no cover
         backend.transpiler.passes[1] = Custom(
             map=qubit_map, connectivity=backend.platform.topology
         )
diff --git a/tests/test_models_error_mitigation.py b/tests/test_models_error_mitigation.py
index 0d19979e60..78a9a1d830 100644
--- a/tests/test_models_error_mitigation.py
+++ b/tests/test_models_error_mitigation.py
@@ -8,8 +8,7 @@
     CDR,
     ICS,
     ZNE,
-    apply_cal_mat_readout_mitigation,
-    apply_randomized_readout_mitigation,
+    apply_readout_mitigation,
     get_calibration_matrix,
     sample_clifford_training_circuit,
     sample_training_circuit_cdr,
@@ -287,6 +286,9 @@ def test_readout_mitigation(backend, nqubits, method, ibu_iters):
         calibration = get_calibration_matrix(
             nqubits, None, noise, nshots=nshots, backend=backend
         )
+        readout = {"calibration_matrix": calibration, "ibu_iters": ibu_iters}
+    elif method == "randomized":
+        readout = {"ncircuits": 10}
     # Define the observable
     obs = np.prod([Z(i) for i in range(nqubits)])
     obs = SymbolicHamiltonian(obs, backend=backend)
@@ -297,17 +299,8 @@ def test_readout_mitigation(backend, nqubits, method, ibu_iters):
     # get noisy expected val
     state = backend.execute_circuit(noise.apply(c), nshots=nshots)
     noisy_val = state.expectation_from_samples(obs)
-    if method == "cal_matrix":
-        mit_state = apply_cal_mat_readout_mitigation(state, calibration, ibu_iters)
-        mit_val = mit_state.expectation_from_samples(obs)
-    elif method == "randomized":
-        ncircuits = 10
-        result, result_cal = apply_randomized_readout_mitigation(
-            c, noise, nshots, ncircuits, backend
-        )
-        mit_val = result.expectation_from_samples(
-            obs
-        ) / result_cal.expectation_from_samples(obs)
+
+    mit_val = apply_readout_mitigation(c, obs, noise, nshots, readout, backend=backend)
 
     assert np.abs(true_val - mit_val) <= np.abs(true_val - noisy_val)
 

From 1443cb67b45542d49a2a0fdf6d03562795db8d15 Mon Sep 17 00:00:00 2001
From: Alejandro Sopena <alexsopena15@gmail.com>
Date: Fri, 15 Dec 2023 02:38:55 +0100
Subject: [PATCH 31/47] fix coverage (raise error)

---
 tests/test_models_error_mitigation.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/tests/test_models_error_mitigation.py b/tests/test_models_error_mitigation.py
index 78a9a1d830..0dc5fd395c 100644
--- a/tests/test_models_error_mitigation.py
+++ b/tests/test_models_error_mitigation.py
@@ -208,6 +208,7 @@ def test_sample_training_circuit(nqubits):
         c.queue[j].clifford = True
     with pytest.raises(ValueError):
         sample_training_circuit_cdr(c)
+    with pytest.raises(ValueError):
         sample_clifford_training_circuit(c)
 
 

From 2b27317bb8ebe8305d56c807609a3608dfdfe443 Mon Sep 17 00:00:00 2001
From: Alejandro Sopena <alexsopena15@gmail.com>
Date: Wed, 20 Dec 2023 07:44:47 +0100
Subject: [PATCH 32/47] updating tests

---
 src/qibo/models/error_mitigation.py   | 39 +++-----------
 tests/test_models_error_mitigation.py | 77 ++++++++++++++++-----------
 2 files changed, 54 insertions(+), 62 deletions(-)

diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index 1f66cf2a12..4731d43a2c 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -706,15 +706,14 @@ def sample_clifford_training_circuit(
 
     for i, gate in enumerate(circuit.queue):
         if isinstance(gate, gates.M):
-            gate_rand = gates.Unitary(
-                random_clifford(1, backend=backend, return_circuit=False),
-                gate.qubits[0],
-            )
-            gate_rand.clifford = True
-
-            sampled_circuit.add(gate_rand)
+            for q in gate.qubits:
+                gate_rand = gates.Unitary(
+                    random_clifford(1, backend=backend, return_circuit=False),
+                    q,
+                )
+                gate_rand.clifford = True
+                sampled_circuit.add(gate_rand)
             sampled_circuit.add(gate)
-
         else:
             if i in [index for index, _ in non_clifford_gates]:
                 gate = gates.Unitary(
@@ -762,6 +761,7 @@ def error_sensitive_circuit(circuit, observable, backend=None):
     num_qubits = sampled_circuit.nqubits
 
     comp_to_pauli = comp_basis_to_pauli(num_qubits, backend=backend)
+    observable.nqubits = num_qubits
     observable_liouville = vectorization(
         np.transpose(np.conjugate(unitary_matrix)) @ observable.matrix @ unitary_matrix,
         order="row",
@@ -916,29 +916,6 @@ def ICS(
     return mitigated_expectation
 
 
-# def transpile_circ(circuit, qubit_map, backend):
-#     from qibo.transpiler.unitary_decompositions import u3_decomposition
-
-#     if backend.name == "qibolab":
-#         if qubit_map is None:
-#             qubit_map = list(range(circuit.nqubits))
-#         new_c = circuit.__class__(backend.platform.nqubits)
-#         for gate in circuit.queue:
-#             qubits = [qubit_map[j] for j in gate.qubits]
-#             if isinstance(gate, gates.M):
-#                 new_gate = gates.M(*tuple(qubits), **gate.init_kwargs)
-#                 new_gate.result = gate.result
-#                 new_c.add(new_gate)
-#             elif isinstance(gate, gates.I):
-#                 new_c.add(gate.__class__(*tuple(qubits), **gate.init_kwargs))
-#             else:
-#                 matrix = gate.matrix()
-#                 new_c.add(gates.U3(qubits[0], *u3_decomposition(matrix)))
-#         return new_c
-
-#     return circuit
-
-
 def _circuit_conf(circuit, qubit_map, noise_model, nshots, backend=None):
     from qibo.transpiler.placer import Custom
 
diff --git a/tests/test_models_error_mitigation.py b/tests/test_models_error_mitigation.py
index 0dc5fd395c..6951b4e6e0 100644
--- a/tests/test_models_error_mitigation.py
+++ b/tests/test_models_error_mitigation.py
@@ -28,7 +28,9 @@ def get_noise_model(error, gate, cal_matrix=[False, None]):
     return noise
 
 
-def get_circuit(nqubits):
+def get_circuit(nqubits, nmeas=None):
+    if nmeas is None:
+        nmeas = nqubits
     # Define the circuit
     hz = 0.5
     hx = 0.5
@@ -45,7 +47,7 @@ def get_circuit(nqubits):
     c.add(gates.CNOT(q, q + 1) for q in range(1, nqubits, 2))
     c.add(gates.RZ(q + 1, theta=-2 * dt) for q in range(1, nqubits, 2))
     c.add(gates.CNOT(q, q + 1) for q in range(1, nqubits, 2))
-    c.add(gates.M(*range(nqubits)))
+    c.add(gates.M(*range(nmeas)))
 
     return c
 
@@ -55,8 +57,8 @@ def get_circuit(nqubits):
 cal_matrix_1q = random_stochastic_matrix(
     2, diagonally_dominant=True, seed=2, backend=backend
 )
-cal_matrix_3q = random_stochastic_matrix(
-    8, diagonally_dominant=True, seed=2, backend=backend
+cal_matrix_2q = random_stochastic_matrix(
+    4, diagonally_dominant=True, seed=2, backend=backend
 )
 
 
@@ -68,7 +70,7 @@ def get_circuit(nqubits):
             3,
             get_noise_model(DepolarizingError(0.1), gates.CNOT),
             "CNOT",
-            {"calibration_matrix": cal_matrix_3q, "ibu_iters": None},
+            {"calibration_matrix": cal_matrix_2q, "ibu_iters": None},
         ),
         (
             3,
@@ -102,13 +104,19 @@ def test_zne(backend, nqubits, noise, solve, insertion_gate, readout):
         tf.config.threading.set_intra_op_parallelism_threads = 1
     else:
         backend.set_threads(1)
+
+    if nqubits == 1:
+        nmeas = 1
+    else:
+        nmeas = nqubits - 1
     # Define the circuit
-    c = get_circuit(nqubits)
+    c = get_circuit(nqubits, nmeas)
     # Define the observable
-    obs = np.prod([Z(i) for i in range(nqubits)])
+    obs = np.prod([Z(i) for i in range(nmeas)])
+    obs_exact = SymbolicHamiltonian(obs, nqubits=nqubits, backend=backend)
     obs = SymbolicHamiltonian(obs, backend=backend)
     # Noise-free expected value
-    exact = obs.expectation(backend.execute_circuit(c).state())
+    exact = obs_exact.expectation(backend.execute_circuit(c).state())
     # Noisy expected value without mitigation
     state = backend.execute_circuit(noise.apply(c), nshots=10000)
     noisy = state.expectation_from_samples(obs)
@@ -135,20 +143,21 @@ def test_zne(backend, nqubits, noise, solve, insertion_gate, readout):
     [
         (get_noise_model(DepolarizingError(0.1), gates.CNOT), {}),
         (
-            get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, cal_matrix_3q]),
-            {"calibration_matrix": cal_matrix_3q, "ibu_iters": None},
+            get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, cal_matrix_2q]),
+            {"calibration_matrix": cal_matrix_2q, "ibu_iters": None},
         ),
         (
-            get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, cal_matrix_3q]),
-            {"calibration_matrix": cal_matrix_3q, "ibu_iters": 10},
+            get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, cal_matrix_2q]),
+            {"calibration_matrix": cal_matrix_2q, "ibu_iters": 10},
         ),
         (
-            get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, cal_matrix_3q]),
+            get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, cal_matrix_2q]),
             {"ncircuits": 2},
         ),
     ],
 )
 def test_cdr(backend, nqubits, noise, full_output, readout):
+    """Test that CDR reduces the noise."""
     if backend.name == "tensorflow":
         import tensorflow as tf
 
@@ -156,14 +165,19 @@ def test_cdr(backend, nqubits, noise, full_output, readout):
         tf.config.threading.set_intra_op_parallelism_threads = 1
     else:
         backend.set_threads(1)
-    """Test that CDR reduces the noise."""
+
+    if nqubits == 1:
+        nmeas = 1
+    else:
+        nmeas = nqubits - 1
     # Define the circuit
-    c = get_circuit(nqubits)
+    c = get_circuit(nqubits, nmeas)
     # Define the observable
-    obs = np.prod([Z(i) for i in range(nqubits)])
+    obs = np.prod([Z(i) for i in range(nmeas)])
+    obs_exact = SymbolicHamiltonian(obs, nqubits=nqubits, backend=backend)
     obs = SymbolicHamiltonian(obs, backend=backend)
     # Noise-free expected value
-    exact = obs.expectation(backend.execute_circuit(c).state())
+    exact = obs_exact.expectation(backend.execute_circuit(c).state())
     # Noisy expected value without mitigation
     state = backend.execute_circuit(noise.apply(c), nshots=10000)
     noisy = state.expectation_from_samples(obs)
@@ -268,10 +282,10 @@ def test_vncdr(backend, nqubits, noise, full_output, insertion_gate, readout):
     assert np.abs(exact - estimate) <= np.abs(exact - noisy)
 
 
-@pytest.mark.parametrize("nqubits", [3])
+@pytest.mark.parametrize("nqubits,nmeas", [(3, 2)])
 @pytest.mark.parametrize("method", ["cal_matrix", "randomized"])
 @pytest.mark.parametrize("ibu_iters", [None, 10])
-def test_readout_mitigation(backend, nqubits, method, ibu_iters):
+def test_readout_mitigation(backend, nqubits, nmeas, method, ibu_iters):
     if backend.name == "tensorflow":
         import tensorflow as tf
 
@@ -280,21 +294,21 @@ def test_readout_mitigation(backend, nqubits, method, ibu_iters):
     else:
         backend.set_threads(1)
     nshots = 10000
-    p = random_stochastic_matrix(2**nqubits, diagonally_dominant=True, seed=5)
+    p = random_stochastic_matrix(2**nmeas, diagonally_dominant=True, seed=5)
     noise = NoiseModel()
     noise.add(ReadoutError(probabilities=p), gate=gates.M)
     if method == "cal_matrix":
         calibration = get_calibration_matrix(
-            nqubits, None, noise, nshots=nshots, backend=backend
+            nmeas, None, noise, nshots=nshots, backend=backend
         )
         readout = {"calibration_matrix": calibration, "ibu_iters": ibu_iters}
     elif method == "randomized":
         readout = {"ncircuits": 10}
     # Define the observable
-    obs = np.prod([Z(i) for i in range(nqubits)])
+    obs = np.prod([Z(i) for i in range(nmeas)])
     obs = SymbolicHamiltonian(obs, backend=backend)
     # get noise free expected val
-    c = get_circuit(nqubits)
+    c = get_circuit(nqubits, nmeas)
     true_state = backend.execute_circuit(c, nshots=nshots)
     true_val = true_state.expectation_from_samples(obs)
     # get noisy expected val
@@ -313,15 +327,15 @@ def test_readout_mitigation(backend, nqubits, method, ibu_iters):
     [
         (get_noise_model(DepolarizingError(0.1), gates.CNOT), {}),
         (
-            get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, cal_matrix_3q]),
-            {"calibration_matrix": cal_matrix_3q, "ibu_iters": None},
+            get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, cal_matrix_2q]),
+            {"calibration_matrix": cal_matrix_2q, "ibu_iters": None},
         ),
         (
-            get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, cal_matrix_3q]),
-            {"calibration_matrix": cal_matrix_3q, "ibu_iters": 10},
+            get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, cal_matrix_2q]),
+            {"calibration_matrix": cal_matrix_2q, "ibu_iters": 10},
         ),
         (
-            get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, cal_matrix_3q]),
+            get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, cal_matrix_2q]),
             {"ncircuits": 2},
         ),
     ],
@@ -336,12 +350,13 @@ def test_ics(backend, nqubits, noise, full_output, readout):
         backend.set_threads(1)
     """Test that ICS reduces the noise."""
     # Define the circuit
-    c = get_circuit(nqubits)
+    c = get_circuit(nqubits, nqubits - 1)
     # Define the observable
-    obs = np.prod([Z(i) for i in range(nqubits)])
+    obs = np.prod([Z(i) for i in range(nqubits - 1)])
+    obs_exact = SymbolicHamiltonian(obs, nqubits=nqubits, backend=backend)
     obs = SymbolicHamiltonian(obs, backend=backend)
     # Noise-free expected value
-    exact = obs.expectation(backend.execute_circuit(c).state())
+    exact = obs_exact.expectation(backend.execute_circuit(c).state())
     # Noisy expected value without mitigation
     state = backend.execute_circuit(noise.apply(c), nshots=10000)
     noisy = state.expectation_from_samples(obs)

From abdd6234e5b250ed85e88478bf525aa43ae31935 Mon Sep 17 00:00:00 2001
From: Alejandro Sopena <alexsopena15@gmail.com>
Date: Wed, 20 Dec 2023 08:35:54 +0100
Subject: [PATCH 33/47] updating doc strings

---
 src/qibo/models/error_mitigation.py   | 108 +++++++++++++++++---------
 tests/test_models_error_mitigation.py |   4 +-
 2 files changed, 72 insertions(+), 40 deletions(-)

diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index 4731d43a2c..93792fac49 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -105,7 +105,7 @@ def ZNE(
     observable,
     noise_levels,
     noise_model=None,
-    nshots=int(1e4),
+    nshots=10000,
     solve_for_gammas=False,
     insertion_gate="CNOT",
     readout: dict = {},
@@ -120,11 +120,11 @@ def ZNE(
 
     Args:
         circuit (:class:`qibo.models.Circuit`): input circuit.
-        observable (numpy.ndarray): Observable to measure.
+        observable (:class:`qibo.hamiltonians.Hamiltonian/:class:`qibo.hamiltonians.SymbolicHamiltonian`): Observable to measure.
         noise_levels (numpy.ndarray): Sequence of noise levels.
         noise_model (:class:`qibo.noise.NoiseModel`, optional): Noise model applied
             to simulate noisy computation.
-        nshots (int, optional): Number of shots. Defauylts to 10000.
+        nshots (int, optional): Number of shots. Defaults to :math:`10000`.
         solve_for_gammas (bool, optional): If ``True``, explicitly solve the
             equations to obtain the ``gamma`` coefficients. Default is ``False``.
         insertion_gate (str, optional): gate to be used in the insertion.
@@ -133,8 +133,9 @@ def ZNE(
         readout (dict, optional): A dictionary that may contain the following keys:
             - 'calibration_matrix': numpy.ndarray, used for applying a pre-computed calibration matrix for readout error mitigation.
             - 'random_ncircuits': int, specifies the number of random circuits to use for the randomized method of readout error mitigation.
-            - 'ibu_iters': int, specifies the number of iterations for the iterative Bayesian update method of readout error mitigation.
+            - 'ibu_iters': int, specifies the number of iterations for the iterative Bayesian unfolding method of readout error mitigation.
             If provided, the corresponding readout error mitigation method is used.
+        qubit_map (list, optional): the qubit map. If None, a list of range of circuit's qubits is used. Defaults to ``None``.
         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``.
@@ -255,7 +256,7 @@ def CDR(
     circuit,
     observable,
     noise_model,
-    nshots: int = int(1e4),
+    nshots: int = 10000,
     model=lambda x, a, b: a * x + b,
     n_training_samples: int = 100,
     full_output: bool = False,
@@ -268,10 +269,10 @@ def CDR(
     Args:
         circuit (:class:`qibo.models.Circuit`): input circuit decomposed in the
             primitive gates ``X``, ``CNOT``, ``RX(pi/2)``, ``RZ(theta)``.
-        observable (numpy.ndarray): observable to be measured.
+        observable (:class:`qibo.hamiltonians.Hamiltonian/:class:`qibo.hamiltonians.SymbolicHamiltonian`): observable to be measured.
         noise_model (:class:`qibo.noise.NoiseModel`): noise model used for simulating
             noisy computation.
-        nshots (int, optional): number of shots. Defaults 10000.
+        nshots (int, optional): number of shots. Defaults :math:`10000`.
         model (callable, optional): model used for fitting. This should be a callable
             function object ``f(x, *params)``, taking as input the predictor variable
             and the parameters. Default is a simple linear model ``f(x,a,b) := a*x + b``.
@@ -281,8 +282,9 @@ def CDR(
         readout (dict, optional): A dictionary that may contain the following keys:
             - 'calibration_matrix': numpy.ndarray, used for applying a pre-computed calibration matrix for readout error mitigation.
             - 'random_ncircuits': int, specifies the number of random circuits to use for the randomized method of readout error mitigation.
-            - 'ibu_iters': int, specifies the number of iterations for the iterative Bayesian update method of readout error mitigation.
+            - 'ibu_iters': int, specifies the number of iterations for the iterative Bayesian unfolding method of readout error mitigation.
             If provided, the corresponding readout error mitigation method is used.
+        qubit_map (list, optional): the qubit map. If None, a list of range of circuit's qubits is used. Defaults to ``None``.
         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``.
@@ -331,7 +333,7 @@ def vnCDR(
     observable,
     noise_levels,
     noise_model,
-    nshots: int = int(1e4),
+    nshots: int = 10000,
     model=lambda x, *params: (x * np.array(params).reshape(-1, 1)).sum(0),
     n_training_samples: int = 100,
     insertion_gate: str = "CNOT",
@@ -345,11 +347,11 @@ def vnCDR(
     Args:
         circuit (:class:`qibo.models.Circuit`): input circuit decomposed in the
             primitive gates ``X``, ``CNOT``, ``RX(pi/2)``, ``RZ(theta)``.
-        observable (numpy.ndarray): observable to be measured.
+        observable (:class:`qibo.hamiltonians.Hamiltonian/:class:`qibo.hamiltonians.SymbolicHamiltonian`): observable to be measured.
         noise_levels (numpy.ndarray): sequence of noise levels.
         noise_model (:class:`qibo.noise.NoiseModel`): noise model used for
             simulating noisy computation.
-        nshots (int, optional): number of shots.
+        nshots (int, optional): number of shots. Defaults to :math:`10000`.
         model (callable, optional): model used for fitting. This should be a callable
             function object ``f(x, *params)``, taking as input the predictor variable
             and the parameters. Default is a simple linear model ``f(x,a,b) := a*x + b``.
@@ -362,8 +364,9 @@ def vnCDR(
         readout (dict, optional): A dictionary that may contain the following keys:
             - 'calibration_matrix': numpy.ndarray, used for applying a pre-computed calibration matrix for readout error mitigation.
             - 'random_ncircuits': int, specifies the number of random circuits to use for the randomized method of readout error mitigation.
-            - 'ibu_iters': int, specifies the number of iterations for the iterative Bayesian update method of readout error mitigation.
+            - 'ibu_iters': int, specifies the number of iterations for the iterative Bayesian unfolding method of readout error mitigation.
             If provided, the corresponding readout error mitigation method is used.
+        qubit_map (list, optional): the qubit map. If None, a list of range of circuit's qubits is used. Defaults to ``None``.
         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``.
@@ -457,17 +460,18 @@ def iterative_bayesian_unfolding(probabilities, response_matrix, iterations=10):
     return unfolded_probabilities
 
 
-def get_calibration_matrix(
+def get_response_matrix(
     nqubits, qubit_map=None, noise_model=None, nshots: int = 10000, backend=None
 ):
     """Computes the calibration matrix for readout mitigation.
 
     Args:
         nqubits (int): Total number of qubits.
+        qubit_map (list, optional): the qubit map. If None, a list of range of circuit's qubits is used. Defaults to ``None``.
         noise_model (:class:`qibo.noise.NoiseModel`, optional): noise model used for simulating
             noisy computation. This matrix can be used to mitigate the effect of
             `qibo.noise.ReadoutError`.
-        nshots (int, optional): number of shots. Defaults to 10000.
+        nshots (int, optional): number of shots. Defaults to :math:`10000`.
         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``.
@@ -481,7 +485,7 @@ def get_calibration_matrix(
     if backend is None:  # pragma: no cover
         backend = GlobalBackend()
 
-    calibration_matrix = np.zeros((2**nqubits, 2**nqubits))
+    response_matrix = np.zeros((2**nqubits, 2**nqubits))
 
     for i in range(2**nqubits):
         binary_state = format(i, f"0{nqubits}b")
@@ -501,19 +505,19 @@ def get_calibration_matrix(
         column = np.zeros(2**nqubits)
         for key, value in frequencies.items():
             column[int(key, 2)] = value / nshots
-        calibration_matrix[:, i] = column
+        response_matrix[:, i] = column
 
-    return calibration_matrix
+    return response_matrix
 
 
-def apply_cal_mat_readout_mitigation(state, calibration_matrix, iterations=None):
+def apply_cal_mat_readout_mitigation(state, response_matrix, iterations=None):
     """
     Applies readout error mitigation to the given state using the provided calibration matrix.
 
     Args:
         state (:class:`qibo.measurements.CircuitResult`): the input state to be updated. This state should contain the
             frequencies that need to be mitigated.
-        calibration_matrix (numpy.ndarray, optional): the calibration matrix for readout mitigation.
+        response_matrix (numpy.ndarray): the calibration matrix for readout mitigation.
         iterations (int, optional): the number of iterations to use for the Iterative Bayesian Unfolding method.
             If ``None`` the 'inverse' method is used. Defaults to ``None``.
 
@@ -527,12 +531,12 @@ def apply_cal_mat_readout_mitigation(state, calibration_matrix, iterations=None)
     frequencies = frequencies.reshape(-1, 1)
 
     if iterations is None:
-        calibration_matrix = np.linalg.inv(calibration_matrix)
+        calibration_matrix = np.linalg.inv(response_matrix)
         for i, value in enumerate(calibration_matrix @ frequencies):
             state._frequencies[i] = float(value)
     else:
         mitigated_probabilities = iterative_bayesian_unfolding(
-            frequencies / np.sum(frequencies), calibration_matrix, iterations
+            frequencies / np.sum(frequencies), response_matrix, iterations
         )
         mitigated_frequencies = np.round(
             mitigated_probabilities * np.sum(frequencies), 0
@@ -550,7 +554,7 @@ def apply_cal_mat_readout_mitigation(state, calibration_matrix, iterations=None)
 def apply_randomized_readout_mitigation(
     circuit,
     noise_model=None,
-    nshots: int = int(1e3),
+    nshots: int = 10000,
     ncircuits: int = 10,
     qubit_map=None,
     backend=None,
@@ -561,9 +565,10 @@ def apply_randomized_readout_mitigation(
         circuit (:class:`qibo.models.Circuit`): input circuit.
         noise_model(:class:`qibo.noise.NoiseModel`, optional): noise model used for
             simulating noisy computation. Defaults to ``None``.
-        nshots (int, optional): number of shots. Defaults to 10000.
+        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``.
         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``.
@@ -633,11 +638,31 @@ def apply_readout_mitigation(
     circuit,
     observable,
     noise_model=None,
-    nshots: int = int(1e3),
+    nshots: int = 10000,
     readout={},
     qubit_map=None,
     backend=None,
 ):
+    """
+    Applies readout error mitigation to the given circuit and observable.
+
+    Args:
+        circuit (qibo.models.Circuit): input circuit.
+        observable (:class:`qibo.hamiltonians.Hamiltonian/:class:`qibo.hamiltonians.SymbolicHamiltonian`): The observable to be measured.
+        noise_model (qibo.models.noise.Noise, optional): the noise model to be applied. Defaults to ``None``.
+        nshots (int, optional): the number of shots for the circuit execution. Defaults to :math:`10000`.
+        readout (dict, optional): a dictionary that may contain the following keys:
+            - 'ncircuits': int, specifies the number of random circuits to use for the randomized method of readout error mitigation.
+            - 'calibration_matrix': numpy.ndarray, used for applying a pre-computed calibration matrix for readout error mitigation.
+            - 'ibu_iters': int, specifies the number of iterations for the iterative Bayesian unfolding method of readout error mitigation.
+            If provided, the corresponding readout error mitigation method is used. Defaults to {}.
+        qubit_map (list, optional): the qubit map. If None, a list of range of circuit's qubits is used. Defaults to ``None``.
+        backend (qibo.backends.abstract.Backend, optional): the backend to be used in the execution.
+            If None, it uses the global backend. Defaults to ``None``.
+
+    Returns:
+        float: The mitigated expectation value of the observable.
+    """
     if backend is None:  # pragma: no cover
         backend = GlobalBackend()
 
@@ -672,13 +697,7 @@ def sample_clifford_training_circuit(
     """Samples a training circuit for CDR by susbtituting some of the non-Clifford gates.
 
     Args:
-        circuit (:class:`qibo.models.Circuit`): circuit to sample from,
-            decomposed in ``RX(pi/2)``, ``X``, ``CNOT`` and ``RZ`` gates.
-        replacement_gates (list, optional): candidates for the substitution of the
-            non-Clifford gates. The ``list`` should be composed by ``tuples`` of the
-            form (``gates.XYZ``, ``kwargs``). For example, phase gates are used by default:
-            ``list((RZ, {'theta':0}), (RZ, {'theta':pi/2}), (RZ, {'theta':pi}), (RZ, {'theta':3*pi/2}))``.
-        sigma (float, optional): standard devation of the Gaussian distribution used for sampling.
+        circuit (:class:`qibo.models.Circuit`): circuit to sample from.
         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``.
@@ -732,8 +751,7 @@ def error_sensitive_circuit(circuit, observable, backend=None):
 
     Args:
         circuit (:class:`qibo.models.Circuit`): input circuit.
-        observable (numpy.ndarray): Pauli observable to be measured.
-        fuse (bool, optional): if True, fuse the single qubits gates in the circuit.
+        observable (:class:`qibo.hamiltonians.Hamiltonian/:class:`qibo.hamiltonians.SymbolicHamiltonian`): Pauli observable to be measured.
         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``.
@@ -824,16 +842,16 @@ def ICS(
 
     Args:
         circuit (:class:`qibo.models.Circuit`): input circuit.
-        observable (numpy.ndarray): the observable to be measured.
+        observable (:class:`qibo.hamiltonians.Hamiltonian/:class:`qibo.hamiltonians.SymbolicHamiltonian`): the observable to be measured.
         readout (dict, optional): A dictionary that may contain the following keys:
             - 'calibration_matrix': numpy.ndarray, used for applying a pre-computed calibration matrix for readout error mitigation.
             - 'random_ncircuits': int, specifies the number of random circuits to use for the randomized method of readout error mitigation.
-            - 'ibu_iters': int, specifies the number of iterations for the iterative Bayesian update method of readout error mitigation.
+            - 'ibu_iters': int, specifies the number of iterations for the iterative Bayesian unfolding method of readout error mitigation.
                             If ``None`` the 'inverse' method is used.
             If provided, the corresponding readout error mitigation method is used.
-        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``.
         noise_model (qibo.models.noise.Noise, optional): the noise model to be applied. Defaults to ``None``.
-        nshots (int, optional): the number of shots for the circuit execution. Defaults to 10000.
+        nshots (int, optional): the number of shots for the circuit execution. Defaults to :math:`10000`.
         n_training_samples (int, optional): the number of training samples. Defaults to 10.
         full_output (bool, optional): if ``True``, this function returns additional
             information: ``val``, ``optimal_params``, ``train_val``. Defaults to ``False``.
@@ -916,7 +934,21 @@ def ICS(
     return mitigated_expectation
 
 
-def _circuit_conf(circuit, qubit_map, noise_model, nshots, backend=None):
+def _circuit_conf(circuit, qubit_map, noise_model=None, nshots=10000, backend=None):
+    """
+    Helper function to execute the given circuit with the specified parameters.
+
+    Args:
+        circuit (qibo.models.Circuit): input circuit.
+        qubit_map (list): the qubit map. If ``None``, a list of range of circuit's qubits is used. Defaults to ``None``.
+        noise_model (qibo.models.noise.Noise, optional): The noise model to be applied. Defaults to ``None``.
+        nshots (int): the number of shots for the circuit execution. Defaults to :math:`10000`..
+        backend (qibo.backends.abstract.Backend, optional): the backend to be used in the execution.
+            If None, it uses the global backend. Defaults to ``None``.
+
+    Returns:
+        qibo.states.CircuitResult: The result of the circuit execution.
+    """
     from qibo.transpiler.placer import Custom
 
     if backend is None:  # pragma: no cover
diff --git a/tests/test_models_error_mitigation.py b/tests/test_models_error_mitigation.py
index 6951b4e6e0..d8c0066793 100644
--- a/tests/test_models_error_mitigation.py
+++ b/tests/test_models_error_mitigation.py
@@ -9,7 +9,7 @@
     ICS,
     ZNE,
     apply_readout_mitigation,
-    get_calibration_matrix,
+    get_response_matrix,
     sample_clifford_training_circuit,
     sample_training_circuit_cdr,
     vnCDR,
@@ -298,7 +298,7 @@ def test_readout_mitigation(backend, nqubits, nmeas, method, ibu_iters):
     noise = NoiseModel()
     noise.add(ReadoutError(probabilities=p), gate=gates.M)
     if method == "cal_matrix":
-        calibration = get_calibration_matrix(
+        calibration = get_response_matrix(
             nmeas, None, noise, nshots=nshots, backend=backend
         )
         readout = {"calibration_matrix": calibration, "ibu_iters": ibu_iters}

From d9d0e76cc18d9a77c2f6930d609543c6c2e8f6af Mon Sep 17 00:00:00 2001
From: Alejandro Sopena <alexsopena15@gmail.com>
Date: Wed, 20 Dec 2023 09:05:36 +0100
Subject: [PATCH 34/47] updating documentation

---
 doc/source/api-reference/qibo.rst   | 48 +++++++++++++++++++++----
 src/qibo/models/error_mitigation.py | 54 ++++++++++++++++-------------
 2 files changed, 71 insertions(+), 31 deletions(-)

diff --git a/doc/source/api-reference/qibo.rst b/doc/source/api-reference/qibo.rst
index 996586c4ab..20cbdfe037 100644
--- a/doc/source/api-reference/qibo.rst
+++ b/doc/source/api-reference/qibo.rst
@@ -279,7 +279,7 @@ and both can be used as standalone functions or in combination with the other
 general mitigation methods by setting the paramter `readout`.
 
 
-Calibration Matrix
+Response Matrix
 """"""""""""""""""
 Given :math:`n` qubits, all the possible :math:`2^n` states are constructed via the
 application of the corresponding sequence of :math:`X` gates
@@ -297,14 +297,22 @@ columns :math:`M=\big(F_0^{noisy},...,F_{n-1}^{noisy}\big)`. We have indeed that
 and, therefore, the calibration matrix obtained as :math:`M_{\text{cal}}=M^{-1}`
 can be used to recover the noise-free frequencies.
 
-.. autofunction:: qibo.models.error_mitigation.calibration_matrix
+The calibration matrix :math:`M_{\text{cal}}` lacks stochasticity, resulting in a 'negative probability' issue.
+The distributions that arise after applying :math:`M_{\text{cal}}` are quasiprobabilities;
+the individual elements can be negative surpass 1, provided they sum to 1.
+It is posible to use Iterative Bayesian Unfolding (IBU) to preserve non-negativity.
+See `Srinivasan et al <https://arxiv.org/abs/2210.12284>`_ for more details.
+
+
+
+.. autofunction:: qibo.models.error_mitigation.get_response_matrix
 
 
 .. autofunction:: qibo.models.error_mitigation.apply_readout_mitigation
 
 
-Randomized
-""""""""""
+Randomized readout mitigation
+""""""""""""""""""""""""""""""
 This approach converts the effect of any noise map :math:`A` into a single multiplication
 factor for each Pauli observable, that is, diagonalizes the measurement channel.
 The multiplication factor :math:`\lambda` can be directly measured even without
@@ -320,7 +328,7 @@ factor results in the mitigated Pauli expectation value :math:`\langle O\rangle_
 Zero Noise Extrapolation (ZNE)
 """"""""""""""""""""""""""""""
 
-Given a noisy circuit :math:`C` and an observable :math:`A`,  Zero Noise Extrapolation (ZNE)
+Given a noisy circuit :math:`C` and an observable :math:`A`, Zero Noise Extrapolation (ZNE)
 consists in running :math:`n+1` versions of the circuit with different noise levels
 :math:`\{c_j\}_{j=0..n}` and, for each of them, measuring the expected value of the observable
 :math:`E_j=\langle A\rangle_j`.
@@ -381,7 +389,7 @@ See `Sopena et al <https://arxiv.org/abs/2103.12680>`_ for more details.
 .. autofunction:: qibo.models.error_mitigation.CDR
 
 
-.. autofunction:: qibo.models.error_mitigation.sample_training_circuit
+.. autofunction:: qibo.models.error_mitigation.sample_training_circuit_cdr
 
 
 Variable Noise CDR (vnCDR)
@@ -416,6 +424,34 @@ See `Sopena et al <https://arxiv.org/abs/2103.12680>`_ for all the details.
 .. autofunction:: qibo.models.error_mitigation.vnCDR
 
 
+Importance Clifford Sampling (ICS)
+""""""""""""""""""""""""""""""
+
+In the Importance Clifford Sampling (ICS) method, a set of :math:`n` circuits
+:math:`S_n=\{C_i\}_{i=1,..,n}` that stabilizes a given Pauli observable is generated starting from the original circuit
+:math:`C_0` by replacing all the non-Clifford gates with Clifford ones.
+Given an observable :math:`A`, all the circuits of :math:`S_n` are both simulated
+to obtain the correspondent expected values of :math:`A` in noise-free condition
+:math:`\{a_i^{exact}\}_{i=1,..,n}`, and run in noisy conditions to obtain the noisy
+expected values :math:`\{a_i^{noisy}\}_{i=1,..,n}`.
+
+Finally, a theoretically inspired model :math:`f` is learned using the training data.
+
+The mitigated expected value of :math:`A` at the end of :math:`C_0` is then
+obtained simply with :math:`f(a_0^{noisy})`.
+
+In this implementation the initial circuit is expected to be decomposed in the three
+Clifford gates :math:`RX(\frac{\pi}{2})`, :math:`CNOT`, :math:`X` and in :math:`RZ(\theta)`
+(which is Clifford only for :math:`\theta=\frac{n\pi}{2}`).
+By default the set of Clifford gates used for substitution is
+:math:`\{RZ(0),RZ(\frac{\pi}{2}),RZ(\pi),RZ(\frac{3}{2}\pi)\}`.
+See `Sopena et al <https://arxiv.org/abs/2103.12680>`_ for more details.
+
+.. autofunction:: qibo.models.error_mitigation.ICS
+
+
+.. autofunction:: qibo.models.error_mitigation.sample_clifford_training_circuit
+
 _______________________
 
 .. _Gates:
diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index 93792fac49..6d1bb4cf64 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -130,11 +130,12 @@ def ZNE(
         insertion_gate (str, optional): gate to be used in the insertion.
             If ``"RX"``, the gate used is :math:``RX(\\pi / 2)``.
             Defaults to ``"CNOT"``.
-        readout (dict, optional): A dictionary that may contain the following keys:
-            - 'calibration_matrix': numpy.ndarray, used for applying a pre-computed calibration matrix for readout error mitigation.
-            - 'random_ncircuits': int, specifies the number of random circuits to use for the randomized method of readout error mitigation.
-            - 'ibu_iters': int, specifies the number of iterations for the iterative Bayesian unfolding method of readout error mitigation.
-            If provided, the corresponding readout error mitigation method is used.
+        readout (dict, optional): a dictionary that may contain the following keys:
+
+            *    ncircuits: int, specifies the number of random circuits to use for the randomized method of readout error mitigation.
+            *    calibration_matrix: numpy.ndarray, used for applying a pre-computed calibration matrix for readout error mitigation.
+            *    ibu_iters: int, specifies the number of iterations for the iterative Bayesian unfolding method of readout error mitigation. If provided, the corresponding readout error mitigation method is used. Defaults to {}.
+
         qubit_map (list, optional): the qubit map. If None, a list of range of circuit's qubits is used. Defaults to ``None``.
         backend (:class:`qibo.backends.abstract.Backend`, optional): backend to be used
             in the execution. If ``None``, it uses :class:`qibo.backends.GlobalBackend`.
@@ -279,11 +280,12 @@ def CDR(
         n_training_samples (int, optional): number of training circuits to sample. Defaults to 100.
         full_output (bool, optional): if ``True``, this function returns additional
             information: ``val``, ``optimal_params``, ``train_val``. Defaults to ``False``.
-        readout (dict, optional): A dictionary that may contain the following keys:
-            - 'calibration_matrix': numpy.ndarray, used for applying a pre-computed calibration matrix for readout error mitigation.
-            - 'random_ncircuits': int, specifies the number of random circuits to use for the randomized method of readout error mitigation.
-            - 'ibu_iters': int, specifies the number of iterations for the iterative Bayesian unfolding method of readout error mitigation.
-            If provided, the corresponding readout error mitigation method is used.
+        readout (dict, optional): a dictionary that may contain the following keys:
+
+            *    ncircuits: int, specifies the number of random circuits to use for the randomized method of readout error mitigation.
+            *    calibration_matrix: numpy.ndarray, used for applying a pre-computed calibration matrix for readout error mitigation.
+            *    ibu_iters: int, specifies the number of iterations for the iterative Bayesian unfolding method of readout error mitigation. If provided, the corresponding readout error mitigation method is used. Defaults to {}.
+
         qubit_map (list, optional): the qubit map. If None, a list of range of circuit's qubits is used. Defaults to ``None``.
         backend (:class:`qibo.backends.abstract.Backend`, optional): backend to be used
             in the execution. If ``None``, it uses :class:`qibo.backends.GlobalBackend`.
@@ -361,11 +363,12 @@ def vnCDR(
             Default is ``"CNOT"``.
         full_output (bool, optional): if ``True``, this function returns additional
             information: ``val``, ``optimal_params``, ``train_val``. Defaults to ``False``.
-        readout (dict, optional): A dictionary that may contain the following keys:
-            - 'calibration_matrix': numpy.ndarray, used for applying a pre-computed calibration matrix for readout error mitigation.
-            - 'random_ncircuits': int, specifies the number of random circuits to use for the randomized method of readout error mitigation.
-            - 'ibu_iters': int, specifies the number of iterations for the iterative Bayesian unfolding method of readout error mitigation.
-            If provided, the corresponding readout error mitigation method is used.
+        readout (dict, optional): a dictionary that may contain the following keys:
+
+            *    ncircuits: int, specifies the number of random circuits to use for the randomized method of readout error mitigation.
+            *    calibration_matrix: numpy.ndarray, used for applying a pre-computed calibration matrix for readout error mitigation.
+            *    ibu_iters: int, specifies the number of iterations for the iterative Bayesian unfolding method of readout error mitigation. If provided, the corresponding readout error mitigation method is used. Defaults to {}.
+
         qubit_map (list, optional): the qubit map. If None, a list of range of circuit's qubits is used. Defaults to ``None``.
         backend (:class:`qibo.backends.abstract.Backend`, optional): backend to be used
             in the execution. If ``None``, it uses :class:`qibo.backends.GlobalBackend`.
@@ -652,10 +655,11 @@ def apply_readout_mitigation(
         noise_model (qibo.models.noise.Noise, optional): the noise model to be applied. Defaults to ``None``.
         nshots (int, optional): the number of shots for the circuit execution. Defaults to :math:`10000`.
         readout (dict, optional): a dictionary that may contain the following keys:
-            - 'ncircuits': int, specifies the number of random circuits to use for the randomized method of readout error mitigation.
-            - 'calibration_matrix': numpy.ndarray, used for applying a pre-computed calibration matrix for readout error mitigation.
-            - 'ibu_iters': int, specifies the number of iterations for the iterative Bayesian unfolding method of readout error mitigation.
-            If provided, the corresponding readout error mitigation method is used. Defaults to {}.
+
+            *    ncircuits: int, specifies the number of random circuits to use for the randomized method of readout error mitigation.
+            *    calibration_matrix: numpy.ndarray, used for applying a pre-computed calibration matrix for readout error mitigation.
+            *    ibu_iters: int, specifies the number of iterations for the iterative Bayesian unfolding method of readout error mitigation. If provided, the corresponding readout error mitigation method is used. Defaults to {}.
+
         qubit_map (list, optional): the qubit map. If None, a list of range of circuit's qubits is used. Defaults to ``None``.
         backend (qibo.backends.abstract.Backend, optional): the backend to be used in the execution.
             If None, it uses the global backend. Defaults to ``None``.
@@ -843,12 +847,12 @@ def ICS(
     Args:
         circuit (:class:`qibo.models.Circuit`): input circuit.
         observable (:class:`qibo.hamiltonians.Hamiltonian/:class:`qibo.hamiltonians.SymbolicHamiltonian`): the observable to be measured.
-        readout (dict, optional): A dictionary that may contain the following keys:
-            - 'calibration_matrix': numpy.ndarray, used for applying a pre-computed calibration matrix for readout error mitigation.
-            - 'random_ncircuits': int, specifies the number of random circuits to use for the randomized method of readout error mitigation.
-            - 'ibu_iters': int, specifies the number of iterations for the iterative Bayesian unfolding method of readout error mitigation.
-                            If ``None`` the 'inverse' method is used.
-            If provided, the corresponding readout error mitigation method is used.
+        readout (dict, optional): a dictionary that may contain the following keys:
+
+            *    ncircuits: int, specifies the number of random circuits to use for the randomized method of readout error mitigation.
+            *    calibration_matrix: numpy.ndarray, used for applying a pre-computed calibration matrix for readout error mitigation.
+            *    ibu_iters: int, specifies the number of iterations for the iterative Bayesian unfolding method of readout error mitigation. If provided, the corresponding readout error mitigation method is used. Defaults to {}.
+
         qubit_map (list, optional): the qubit map. If ``None``, a list of range of circuit's qubits is used. Defaults to ``None``.
         noise_model (qibo.models.noise.Noise, optional): the noise model to be applied. Defaults to ``None``.
         nshots (int, optional): the number of shots for the circuit execution. Defaults to :math:`10000`.

From caf47847235772b90ea552077feec620c8489458 Mon Sep 17 00:00:00 2001
From: Alejandro Sopena <alexsopena15@gmail.com>
Date: Wed, 20 Dec 2023 09:28:23 +0100
Subject: [PATCH 35/47] fix example

---
 doc/source/code-examples/advancedexamples.rst | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/doc/source/code-examples/advancedexamples.rst b/doc/source/code-examples/advancedexamples.rst
index 4a3fd766f8..53c166ca07 100644
--- a/doc/source/code-examples/advancedexamples.rst
+++ b/doc/source/code-examples/advancedexamples.rst
@@ -1293,11 +1293,11 @@ calibration matrix and use it modify the final state after the circuit execution
 
 .. testcode::
 
-   from qibo.models.error_mitigation import apply_cal_mat_readout_mitigation, get_calibration_matrix
+   from qibo.models.error_mitigation import apply_cal_mat_readout_mitigation, get_response_matrix
 
    nshots = 10000
    # compute the calibration matrix
-   calibration = get_calibration_matrix(
+   calibration = get_response_matrix(
        nqubits, backend=backend, noise_model=noise, nshots=nshots
    )
    # execute the circuit

From 4af1eabb1a9c57fe592128b183abde3db1f76a8a Mon Sep 17 00:00:00 2001
From: Alejandro Sopena <alexsopena15@gmail.com>
Date: Wed, 20 Dec 2023 09:52:49 +0100
Subject: [PATCH 36/47] add example for ICS method

---
 doc/source/code-examples/advancedexamples.rst | 103 +++++++++++-------
 1 file changed, 65 insertions(+), 38 deletions(-)

diff --git a/doc/source/code-examples/advancedexamples.rst b/doc/source/code-examples/advancedexamples.rst
index 53c166ca07..cbea186c4d 100644
--- a/doc/source/code-examples/advancedexamples.rst
+++ b/doc/source/code-examples/advancedexamples.rst
@@ -1205,23 +1205,23 @@ Let's see how to use them. For starters, let's define a dummy circuit with some
    hz = 0.5
    hx = 0.5
    dt = 0.25
-   c = Circuit(nqubits, density_matrix=True)
-   c.add(gates.RZ(q, theta=-2 * hz * dt - np.pi / 2) for q in range(nqubits))
-   c.add(gates.RX(q, theta=np.pi / 2) for q in range(nqubits))
-   c.add(gates.RZ(q, theta=-2 * hx * dt + np.pi) for q in range(nqubits))
-   c.add(gates.RX(q, theta=np.pi / 2) for q in range(nqubits))
-   c.add(gates.RZ(q, theta=-np.pi / 2) for q in range(nqubits))
-   c.add(gates.CNOT(q, q + 1) for q in range(0, nqubits - 1, 2))
-   c.add(gates.RZ(q + 1, theta=-2 * dt) for q in range(0, nqubits - 1, 2))
-   c.add(gates.CNOT(q, q + 1) for q in range(0, nqubits - 1, 2))
-   c.add(gates.CNOT(q, q + 1) for q in range(1, nqubits, 2))
-   c.add(gates.RZ(q + 1, theta=-2 * dt) for q in range(1, nqubits, 2))
-   c.add(gates.CNOT(q, q + 1) for q in range(1, nqubits, 2))
+   circ = Circuit(nqubits, density_matrix=True)
+   circ.add(gates.RZ(q, theta=-2 * hz * dt - np.pi / 2) for q in range(nqubits))
+   circ.add(gates.RX(q, theta=np.pi / 2) for q in range(nqubits))
+   circ.add(gates.RZ(q, theta=-2 * hx * dt + np.pi) for q in range(nqubits))
+   circ.add(gates.RX(q, theta=np.pi / 2) for q in range(nqubits))
+   circ.add(gates.RZ(q, theta=-np.pi / 2) for q in range(nqubits))
+   circ.add(gates.CNOT(q, q + 1) for q in range(0, nqubits - 1, 2))
+   circ.add(gates.RZ(q + 1, theta=-2 * dt) for q in range(0, nqubits - 1, 2))
+   circ.add(gates.CNOT(q, q + 1) for q in range(0, nqubits - 1, 2))
+   circ.add(gates.CNOT(q, q + 1) for q in range(1, nqubits, 2))
+   circ.add(gates.RZ(q + 1, theta=-2 * dt) for q in range(1, nqubits, 2))
+   circ.add(gates.CNOT(q, q + 1) for q in range(1, nqubits, 2))
    # Include the measurements
-   c.add(gates.M(*range(nqubits)))
+   circ.add(gates.M(*range(nqubits)))
 
    # visualize the circuit
-   print(c.draw())
+   print(circ.draw())
 
    #  q0: ─RZ─RX─RZ─RX─RZ─o────o────────M─
    #  q1: ─RZ─RX─RZ─RX─RZ─X─RZ─X─o────o─M─
@@ -1254,7 +1254,7 @@ the real quantum hardware, instead, we can use a noise model:
 .. testcode::
 
    # Noise-free expected value
-   exact = obs.expectation(backend.execute_circuit(c).state())
+   exact = obs.expectation(backend.execute_circuit(circ).state())
    print(exact)
    # 0.9096065335014379
 
@@ -1272,7 +1272,7 @@ the real quantum hardware, instead, we can use a noise model:
    )
    noise.add(ReadoutError(probabilities=prob), gate=gates.M)
    # Noisy expected value without mitigation
-   noisy = obs.expectation(backend.execute_circuit(noise.apply(c)).state())
+   noisy = obs.expectation(backend.execute_circuit(noise.apply(circ)).state())
    print(noisy)
    # 0.5647937721701448
 
@@ -1293,18 +1293,17 @@ calibration matrix and use it modify the final state after the circuit execution
 
 .. testcode::
 
-   from qibo.models.error_mitigation import apply_cal_mat_readout_mitigation, get_response_matrix
+   from qibo.models.error_mitigation import apply_readout_mitigation, get_response_matrix
 
    nshots = 10000
-   # compute the calibration matrix
-   calibration = get_response_matrix(
+   # compute the response matrix
+   response_matrix = get_response_matrix(
        nqubits, backend=backend, noise_model=noise, nshots=nshots
    )
-   # execute the circuit
-   state = backend.execute_circuit(noise.apply(c), nshots=nshots)
+   # define mitigation options
+   readout = {"calibration_matrix": response_matrix}
    # mitigate the readout errors
-   mit_state = apply_cal_mat_readout_mitigation(state, calibration)
-   mit_val = mit_state.expectation_from_samples(obs)
+   mit_val = apply_readout_mitigation(circ, obs, noise, readout=readout)
    print(mit_val)
    # 0.5945794816381054
 
@@ -1319,14 +1318,10 @@ Or use the randomized readout mitigation:
 
    from qibo.models.error_mitigation import apply_randomized_readout_mitigation
 
-   ncircuits = 10
-   result, result_cal = apply_randomized_readout_mitigation(
-       c, backend=backend, noise_model=noise, nshots=nshots, ncircuits=ncircuits
-   )
-   mit_val = result.expectation_from_samples(
-       obs
-   ) / result_cal.expectation_from_samples(obs)
-   print(mit_val)
+   # define mitigation options
+   readout = {"ncircuits": 10}
+   # mitigate the readout errors
+   mit_val = apply_readout_mitigation(circ, obs, noise, readout=readout)
    # 0.5860884499785314
 
 .. testoutput::
@@ -1365,7 +1360,7 @@ For example if we use the five levels ``[0,1,2,3,4]`` :
 
    # Mitigated expected value
    estimate = ZNE(
-       circuit=c,
+       circuit=circ,
        observable=obs,
        noise_levels=np.arange(5),
        noise_model=noise,
@@ -1388,13 +1383,13 @@ combined with the readout mitigation:
 
    # we can either use
    # the calibration matrix computed earlier
-   readout = {'calibration_matrix': calibration}
+   readout = {'calibration_matrix': response_matrix}
    # or the randomized readout
    readout = {'ncircuits': 10}
 
    # Mitigated expected value
    estimate = ZNE(
-       circuit=c,
+       circuit=circ,
        observable=obs,
        backend=backend,
        noise_levels=np.arange(5),
@@ -1424,15 +1419,16 @@ circuit is expected to be decomposed in the set of primitive gates :math:`RX(\fr
 
    # Mitigated expected value
    estimate = CDR(
-       circuit=c,
+       circuit=circ,
        observable=obs,
+       n_training_samples=10,
        backend=backend,
        noise_model=noise,
        nshots=10000,
        readout=readout,
    )
    print(estimate)
-   # 0.9090604794014961
+   # 0.8983676333969615
 
 .. testoutput::
    :hide:
@@ -1441,6 +1437,7 @@ circuit is expected to be decomposed in the set of primitive gates :math:`RX(\fr
 
 Again, the mitigated expected value improves over the noisy one and is also slightly better compared to ZNE.
 
+
 Variable Noise CDR (vnCDR)
 ^^^^^^^^^^^^^^^^^^^^^^^^^^
 
@@ -1453,8 +1450,9 @@ caveat about the input circuit for CDR is valid here as well.
 
    # Mitigated expected value
    estimate = vnCDR(
-       circuit=c,
+       circuit=circ,
        observable=obs,
+       n_training_samples=10,
        backend=backend,
        noise_levels=np.arange(3),
        noise_model=noise,
@@ -1463,7 +1461,7 @@ caveat about the input circuit for CDR is valid here as well.
        readout=readout,
    )
    print(estimate)
-   # 0.9085991439303123
+   # 0.8998376314644383
 
 .. testoutput::
    :hide:
@@ -1474,6 +1472,35 @@ The result is similar to the one obtained by CDR. Usually, one would expect slig
 however, this can substantially vary depending on the circuit and the observable considered and, therefore, it is hard to tell
 a priori.
 
+
+Importance Clifford Sampling (ICS)
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+The use of iCS is straightforward, analogous to CDR and vnCDR.
+
+.. testcode::
+
+   from qibo.models.error_mitigation import ICS
+
+   # Mitigated expected value
+   estimate = ICS(
+       circuit=circ,
+       observable=obs,
+       n_training_samples=10,
+       backend=backend,
+       noise_model=noise,
+       nshots=10000,
+       readout=readout,
+   )
+   print(estimate)
+   # 0.9183495097398502
+
+.. testoutput::
+   :hide:
+
+   ...
+
+Again, the mitigated expected value improves over the noisy one and is also slightly better compared to ZNE.
 This was just a basic example usage of the three methods, for all the details about them you should check the API-reference page :ref:`Error Mitigation <error-mitigation>`.
 
 .. _timeevol-example:

From 2b39afebb5c8ad1800caa10fb89dfba7a9833744 Mon Sep 17 00:00:00 2001
From: Alejandro Sopena <alexsopena15@gmail.com>
Date: Wed, 20 Dec 2023 10:51:28 +0100
Subject: [PATCH 37/47] adding print to example

---
 doc/source/code-examples/advancedexamples.rst | 1 +
 1 file changed, 1 insertion(+)

diff --git a/doc/source/code-examples/advancedexamples.rst b/doc/source/code-examples/advancedexamples.rst
index cbea186c4d..87f7397c24 100644
--- a/doc/source/code-examples/advancedexamples.rst
+++ b/doc/source/code-examples/advancedexamples.rst
@@ -1322,6 +1322,7 @@ Or use the randomized readout mitigation:
    readout = {"ncircuits": 10}
    # mitigate the readout errors
    mit_val = apply_readout_mitigation(circ, obs, noise, readout=readout)
+   print(mit_val)
    # 0.5860884499785314
 
 .. testoutput::

From 6fc230db703742a0436e46d5deae4525cb5a6806 Mon Sep 17 00:00:00 2001
From: Alejandro Sopena <alexsopena15@gmail.com>
Date: Wed, 20 Dec 2023 12:10:52 +0100
Subject: [PATCH 38/47] Add reference

---
 doc/source/api-reference/qibo.rst   | 11 ++++++++++-
 src/qibo/models/error_mitigation.py |  4 +++-
 2 files changed, 13 insertions(+), 2 deletions(-)

diff --git a/doc/source/api-reference/qibo.rst b/doc/source/api-reference/qibo.rst
index 20cbdfe037..c472331424 100644
--- a/doc/source/api-reference/qibo.rst
+++ b/doc/source/api-reference/qibo.rst
@@ -301,13 +301,22 @@ The calibration matrix :math:`M_{\text{cal}}` lacks stochasticity, resulting in
 The distributions that arise after applying :math:`M_{\text{cal}}` are quasiprobabilities;
 the individual elements can be negative surpass 1, provided they sum to 1.
 It is posible to use Iterative Bayesian Unfolding (IBU) to preserve non-negativity.
-See `Srinivasan et al <https://arxiv.org/abs/2210.12284>`_ for more details.
+See `Nachman et al <https://arxiv.org/abs/1910.01969>`_ for more details.
 
 
 
 .. autofunction:: qibo.models.error_mitigation.get_response_matrix
 
 
+.. autofunction:: qibo.models.error_mitigation.iterative_bayesian_unfolding
+
+
+.. autofunction:: qibo.models.error_mitigation.apply_cal_mat_readout_mitigation
+
+
+.. autofunction:: qibo.models.error_mitigation.apply_randomized_readout_mitigation
+
+
 .. autofunction:: qibo.models.error_mitigation.apply_readout_mitigation
 
 
diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index 6d1bb4cf64..fdf5ffd5d9 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -449,7 +449,9 @@ def iterative_bayesian_unfolding(probabilities, response_matrix, iterations=10):
         numpy.ndarray: the unfolded probabilities.
 
     Reference:
-        1. S. Srinivasan, B. Pokharel et al, *Scalable Measurement Error Mitigation via Iterative Bayesian Unfolding*.
+        1. B. Nachman, M. Urbanek et al, *Unfolding Quantum Computer Readout Noise*.
+           `arXiv:1910.01969 [quant-ph] <https://arxiv.org/abs/1910.01969>`_.
+        2. S. Srinivasan, B. Pokharel et al, *Scalable Measurement Error Mitigation via Iterative Bayesian Unfolding*.
            `arXiv:2210.12284 [quant-ph] <https://arxiv.org/abs/2210.12284>`_.
     """
     unfolded_probabilities = np.ones((len(probabilities), 1)) / len(probabilities)

From 93024721011020752d34d0e4a8e34ff7e58abe91 Mon Sep 17 00:00:00 2001
From: Alejandro Sopena <alexsopena15@gmail.com>
Date: Sun, 7 Jan 2024 03:11:20 +0100
Subject: [PATCH 39/47] matteo's suggestions

---
 doc/source/api-reference/qibo.rst             |  8 +--
 doc/source/code-examples/advancedexamples.rst | 14 ++---
 src/qibo/models/error_mitigation.py           | 42 +++++++-------
 tests/test_models_error_mitigation.py         | 58 ++++++++++---------
 4 files changed, 62 insertions(+), 60 deletions(-)

diff --git a/doc/source/api-reference/qibo.rst b/doc/source/api-reference/qibo.rst
index dc6122373f..3725e143b9 100644
--- a/doc/source/api-reference/qibo.rst
+++ b/doc/source/api-reference/qibo.rst
@@ -288,7 +288,7 @@ In the presence of readout errors, we will measure for each state :math:`i` some
 frequencies :math:`F_i^{noisy}` different from the ideal ones
 :math:`F_i^{ideal}=\delta_{i,j}`.
 
-The effect of the error is modeled by the matrix composed of the noisy frequencies as
+The effect of the error is modeled by the response matrix composed of the noisy frequencies as
 columns :math:`M=\big(F_0^{noisy},...,F_{n-1}^{noisy}\big)`. We have indeed that:
 
 .. math::
@@ -311,13 +311,13 @@ See `Nachman et al <https://arxiv.org/abs/1910.01969>`_ for more details.
 .. autofunction:: qibo.models.error_mitigation.iterative_bayesian_unfolding
 
 
-.. autofunction:: qibo.models.error_mitigation.apply_cal_mat_readout_mitigation
+.. autofunction:: qibo.models.error_mitigation.apply_resp_mat_readout_mitigation
 
 
 .. autofunction:: qibo.models.error_mitigation.apply_randomized_readout_mitigation
 
 
-.. autofunction:: qibo.models.error_mitigation.apply_readout_mitigation
+.. autofunction:: qibo.models.error_mitigation.apply_problem_with_readout_conf
 
 
 Randomized readout mitigation
@@ -434,7 +434,7 @@ See `Sopena et al <https://arxiv.org/abs/2103.12680>`_ for all the details.
 
 
 Importance Clifford Sampling (ICS)
-""""""""""""""""""""""""""""""
+""""""""""""""""""""""""""""""""""
 
 In the Importance Clifford Sampling (ICS) method, a set of :math:`n` circuits
 :math:`S_n=\{C_i\}_{i=1,..,n}` that stabilizes a given Pauli observable is generated starting from the original circuit
diff --git a/doc/source/code-examples/advancedexamples.rst b/doc/source/code-examples/advancedexamples.rst
index 87f7397c24..73176acd65 100644
--- a/doc/source/code-examples/advancedexamples.rst
+++ b/doc/source/code-examples/advancedexamples.rst
@@ -1289,11 +1289,11 @@ Now let's check that error mitigation produces better estimates of the exact exp
 Readout Mitigation
 ^^^^^^^^^^^^^^^^^^
 Firstly, let's try to mitigate the readout errors. To do this, we can either compute the
-calibration matrix and use it modify the final state after the circuit execution:
+response matrix and use it modify the final state after the circuit execution:
 
 .. testcode::
 
-   from qibo.models.error_mitigation import apply_readout_mitigation, get_response_matrix
+   from qibo.models.error_mitigation import apply_problem_with_readout_conf, get_response_matrix
 
    nshots = 10000
    # compute the response matrix
@@ -1301,9 +1301,9 @@ calibration matrix and use it modify the final state after the circuit execution
        nqubits, backend=backend, noise_model=noise, nshots=nshots
    )
    # define mitigation options
-   readout = {"calibration_matrix": response_matrix}
+   readout = {"response_matrix": response_matrix}
    # mitigate the readout errors
-   mit_val = apply_readout_mitigation(circ, obs, noise, readout=readout)
+   mit_val = apply_problem_with_readout_conf(circ, obs, noise, readout=readout)
    print(mit_val)
    # 0.5945794816381054
 
@@ -1321,7 +1321,7 @@ Or use the randomized readout mitigation:
    # define mitigation options
    readout = {"ncircuits": 10}
    # mitigate the readout errors
-   mit_val = apply_readout_mitigation(circ, obs, noise, readout=readout)
+   mit_val = apply_problem_with_readout_conf(circ, obs, noise, readout=readout)
    print(mit_val)
    # 0.5860884499785314
 
@@ -1383,8 +1383,8 @@ combined with the readout mitigation:
 .. testcode::
 
    # we can either use
-   # the calibration matrix computed earlier
-   readout = {'calibration_matrix': response_matrix}
+   # the response matrix computed earlier
+   readout = {'response_matrix': response_matrix}
    # or the randomized readout
    readout = {'ncircuits': 10}
 
diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index fdf5ffd5d9..5f5f1c44ca 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -133,7 +133,7 @@ def ZNE(
         readout (dict, optional): a dictionary that may contain the following keys:
 
             *    ncircuits: int, specifies the number of random circuits to use for the randomized method of readout error mitigation.
-            *    calibration_matrix: numpy.ndarray, used for applying a pre-computed calibration matrix for readout error mitigation.
+            *    response_matrix: numpy.ndarray, used for applying a pre-computed response matrix for readout error mitigation.
             *    ibu_iters: int, specifies the number of iterations for the iterative Bayesian unfolding method of readout error mitigation. If provided, the corresponding readout error mitigation method is used. Defaults to {}.
 
         qubit_map (list, optional): the qubit map. If None, a list of range of circuit's qubits is used. Defaults to ``None``.
@@ -156,7 +156,7 @@ def ZNE(
         noisy_circuit = get_noisy_circuit(
             circuit, num_insertions, insertion_gate=insertion_gate
         )
-        val = apply_readout_mitigation(
+        val = apply_problem_with_readout_conf(
             noisy_circuit,
             observable,
             noise_model,
@@ -283,7 +283,7 @@ def CDR(
         readout (dict, optional): a dictionary that may contain the following keys:
 
             *    ncircuits: int, specifies the number of random circuits to use for the randomized method of readout error mitigation.
-            *    calibration_matrix: numpy.ndarray, used for applying a pre-computed calibration matrix for readout error mitigation.
+            *    response_matrix: numpy.ndarray, used for applying a pre-computed response matrix for readout error mitigation.
             *    ibu_iters: int, specifies the number of iterations for the iterative Bayesian unfolding method of readout error mitigation. If provided, the corresponding readout error mitigation method is used. Defaults to {}.
 
         qubit_map (list, optional): the qubit map. If None, a list of range of circuit's qubits is used. Defaults to ``None``.
@@ -312,14 +312,14 @@ def CDR(
     for circ in training_circuits:
         val = circ(nshots=nshots).expectation_from_samples(observable)
         train_val["noise-free"].append(val)
-        val = apply_readout_mitigation(
+        val = apply_problem_with_readout_conf(
             circ, observable, noise_model, nshots, readout, qubit_map, backend=backend
         )
         train_val["noisy"].append(val)
 
     optimal_params = curve_fit(model, train_val["noisy"], train_val["noise-free"])[0]
 
-    val = apply_readout_mitigation(
+    val = apply_problem_with_readout_conf(
         circuit, observable, noise_model, nshots, readout, qubit_map, backend=backend
     )
     mit_val = model(val, *optimal_params)
@@ -366,7 +366,7 @@ def vnCDR(
         readout (dict, optional): a dictionary that may contain the following keys:
 
             *    ncircuits: int, specifies the number of random circuits to use for the randomized method of readout error mitigation.
-            *    calibration_matrix: numpy.ndarray, used for applying a pre-computed calibration matrix for readout error mitigation.
+            *    response_matrix: numpy.ndarray, used for applying a pre-computed response matrix for readout error mitigation.
             *    ibu_iters: int, specifies the number of iterations for the iterative Bayesian unfolding method of readout error mitigation. If provided, the corresponding readout error mitigation method is used. Defaults to {}.
 
         qubit_map (list, optional): the qubit map. If None, a list of range of circuit's qubits is used. Defaults to ``None``.
@@ -398,7 +398,7 @@ def vnCDR(
         train_val["noise-free"].append(val)
         for level in noise_levels:
             noisy_c = get_noisy_circuit(circ, level, insertion_gate=insertion_gate)
-            val = apply_readout_mitigation(
+            val = apply_problem_with_readout_conf(
                 noisy_c,
                 observable,
                 noise_model,
@@ -417,7 +417,7 @@ def vnCDR(
     val = []
     for level in noise_levels:
         noisy_c = get_noisy_circuit(circuit, level, insertion_gate=insertion_gate)
-        expval = apply_readout_mitigation(
+        expval = apply_problem_with_readout_conf(
             noisy_c,
             observable,
             noise_model,
@@ -468,7 +468,7 @@ def iterative_bayesian_unfolding(probabilities, response_matrix, iterations=10):
 def get_response_matrix(
     nqubits, qubit_map=None, noise_model=None, nshots: int = 10000, backend=None
 ):
-    """Computes the calibration matrix for readout mitigation.
+    """Computes the response matrix for readout mitigation.
 
     Args:
         nqubits (int): Total number of qubits.
@@ -482,7 +482,7 @@ def get_response_matrix(
             Defaults to ``None``.
 
     Returns:
-        numpy.ndarray : The computed (`nqubits`, `nqubits`) calibration matrix for
+        numpy.ndarray : The computed (`nqubits`, `nqubits`) response matrix for
             readout mitigation.
     """
     from qibo import Circuit  # pylint: disable=import-outside-toplevel
@@ -515,14 +515,14 @@ def get_response_matrix(
     return response_matrix
 
 
-def apply_cal_mat_readout_mitigation(state, response_matrix, iterations=None):
+def apply_resp_mat_readout_mitigation(state, response_matrix, iterations=None):
     """
-    Applies readout error mitigation to the given state using the provided calibration matrix.
+    Applies readout error mitigation to the given state using the provided response matrix.
 
     Args:
         state (:class:`qibo.measurements.CircuitResult`): the input state to be updated. This state should contain the
             frequencies that need to be mitigated.
-        response_matrix (numpy.ndarray): the calibration matrix for readout mitigation.
+        response_matrix (numpy.ndarray): the response matrix for readout mitigation.
         iterations (int, optional): the number of iterations to use for the Iterative Bayesian Unfolding method.
             If ``None`` the 'inverse' method is used. Defaults to ``None``.
 
@@ -639,7 +639,7 @@ def apply_randomized_readout_mitigation(
     return results
 
 
-def apply_readout_mitigation(
+def apply_problem_with_readout_conf(
     circuit,
     observable,
     noise_model=None,
@@ -659,7 +659,7 @@ def apply_readout_mitigation(
         readout (dict, optional): a dictionary that may contain the following keys:
 
             *    ncircuits: int, specifies the number of random circuits to use for the randomized method of readout error mitigation.
-            *    calibration_matrix: numpy.ndarray, used for applying a pre-computed calibration matrix for readout error mitigation.
+            *    response_matrix: numpy.ndarray, used for applying a pre-computed response matrix for readout error mitigation.
             *    ibu_iters: int, specifies the number of iterations for the iterative Bayesian unfolding method of readout error mitigation. If provided, the corresponding readout error mitigation method is used. Defaults to {}.
 
         qubit_map (list, optional): the qubit map. If None, a list of range of circuit's qubits is used. Defaults to ``None``.
@@ -681,10 +681,10 @@ def apply_readout_mitigation(
             circuit, qubit_map, noise_model, nshots, backend=backend
         )
 
-    if "calibration_matrix" in readout:
-        circuit_result = apply_cal_mat_readout_mitigation(
+    if "response_matrix" in readout:
+        circuit_result = apply_resp_mat_readout_mitigation(
             circuit_result,
-            readout["calibration_matrix"],
+            readout["response_matrix"],
             readout.get("ibu_iters", None),
         )
 
@@ -852,7 +852,7 @@ def ICS(
         readout (dict, optional): a dictionary that may contain the following keys:
 
             *    ncircuits: int, specifies the number of random circuits to use for the randomized method of readout error mitigation.
-            *    calibration_matrix: numpy.ndarray, used for applying a pre-computed calibration matrix for readout error mitigation.
+            *    response_matrix: numpy.ndarray, used for applying a pre-computed response matrix for readout error mitigation.
             *    ibu_iters: int, specifies the number of iterations for the iterative Bayesian unfolding method of readout error mitigation. If provided, the corresponding readout error mitigation method is used. Defaults to {}.
 
         qubit_map (list, optional): the qubit map. If ``None``, a list of range of circuit's qubits is used. Defaults to ``None``.
@@ -894,7 +894,7 @@ def ICS(
         circuit_result = training_circuit(nshots=nshots)
         expectation = observable.expectation_from_samples(circuit_result.frequencies())
 
-        noisy_expectation = apply_readout_mitigation(
+        noisy_expectation = apply_problem_with_readout_conf(
             training_circuit,
             observable,
             noise_model,
@@ -911,7 +911,7 @@ def ICS(
     dep_param = np.mean(lambda_list)
     dep_param_std = np.std(lambda_list)
 
-    noisy_expectation = apply_readout_mitigation(
+    noisy_expectation = apply_problem_with_readout_conf(
         circuit, observable, noise_model, nshots, readout, qubit_map, backend=backend
     )
 
diff --git a/tests/test_models_error_mitigation.py b/tests/test_models_error_mitigation.py
index d8c0066793..22bccd5d48 100644
--- a/tests/test_models_error_mitigation.py
+++ b/tests/test_models_error_mitigation.py
@@ -8,7 +8,7 @@
     CDR,
     ICS,
     ZNE,
-    apply_readout_mitigation,
+    apply_problem_with_readout_conf,
     get_response_matrix,
     sample_clifford_training_circuit,
     sample_training_circuit_cdr,
@@ -19,11 +19,11 @@
 from qibo.symbols import Z
 
 
-def get_noise_model(error, gate, cal_matrix=[False, None]):
+def get_noise_model(error, gate, resp_matrix=[False, None]):
     noise = NoiseModel()
     noise.add(error, gate)
-    if cal_matrix[0]:
-        noise.add(ReadoutError(probabilities=cal_matrix[1]), gate=gates.M)
+    if resp_matrix[0]:
+        noise.add(ReadoutError(probabilities=resp_matrix[1]), gate=gates.M)
 
     return noise
 
@@ -53,11 +53,11 @@ def get_circuit(nqubits, nmeas=None):
 
 
 backend = construct_backend("numpy")
-# # Generate random calibration matrices
-cal_matrix_1q = random_stochastic_matrix(
+# # Generate random response matrices
+resp_matrix_1q = random_stochastic_matrix(
     2, diagonally_dominant=True, seed=2, backend=backend
 )
-cal_matrix_2q = random_stochastic_matrix(
+resp_matrix_2q = random_stochastic_matrix(
     4, diagonally_dominant=True, seed=2, backend=backend
 )
 
@@ -70,7 +70,7 @@ def get_circuit(nqubits, nmeas=None):
             3,
             get_noise_model(DepolarizingError(0.1), gates.CNOT),
             "CNOT",
-            {"calibration_matrix": cal_matrix_2q, "ibu_iters": None},
+            {"response_matrix": resp_matrix_2q, "ibu_iters": None},
         ),
         (
             3,
@@ -83,13 +83,13 @@ def get_circuit(nqubits, nmeas=None):
             1,
             get_noise_model(DepolarizingError(0.3), gates.RX),
             "RX",
-            {"calibration_matrix": cal_matrix_1q, "ibu_iters": None},
+            {"response_matrix": resp_matrix_1q, "ibu_iters": None},
         ),
         (
             1,
             get_noise_model(DepolarizingError(0.3), gates.RX),
             "RX",
-            {"calibration_matrix": cal_matrix_1q, "ibu_iters": 10},
+            {"response_matrix": resp_matrix_1q, "ibu_iters": 10},
         ),
         (1, get_noise_model(DepolarizingError(0.1), gates.RX), "RX", {"ncircuits": 2}),
     ],
@@ -143,15 +143,15 @@ def test_zne(backend, nqubits, noise, solve, insertion_gate, readout):
     [
         (get_noise_model(DepolarizingError(0.1), gates.CNOT), {}),
         (
-            get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, cal_matrix_2q]),
-            {"calibration_matrix": cal_matrix_2q, "ibu_iters": None},
+            get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, resp_matrix_2q]),
+            {"response_matrix": resp_matrix_2q, "ibu_iters": None},
         ),
         (
-            get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, cal_matrix_2q]),
-            {"calibration_matrix": cal_matrix_2q, "ibu_iters": 10},
+            get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, resp_matrix_2q]),
+            {"response_matrix": resp_matrix_2q, "ibu_iters": 10},
         ),
         (
-            get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, cal_matrix_2q]),
+            get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, resp_matrix_2q]),
             {"ncircuits": 2},
         ),
     ],
@@ -231,13 +231,13 @@ def test_sample_training_circuit(nqubits):
     [
         (
             1,
-            get_noise_model(DepolarizingError(0.1), gates.RX, [True, cal_matrix_1q]),
+            get_noise_model(DepolarizingError(0.1), gates.RX, [True, resp_matrix_1q]),
             "RX",
-            {"calibration_matrix": cal_matrix_1q, "ibu_iters": 10},
+            {"response_matrix": resp_matrix_1q, "ibu_iters": 10},
         ),
         (
             1,
-            get_noise_model(DepolarizingError(0.1), gates.RX, [True, cal_matrix_1q]),
+            get_noise_model(DepolarizingError(0.1), gates.RX, [True, resp_matrix_1q]),
             "RX",
             {"ncircuits": 2},
         ),
@@ -283,7 +283,7 @@ def test_vncdr(backend, nqubits, noise, full_output, insertion_gate, readout):
 
 
 @pytest.mark.parametrize("nqubits,nmeas", [(3, 2)])
-@pytest.mark.parametrize("method", ["cal_matrix", "randomized"])
+@pytest.mark.parametrize("method", ["response_matrix", "randomized"])
 @pytest.mark.parametrize("ibu_iters", [None, 10])
 def test_readout_mitigation(backend, nqubits, nmeas, method, ibu_iters):
     if backend.name == "tensorflow":
@@ -297,11 +297,11 @@ def test_readout_mitigation(backend, nqubits, nmeas, method, ibu_iters):
     p = random_stochastic_matrix(2**nmeas, diagonally_dominant=True, seed=5)
     noise = NoiseModel()
     noise.add(ReadoutError(probabilities=p), gate=gates.M)
-    if method == "cal_matrix":
-        calibration = get_response_matrix(
+    if method == "response_matrix":
+        response = get_response_matrix(
             nmeas, None, noise, nshots=nshots, backend=backend
         )
-        readout = {"calibration_matrix": calibration, "ibu_iters": ibu_iters}
+        readout = {"response_matrix": response, "ibu_iters": ibu_iters}
     elif method == "randomized":
         readout = {"ncircuits": 10}
     # Define the observable
@@ -315,7 +315,9 @@ def test_readout_mitigation(backend, nqubits, nmeas, method, ibu_iters):
     state = backend.execute_circuit(noise.apply(c), nshots=nshots)
     noisy_val = state.expectation_from_samples(obs)
 
-    mit_val = apply_readout_mitigation(c, obs, noise, nshots, readout, backend=backend)
+    mit_val = apply_problem_with_readout_conf(
+        c, obs, noise, nshots, readout, backend=backend
+    )
 
     assert np.abs(true_val - mit_val) <= np.abs(true_val - noisy_val)
 
@@ -327,15 +329,15 @@ def test_readout_mitigation(backend, nqubits, nmeas, method, ibu_iters):
     [
         (get_noise_model(DepolarizingError(0.1), gates.CNOT), {}),
         (
-            get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, cal_matrix_2q]),
-            {"calibration_matrix": cal_matrix_2q, "ibu_iters": None},
+            get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, resp_matrix_2q]),
+            {"response_matrix": resp_matrix_2q, "ibu_iters": None},
         ),
         (
-            get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, cal_matrix_2q]),
-            {"calibration_matrix": cal_matrix_2q, "ibu_iters": 10},
+            get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, resp_matrix_2q]),
+            {"response_matrix": resp_matrix_2q, "ibu_iters": 10},
         ),
         (
-            get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, cal_matrix_2q]),
+            get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, resp_matrix_2q]),
             {"ncircuits": 2},
         ),
     ],

From 94bbd4f1c35cb3792b9c801419699f5b2d06b613 Mon Sep 17 00:00:00 2001
From: Alejandro Sopena <alexsopena15@gmail.com>
Date: Wed, 10 Jan 2024 11:18:35 +0100
Subject: [PATCH 40/47] code review

---
 doc/source/api-reference/qibo.rst             |  6 +-
 doc/source/code-examples/advancedexamples.rst |  6 +-
 src/qibo/gates/special.py                     |  3 +-
 src/qibo/models/circuit.py                    |  2 +-
 src/qibo/models/error_mitigation.py           | 66 +++++++++----------
 tests/test_models_error_mitigation.py         |  4 +-
 6 files changed, 38 insertions(+), 49 deletions(-)

diff --git a/doc/source/api-reference/qibo.rst b/doc/source/api-reference/qibo.rst
index e63d2ec397..bbeeda603c 100644
--- a/doc/source/api-reference/qibo.rst
+++ b/doc/source/api-reference/qibo.rst
@@ -317,7 +317,7 @@ See `Nachman et al <https://arxiv.org/abs/1910.01969>`_ for more details.
 .. autofunction:: qibo.models.error_mitigation.apply_randomized_readout_mitigation
 
 
-.. autofunction:: qibo.models.error_mitigation.apply_problem_with_readout_conf
+.. autofunction:: qibo.models.error_mitigation.get_expectation_val_with_readout_mitigation
 
 
 Randomized readout mitigation
@@ -1465,10 +1465,6 @@ passing a symplectic matrix to the constructor.
    symplectic_matrix = backend.zero_state(nqubits=3)
    clifford = Clifford(symplectic_matrix, engine=NumpyBackend())
 
-   # The initialization above is equivalent to the initialization below
-   circuit = Circuit(nqubits=3)
-   clifford = Clifford(circuit, engine=NumpyBackend())
-
 The generators of the stabilizers can be extracted with the
 :meth:`qibo.quantum_info.clifford.Clifford.generators` method,
 or the complete set of :math:`d = 2^{n}` stabilizers operators can be extracted through the
diff --git a/doc/source/code-examples/advancedexamples.rst b/doc/source/code-examples/advancedexamples.rst
index 7325319de8..ca85dd2e25 100644
--- a/doc/source/code-examples/advancedexamples.rst
+++ b/doc/source/code-examples/advancedexamples.rst
@@ -1291,7 +1291,7 @@ response matrix and use it modify the final state after the circuit execution:
 
 .. testcode::
 
-   from qibo.models.error_mitigation import apply_problem_with_readout_conf, get_response_matrix
+   from qibo.models.error_mitigation import get_expectation_val_with_readout_mitigation, get_response_matrix
 
    nshots = 10000
    # compute the response matrix
@@ -1301,7 +1301,7 @@ response matrix and use it modify the final state after the circuit execution:
    # define mitigation options
    readout = {"response_matrix": response_matrix}
    # mitigate the readout errors
-   mit_val = apply_problem_with_readout_conf(circ, obs, noise, readout=readout)
+   mit_val = get_expectation_val_with_readout_mitigation(circ, obs, noise, readout=readout)
    print(mit_val)
    # 0.5945794816381054
 
@@ -1319,7 +1319,7 @@ Or use the randomized readout mitigation:
    # define mitigation options
    readout = {"ncircuits": 10}
    # mitigate the readout errors
-   mit_val = apply_problem_with_readout_conf(circ, obs, noise, readout=readout)
+   mit_val = get_expectation_val_with_readout_mitigation(circ, obs, noise, readout=readout)
    print(mit_val)
    # 0.5860884499785314
 
diff --git a/src/qibo/gates/special.py b/src/qibo/gates/special.py
index 4f1cb16cd5..d5027bb5ea 100644
--- a/src/qibo/gates/special.py
+++ b/src/qibo/gates/special.py
@@ -1,6 +1,5 @@
 from qibo.backends import GlobalBackend
 from qibo.gates.abstract import SpecialGate
-from qibo.gates.gates import I
 from qibo.gates.measurements import M
 
 
@@ -55,7 +54,7 @@ def __init__(self, *q):
     def from_gate(cls, gate):
         fgate = cls(*gate.qubits)
         fgate.append(gate)
-        if isinstance(gate, (M, SpecialGate, I)):
+        if isinstance(gate, (M, SpecialGate)):
             # special gates do not participate in fusion
             fgate.marked = True
         return fgate
diff --git a/src/qibo/models/circuit.py b/src/qibo/models/circuit.py
index e52398af06..419f5635d2 100644
--- a/src/qibo/models/circuit.py
+++ b/src/qibo/models/circuit.py
@@ -81,7 +81,7 @@ def from_fused(self):
                     queue.append(gate.gates[0])
                 else:
                     queue.append(gate)
-            elif isinstance(gate.gates[0], (gates.SpecialGate, gates.M, gates.I)):
+            elif isinstance(gate.gates[0], (gates.SpecialGate, gates.M)):
                 # special gates are marked by default so we need
                 # to add them manually
                 queue.append(gate.gates[0])
diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index 5f5f1c44ca..06c799262b 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -156,7 +156,7 @@ def ZNE(
         noisy_circuit = get_noisy_circuit(
             circuit, num_insertions, insertion_gate=insertion_gate
         )
-        val = apply_problem_with_readout_conf(
+        val = get_expectation_val_with_readout_mitigation(
             noisy_circuit,
             observable,
             noise_model,
@@ -312,14 +312,14 @@ def CDR(
     for circ in training_circuits:
         val = circ(nshots=nshots).expectation_from_samples(observable)
         train_val["noise-free"].append(val)
-        val = apply_problem_with_readout_conf(
+        val = get_expectation_val_with_readout_mitigation(
             circ, observable, noise_model, nshots, readout, qubit_map, backend=backend
         )
         train_val["noisy"].append(val)
 
     optimal_params = curve_fit(model, train_val["noisy"], train_val["noise-free"])[0]
 
-    val = apply_problem_with_readout_conf(
+    val = get_expectation_val_with_readout_mitigation(
         circuit, observable, noise_model, nshots, readout, qubit_map, backend=backend
     )
     mit_val = model(val, *optimal_params)
@@ -398,7 +398,7 @@ def vnCDR(
         train_val["noise-free"].append(val)
         for level in noise_levels:
             noisy_c = get_noisy_circuit(circ, level, insertion_gate=insertion_gate)
-            val = apply_problem_with_readout_conf(
+            val = get_expectation_val_with_readout_mitigation(
                 noisy_c,
                 observable,
                 noise_model,
@@ -417,7 +417,7 @@ def vnCDR(
     val = []
     for level in noise_levels:
         noisy_c = get_noisy_circuit(circuit, level, insertion_gate=insertion_gate)
-        expval = apply_problem_with_readout_conf(
+        expval = get_expectation_val_with_readout_mitigation(
             noisy_c,
             observable,
             noise_model,
@@ -501,7 +501,7 @@ def get_response_matrix(
                 circuit.add(gates.X(qubit))
         circuit.add(gates.M(*range(nqubits)))
 
-        circuit_result = _circuit_conf(
+        circuit_result = _execute_circuit(
             circuit, qubit_map, noise_model, nshots, backend=backend
         )
 
@@ -537,8 +537,7 @@ def apply_resp_mat_readout_mitigation(state, response_matrix, iterations=None):
 
     if iterations is None:
         calibration_matrix = np.linalg.inv(response_matrix)
-        for i, value in enumerate(calibration_matrix @ frequencies):
-            state._frequencies[i] = float(value)
+        mitigated_frequencies = calibration_matrix @ frequencies
     else:
         mitigated_probabilities = iterative_bayesian_unfolding(
             frequencies / np.sum(frequencies), response_matrix, iterations
@@ -550,8 +549,8 @@ def apply_resp_mat_readout_mitigation(state, response_matrix, iterations=None):
             mitigated_frequencies / np.sum(mitigated_frequencies)
         ) * np.sum(frequencies)
 
-        for i, value in enumerate(mitigated_frequencies):
-            state._frequencies[i] = float(value)
+    for i, value in enumerate(mitigated_frequencies):
+        state._frequencies[i] = float(value)
 
     return state
 
@@ -609,10 +608,6 @@ def apply_randomized_readout_mitigation(
         for j, gate in enumerate(x_gate):
             if isinstance(gate, gates.X) and gate.qubits[0] in meas_qubits:
                 error_map[gate.qubits[0]] = 1
-                # if gate.qubits[0] in meas_qubits:
-                #     error_map[gate.qubits[0]] = 1
-                # else:
-                #     x_gate.queue[j] = gates.I(gate.qubits[0])
 
         circuits = [circuit_c, cal_circuit]
         results = []
@@ -621,7 +616,7 @@ def apply_randomized_readout_mitigation(
             circ.add(x_gate)
             circ.add(gates.M(*meas_qubits))
 
-            result = _circuit_conf(
+            result = _execute_circuit(
                 circ, qubit_map, noise_model, nshots_r, backend=backend
             )
             result._samples = result.apply_bitflips(error_map)
@@ -639,7 +634,7 @@ def apply_randomized_readout_mitigation(
     return results
 
 
-def apply_problem_with_readout_conf(
+def get_expectation_val_with_readout_mitigation(
     circuit,
     observable,
     noise_model=None,
@@ -677,16 +672,15 @@ def apply_problem_with_readout_conf(
             circuit, noise_model, nshots, readout["ncircuits"], backend
         )
     else:
-        circuit_result = _circuit_conf(
+        circuit_result = _execute_circuit(
             circuit, qubit_map, noise_model, nshots, backend=backend
         )
-
-    if "response_matrix" in readout:
-        circuit_result = apply_resp_mat_readout_mitigation(
-            circuit_result,
-            readout["response_matrix"],
-            readout.get("ibu_iters", None),
-        )
+        if "response_matrix" in readout:
+            circuit_result = apply_resp_mat_readout_mitigation(
+                circuit_result,
+                readout["response_matrix"],
+                readout.get("ibu_iters", None),
+            )
 
     exp_val = circuit_result.expectation_from_samples(observable)
 
@@ -718,13 +712,13 @@ def sample_clifford_training_circuit(
     if backend is None:  # pragma: no cover
         backend = GlobalBackend()
 
-    non_clifford_gates = [
-        (i, gate)
+    non_clifford_gates_indices = [
+        i
         for i, gate in enumerate(circuit.queue)
         if not gate.clifford and not isinstance(gate, gates.M)
     ]
 
-    if not non_clifford_gates:
+    if not non_clifford_gates_indices:
         raise_error(ValueError, "No non-Clifford gate found, no circuit sampled.")
 
     sampled_circuit = circuit.__class__(**circuit.init_kwargs)
@@ -740,7 +734,7 @@ def sample_clifford_training_circuit(
                 sampled_circuit.add(gate_rand)
             sampled_circuit.add(gate)
         else:
-            if i in [index for index, _ in non_clifford_gates]:
+            if i in non_clifford_gates_indices:
                 gate = gates.Unitary(
                     random_clifford(1, backend=backend, return_circuit=False),
                     gate.qubits[0],
@@ -894,7 +888,7 @@ def ICS(
         circuit_result = training_circuit(nshots=nshots)
         expectation = observable.expectation_from_samples(circuit_result.frequencies())
 
-        noisy_expectation = apply_problem_with_readout_conf(
+        noisy_expectation = get_expectation_val_with_readout_mitigation(
             training_circuit,
             observable,
             noise_model,
@@ -911,20 +905,20 @@ def ICS(
     dep_param = np.mean(lambda_list)
     dep_param_std = np.std(lambda_list)
 
-    noisy_expectation = apply_problem_with_readout_conf(
+    noisy_expectation = get_expectation_val_with_readout_mitigation(
         circuit, observable, noise_model, nshots, readout, qubit_map, backend=backend
     )
+    one_dep_squared = (1 - dep_param) ** 2
+    dep_std_squared = dep_param_std**2
 
     mitigated_expectation = (
-        (1 - dep_param)
-        * noisy_expectation
-        / ((1 - dep_param) ** 2 + dep_param_std**2)
+        (1 - dep_param) * noisy_expectation / (one_dep_squared + dep_std_squared)
     )
     mitigated_expectation_std = (
         dep_param_std
         * abs(noisy_expectation)
-        * abs((1 - dep_param) ** 2 - dep_param_std**2)
-        / ((1 - dep_param) ** 2 + dep_param_std**2) ** 2
+        * abs((1 - dep_param) ** 2 - dep_std_squared)
+        / (one_dep_squared + dep_std_squared) ** 2
     )
 
     if full_output:
@@ -940,7 +934,7 @@ def ICS(
     return mitigated_expectation
 
 
-def _circuit_conf(circuit, qubit_map, noise_model=None, nshots=10000, backend=None):
+def _execute_circuit(circuit, qubit_map, noise_model=None, nshots=10000, backend=None):
     """
     Helper function to execute the given circuit with the specified parameters.
 
diff --git a/tests/test_models_error_mitigation.py b/tests/test_models_error_mitigation.py
index 22bccd5d48..c1132c5a5a 100644
--- a/tests/test_models_error_mitigation.py
+++ b/tests/test_models_error_mitigation.py
@@ -8,7 +8,7 @@
     CDR,
     ICS,
     ZNE,
-    apply_problem_with_readout_conf,
+    get_expectation_val_with_readout_mitigation,
     get_response_matrix,
     sample_clifford_training_circuit,
     sample_training_circuit_cdr,
@@ -315,7 +315,7 @@ def test_readout_mitigation(backend, nqubits, nmeas, method, ibu_iters):
     state = backend.execute_circuit(noise.apply(c), nshots=nshots)
     noisy_val = state.expectation_from_samples(obs)
 
-    mit_val = apply_problem_with_readout_conf(
+    mit_val = get_expectation_val_with_readout_mitigation(
         c, obs, noise, nshots, readout, backend=backend
     )
 

From d851a8efcbb6786f5c0410a94d755e857ef6add0 Mon Sep 17 00:00:00 2001
From: Renato Mello <renato.msf@gmail.com>
Date: Thu, 11 Jan 2024 23:29:21 +0000
Subject: [PATCH 41/47] Update src/qibo/models/error_mitigation.py

Co-authored-by: Matteo Robbiati <62071516+MatteoRobbiati@users.noreply.github.com>
---
 src/qibo/models/error_mitigation.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index 06c799262b..7139377432 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -662,7 +662,7 @@ def get_expectation_val_with_readout_mitigation(
             If None, it uses the global backend. Defaults to ``None``.
 
     Returns:
-        float: The mitigated expectation value of the observable.
+        float: the mitigated expectation value of the observable.
     """
     if backend is None:  # pragma: no cover
         backend = GlobalBackend()

From 7746326e595c61283ed3e3c651648f85dfb1cf7d Mon Sep 17 00:00:00 2001
From: Renato Mello <renato.msf@gmail.com>
Date: Thu, 11 Jan 2024 23:29:33 +0000
Subject: [PATCH 42/47] Update src/qibo/models/error_mitigation.py

Co-authored-by: Matteo Robbiati <62071516+MatteoRobbiati@users.noreply.github.com>
---
 src/qibo/models/error_mitigation.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index 7139377432..0181a0e7dc 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -703,7 +703,7 @@ def sample_clifford_training_circuit(
             Defaults to ``None``.
 
     Returns:
-        :class:`qibo.models.Circuit`: The sampled circuit.
+        :class:`qibo.models.Circuit`: the sampled circuit.
     """
     from qibo.quantum_info import (  # pylint: disable=import-outside-toplevel
         random_clifford,

From 1002a694ea25d2151f2b5988782671c878a936f2 Mon Sep 17 00:00:00 2001
From: Renato Mello <renato.msf@gmail.com>
Date: Thu, 11 Jan 2024 23:29:52 +0000
Subject: [PATCH 43/47] Update src/qibo/models/error_mitigation.py

Co-authored-by: Matteo Robbiati <62071516+MatteoRobbiati@users.noreply.github.com>
---
 src/qibo/models/error_mitigation.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index 0181a0e7dc..8f9eefd0c3 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -527,7 +527,7 @@ def apply_resp_mat_readout_mitigation(state, response_matrix, iterations=None):
             If ``None`` the 'inverse' method is used. Defaults to ``None``.
 
     Returns:
-        :class:`qibo.measurements.CircuitResult`: The input state with the updated (mitigated) frequencies.
+        :class:`qibo.measurements.CircuitResult`: the input state with the updated (mitigated) frequencies.
     """
     frequencies = np.zeros(2 ** len(state.measurements[0].qubits))
     for key, value in state.frequencies().items():

From f05141e83c07db2424510e4cfb18b980975c3448 Mon Sep 17 00:00:00 2001
From: Renato Mello <renato.msf@gmail.com>
Date: Thu, 11 Jan 2024 23:30:13 +0000
Subject: [PATCH 44/47] Update src/qibo/models/error_mitigation.py

Co-authored-by: Matteo Robbiati <62071516+MatteoRobbiati@users.noreply.github.com>
---
 src/qibo/models/error_mitigation.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index 8f9eefd0c3..886351e7e1 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -482,7 +482,7 @@ def get_response_matrix(
             Defaults to ``None``.
 
     Returns:
-        numpy.ndarray : The computed (`nqubits`, `nqubits`) response matrix for
+        numpy.ndarray : the computed (`nqubits`, `nqubits`) response matrix for
             readout mitigation.
     """
     from qibo import Circuit  # pylint: disable=import-outside-toplevel

From 159707bcee4c8e42c9932d015d6531386ab8d68a Mon Sep 17 00:00:00 2001
From: Alejandro Sopena <alexsopena15@gmail.com>
Date: Fri, 12 Jan 2024 01:12:35 +0100
Subject: [PATCH 45/47] Matteo's suggestions

---
 src/qibo/models/error_mitigation.py | 20 ++++++++++++++------
 1 file changed, 14 insertions(+), 6 deletions(-)

diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index 886351e7e1..036bb53c1a 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -108,7 +108,7 @@ def ZNE(
     nshots=10000,
     solve_for_gammas=False,
     insertion_gate="CNOT",
-    readout: dict = {},
+    readout=None,
     qubit_map=None,
     backend=None,
 ):
@@ -150,6 +150,8 @@ def ZNE(
     """
     if backend is None:  # pragma: no cover
         backend = GlobalBackend()
+    if readout is None:
+        readout = {}
 
     expected_values = []
     for num_insertions in noise_levels:
@@ -261,7 +263,7 @@ def CDR(
     model=lambda x, a, b: a * x + b,
     n_training_samples: int = 100,
     full_output: bool = False,
-    readout: dict = {},
+    readout=None,
     qubit_map=None,
     backend=None,
 ):
@@ -303,6 +305,8 @@ def CDR(
     """
     if backend is None:  # pragma: no cover
         backend = GlobalBackend()
+    if readout is None:
+        readout = {}
 
     training_circuits = [
         sample_training_circuit_cdr(circuit) for _ in range(n_training_samples)
@@ -340,7 +344,7 @@ def vnCDR(
     n_training_samples: int = 100,
     insertion_gate: str = "CNOT",
     full_output: bool = False,
-    readout: dict = {},
+    readout=None,
     qubit_map=None,
     backend=None,
 ):
@@ -387,6 +391,8 @@ def vnCDR(
     """
     if backend is None:  # pragma: no cover
         backend = GlobalBackend()
+    if readout is None:
+        readout = {}
 
     training_circuits = [
         sample_training_circuit_cdr(circuit) for _ in range(n_training_samples)
@@ -694,7 +700,7 @@ def sample_clifford_training_circuit(
     circuit,
     backend=None,
 ):
-    """Samples a training circuit for CDR by susbtituting some of the non-Clifford gates.
+    """Samples a training circuit for CDR by susbtituting all the non-Clifford gates.
 
     Args:
         circuit (:class:`qibo.models.Circuit`): circuit to sample from.
@@ -736,8 +742,10 @@ def sample_clifford_training_circuit(
         else:
             if i in non_clifford_gates_indices:
                 gate = gates.Unitary(
-                    random_clifford(1, backend=backend, return_circuit=False),
-                    gate.qubits[0],
+                    random_clifford(
+                        len(gate.qubits), backend=backend, return_circuit=False
+                    ),
+                    *gate.qubits,
                 )
                 gate.clifford = True
             sampled_circuit.add(gate)

From bee350a8053bab97caf24fb1a67be706ffefb37b Mon Sep 17 00:00:00 2001
From: Alejandro Sopena <alexsopena15@gmail.com>
Date: Fri, 12 Jan 2024 02:01:58 +0100
Subject: [PATCH 46/47] fix coverage

---
 src/qibo/models/error_mitigation.py   | 8 ++++++--
 tests/test_models_error_mitigation.py | 8 ++++----
 2 files changed, 10 insertions(+), 6 deletions(-)

diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index 036bb53c1a..a1f2c4ec50 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -645,7 +645,7 @@ def get_expectation_val_with_readout_mitigation(
     observable,
     noise_model=None,
     nshots: int = 10000,
-    readout={},
+    readout=None,
     qubit_map=None,
     backend=None,
 ):
@@ -672,6 +672,8 @@ def get_expectation_val_with_readout_mitigation(
     """
     if backend is None:  # pragma: no cover
         backend = GlobalBackend()
+    if readout is None:
+        readout = {}
 
     if "ncircuits" in readout:
         circuit_result, circuit_result_cal = apply_randomized_readout_mitigation(
@@ -837,7 +839,7 @@ def error_sensitive_circuit(circuit, observable, backend=None):
 def ICS(
     circuit,
     observable,
-    readout={},
+    readout=None,
     qubit_map=None,
     noise_model=None,
     nshots=int(1e4),
@@ -880,6 +882,8 @@ def ICS(
     """
     if backend is None:  # pragma: no cover
         backend = GlobalBackend()
+    if readout is None:
+        readout = {}
 
     if qubit_map is None:
         qubit_map = list(range(circuit.nqubits))
diff --git a/tests/test_models_error_mitigation.py b/tests/test_models_error_mitigation.py
index c1132c5a5a..7d64b2481b 100644
--- a/tests/test_models_error_mitigation.py
+++ b/tests/test_models_error_mitigation.py
@@ -65,7 +65,7 @@ def get_circuit(nqubits, nmeas=None):
 @pytest.mark.parametrize(
     "nqubits,noise,insertion_gate,readout",
     [
-        (3, get_noise_model(DepolarizingError(0.1), gates.CNOT), "CNOT", {}),
+        (3, get_noise_model(DepolarizingError(0.1), gates.CNOT), "CNOT", None),
         (
             3,
             get_noise_model(DepolarizingError(0.1), gates.CNOT),
@@ -78,7 +78,7 @@ def get_circuit(nqubits, nmeas=None):
             "CNOT",
             {"ncircuits": 2},
         ),
-        (1, get_noise_model(DepolarizingError(0.1), gates.RX), "RX", {}),
+        (1, get_noise_model(DepolarizingError(0.1), gates.RX), "RX", None),
         (
             1,
             get_noise_model(DepolarizingError(0.3), gates.RX),
@@ -141,7 +141,7 @@ def test_zne(backend, nqubits, noise, solve, insertion_gate, readout):
 @pytest.mark.parametrize(
     "noise,readout",
     [
-        (get_noise_model(DepolarizingError(0.1), gates.CNOT), {}),
+        (get_noise_model(DepolarizingError(0.1), gates.CNOT), None),
         (
             get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, resp_matrix_2q]),
             {"response_matrix": resp_matrix_2q, "ibu_iters": None},
@@ -327,7 +327,7 @@ def test_readout_mitigation(backend, nqubits, nmeas, method, ibu_iters):
 @pytest.mark.parametrize(
     "noise,readout",
     [
-        (get_noise_model(DepolarizingError(0.1), gates.CNOT), {}),
+        (get_noise_model(DepolarizingError(0.1), gates.CNOT), None),
         (
             get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, resp_matrix_2q]),
             {"response_matrix": resp_matrix_2q, "ibu_iters": None},

From 9cebb486ee1aec92e1bf432d671d8badd3cb30ea Mon Sep 17 00:00:00 2001
From: Alejandro Sopena <alexsopena15@gmail.com>
Date: Fri, 12 Jan 2024 04:01:03 +0100
Subject: [PATCH 47/47] fix coverage

---
 src/qibo/models/error_mitigation.py   | 2 +-
 tests/test_models_error_mitigation.py | 1 +
 2 files changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index a1f2c4ec50..2e2d577d2e 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -672,7 +672,7 @@ def get_expectation_val_with_readout_mitigation(
     """
     if backend is None:  # pragma: no cover
         backend = GlobalBackend()
-    if readout is None:
+    if readout is None:  # pragma: no cover
         readout = {}
 
     if "ncircuits" in readout:
diff --git a/tests/test_models_error_mitigation.py b/tests/test_models_error_mitigation.py
index 7d64b2481b..035daabd1b 100644
--- a/tests/test_models_error_mitigation.py
+++ b/tests/test_models_error_mitigation.py
@@ -229,6 +229,7 @@ def test_sample_training_circuit(nqubits):
 @pytest.mark.parametrize(
     "nqubits,noise,insertion_gate,readout",
     [
+        (1, get_noise_model(DepolarizingError(0.1), gates.RX), "RX", None),
         (
             1,
             get_noise_model(DepolarizingError(0.1), gates.RX, [True, resp_matrix_1q]),