[Feature] Correlated readout

docs/tutorials/realistic_sims/mitigation.md

@@ -45,7 +45,7 @@ observable = hamiltonian_factory(circuit.n_qubits, detuning=Z)
 model = QuantumModel(circuit=circuit, observable=observable)
 
 # Define a noise model to use:
-noise = NoiseHandler(NoiseProtocol.READOUT)
+noise = NoiseHandler(NoiseProtocol.READOUT.INDEPENDENT)
 # Define the mitigation method solving the minimization problem:
 options={"optimization_type": ReadOutOptimization.CONSTRAINED}  # ReadOutOptimization.MLE for the alternative method.
 mitigation = Mitigations(protocol=Mitigations.READOUT, options=options)

docs/tutorials/realistic_sims/noise.md

@@ -13,7 +13,7 @@ from qadence.types import NoiseProtocol
 
 analog_noise = NoiseHandler(protocol=NoiseProtocol.ANALOG.DEPOLARIZING, options={"noise_probs": 0.1})
 digital_noise = NoiseHandler(protocol=NoiseProtocol.DIGITAL.DEPOLARIZING, options={"error_probability": 0.1})
-readout_noise = NoiseHandler(protocol=NoiseProtocol.READOUT, options={"error_probability": 0.1, "seed": 0})
+readout_noise = NoiseHandler(protocol=NoiseProtocol.READOUT.INDEPENDENT, options={"error_probability": 0.1, "seed": 0})
 ```
 
 One can also define a `NoiseHandler` passing a list of protocols and a list of options (careful with the order):
@@ -36,7 +36,7 @@ from qadence import NoiseHandler
 from qadence.types import NoiseProtocol
 
 depo_noise = NoiseHandler(protocol=NoiseProtocol.DIGITAL.DEPOLARIZING, options={"error_probability": 0.1})
-readout_noise = NoiseHandler(protocol=NoiseProtocol.READOUT, options={"error_probability": 0.1, "seed": 0})
+readout_noise = NoiseHandler(protocol=NoiseProtocol.READOUT.INDEPENDENT, options={"error_probability": 0.1, "seed": 0})
 
 noise_combination = NoiseHandler(protocol=NoiseProtocol.DIGITAL.BITFLIP, options={"error_probability": 0.1})
 noise_combination.append([depo_noise, readout_noise])
@@ -49,7 +49,7 @@ Finally, one can add directly a few pre-defined types using several `NoiseHandle
 from qadence import NoiseHandler
 from qadence.types import NoiseProtocol
 noise_combination = NoiseHandler(protocol=NoiseProtocol.DIGITAL.BITFLIP, options={"error_probability": 0.1})
-noise_combination.digital_depolarizing({"error_probability": 0.1}).readout({"error_probability": 0.1, "seed": 0})
+noise_combination.digital_depolarizing({"error_probability": 0.1}).readout_independent({"error_probability": 0.1, "seed": 0})
 print(noise_combination)
 ```
 
@@ -65,6 +65,10 @@ $$
 T(x|x')=\delta_{xx'}
 $$
 
+Two types of readout protocols are available:
+- `NoiseProtocol.READOUT.INDEPENDENT` where each bit can be corrupted independently of each other.
+- `NoiseProtocol.READOUT.CORRELATED` where we can define of confusion matrix of corruption between each
+possible bitstrings.
 
 Qadence offers to simulate readout errors with the `NoiseHandler` to corrupt the output
 samples of a simulation, through execution via a `QuantumModel`:
@@ -82,7 +86,7 @@ observable = hamiltonian_factory(circuit.n_qubits, detuning=Z)
 model = QuantumModel(circuit=circuit, observable=observable)
 
 # Define a noise model to use.
-noise = NoiseHandler(protocol=NoiseProtocol.READOUT)
+noise = NoiseHandler(protocol=NoiseProtocol.READOUT.INDEPENDENT)
 
 # Run noiseless and noisy simulations.
 noiseless_samples = model.sample(n_shots=100)
@@ -93,12 +97,21 @@ print(f"noisy = {noisy_samples}") # markdown-exec: hide
 ```
 
 It is possible to pass options to the noise model. In the previous example, a noise matrix is implicitly computed from a
-uniform distribution. The `option` dictionary argument accepts the following options:
+uniform distribution.
+
+For `NoiseProtocol.READOUT.INDEPENDENT`, the `option` dictionary argument accepts the following options:
 
 - `seed`: defaulted to `None`, for reproducibility purposes
-- `error_probability`: defaulted to 0.1, a bit flip probability
+- `error_probability`: If float, the same probability is applied to every bit. By default, this is 0.1.
+    If a 1D tensor with the number of elements equal to the number of qubits, a different probability can be set for each qubit. If a tensor of shape (n_qubits, 2, 2) is passed, that is a confusion matrix obtained from experiments, we extract the error_probability.
+    and do not compute internally the confusion matrix as in the other cases.
 - `noise_distribution`: defaulted to `WhiteNoise.UNIFORM`, for non-uniform noise distributions
 
+For `NoiseProtocol.READOUT.CORRELATED`, the `option` dictionary argument accepts the following options:
+- `confusion_matrix`: The square matrix representing $T(x|x')$ for each possible bitstring of length `n` qubits. Should be of size (2**n, 2**n).
+- `seed`: defaulted to `None`, for reproducibility purposes
+
+
 Noisy simulations go hand-in-hand with measurement protocols discussed in the previous [section](measurements.md), to assess the impact of noise on expectation values. In this case, both measurement and noise protocols have to be defined appropriately. Please note that a noise protocol without a measurement protocol will be ignored for expectation values computations.
 
 
@@ -107,7 +120,7 @@ from qadence.measurements import Measurements
 
 # Define a noise model with options.
 options = {"error_probability": 0.01}
-noise = NoiseHandler(protocol=NoiseProtocol.READOUT, options=options)
+noise = NoiseHandler(protocol=NoiseProtocol.READOUT.INDEPENDENT, options=options)
 
 # Define a tomographical measurement protocol with options.
 options = {"n_shots": 10000}

qadence/backends/pyqtorch/convert_ops.py

@@ -373,4 +373,8 @@ def convert_readout_noise(n_qubits: int, noise: NoiseHandler) -> pyq.noise.Reado
     readout_part = noise.filter(NoiseProtocol.READOUT)
     if readout_part is None:
         return None
-    return pyq.noise.ReadoutNoise(n_qubits, **readout_part.options[0])
+
+    if readout_part.protocol[0] == NoiseProtocol.READOUT.INDEPENDENT:
+        return pyq.noise.ReadoutNoise(n_qubits, **readout_part.options[0])
+    else:
+        return pyq.noise.CorrelatedReadoutNoise(**readout_part.options[0])

qadence/noise/protocols.py

@@ -53,7 +53,7 @@ def __init__(
         self.verify_all_protocols()
 
     def _verify_single_protocol(self, protocol: NoiseEnum, option: dict) -> None:
-        if protocol != NoiseProtocol.READOUT:
+        if not isinstance(protocol, NoiseProtocol.READOUT):  # type: ignore[arg-type]
             name_mandatory_option = (
                 "noise_probs" if isinstance(protocol, NoiseProtocol.ANALOG) else "error_probability"
             )
@@ -86,10 +86,10 @@ def verify_all_protocols(self) -> None:
             if types.count(NoiseProtocol.ANALOG) > 1:
                 raise ValueError("Multiple Analog Noises are not supported yet.")
 
-        if NoiseProtocol.READOUT in self.protocol:
+        if NoiseProtocol.READOUT in unique_types:
             if (
-                self.protocol[-1] != NoiseProtocol.READOUT
-                or self.protocol.count(NoiseProtocol.READOUT) > 1
+                not isinstance(self.protocol[-1], NoiseProtocol.READOUT)
+                or types.count(NoiseProtocol.READOUT) > 1
             ):
                 raise ValueError("Only define a NoiseHandler with one READOUT as the last Noise.")
 
@@ -149,11 +149,7 @@ def list(cls) -> list:
         return list(filter(lambda el: not el.startswith("__"), dir(cls)))
 
     def filter(self, protocol: NoiseEnum) -> NoiseHandler | None:
-        protocol_matches: list = list()
-        if protocol == NoiseProtocol.READOUT:
-            protocol_matches = [p == protocol for p in self.protocol]
-        else:
-            protocol_matches = [isinstance(p, protocol) for p in self.protocol]  # type: ignore[arg-type]
+        protocol_matches: list = [isinstance(p, protocol) for p in self.protocol]  # type: ignore[arg-type]
 
         # if we have at least a match
         if True in protocol_matches:
@@ -201,6 +197,10 @@ def dephasing(self, *args: Any, **kwargs: Any) -> NoiseHandler:
         self.append(NoiseHandler(NoiseProtocol.ANALOG.DEPHASING, *args, **kwargs))
         return self
 
-    def readout(self, *args: Any, **kwargs: Any) -> NoiseHandler:
-        self.append(NoiseHandler(NoiseProtocol.READOUT, *args, **kwargs))
+    def readout_independent(self, *args: Any, **kwargs: Any) -> NoiseHandler:
+        self.append(NoiseHandler(NoiseProtocol.READOUT.INDEPENDENT, *args, **kwargs))
+        return self
+
+    def readout_correlated(self, *args: Any, **kwargs: Any) -> NoiseHandler:
+        self.append(NoiseHandler(NoiseProtocol.READOUT.CORRELATED, *args, **kwargs))
         return self

qadence/types.py

@@ -470,16 +470,25 @@ class AnalogNoise(StrEnum):
     DEPHASING = "Dephasing"
 
 
+class ReadoutNoise(StrEnum):
+    """Type of readout protocol."""
+
+    INDEPENDENT = "Independent Readout"
+    """Simple readout protocols where each qubit is corrupted independently."""
+    CORRELATED = "Correlated Readout"
+    """Using a confusion matrix (2**n, 2**n) for corrupting bitstrings values."""
+
+
 @dataclass
 class NoiseProtocol:
     """Type of noise protocol."""
 
     ANALOG = AnalogNoise
     """Noise applied in analog blocks."""
-    READOUT = "Readout"
+    READOUT = ReadoutNoise
     """Noise applied on outputs of quantum programs."""
     DIGITAL = DigitalNoise
     """Noise applied to digital blocks."""
 
 
-NoiseEnum = Union[DigitalNoise, AnalogNoise, str]
+NoiseEnum = Union[DigitalNoise, AnalogNoise, ReadoutNoise]

tests/qadence/test_noise/test_readout.py

@@ -32,40 +32,57 @@
 
 @pytest.mark.flaky(max_runs=5)
 @pytest.mark.parametrize(
-    "error_probability, n_shots, block, backend",
+    "error_probability, n_shots, block",
     [
-        (0.1, 100, kron(X(0), X(1)), BackendName.PYQTORCH),
-        (0.1, 200, kron(Z(0), Z(1), Z(2)) + kron(X(0), Y(1), Z(2)), BackendName.PYQTORCH),
-        (0.01, 1000, add(Z(0), Z(1), Z(2)), BackendName.PYQTORCH),
+        (
+            0.1,
+            100,
+            kron(X(0), X(1)),
+        ),
+        (
+            0.1,
+            200,
+            kron(Z(0), Z(1), Z(2)) + kron(X(0), Y(1), Z(2)),
+        ),
+        (
+            0.01,
+            1000,
+            add(Z(0), Z(1), Z(2)),
+        ),
         (
             0.1,
             2000,
             HamEvo(
                 generator=kron(X(0), X(1)) + kron(Z(0), Z(1)) + kron(X(2), X(3)), parameter=0.005
             ),
-            BackendName.PYQTORCH,
         ),
-        (0.1, 500, add(Z(0), Z(1), kron(X(2), X(3))) + add(X(2), X(3)), BackendName.PYQTORCH),
-        (0.05, 10000, add(kron(Z(0), Z(1)), kron(X(2), X(3))), BackendName.PYQTORCH),
-        (0.2, 1000, hamiltonian_factory(4, detuning=Z), BackendName.PYQTORCH),
-        (0.1, 500, kron(Z(0), Z(1)) + CNOT(0, 1), BackendName.PYQTORCH),
+        (
+            0.1,
+            500,
+            add(Z(0), Z(1), kron(X(2), X(3))) + add(X(2), X(3)),
+        ),
+        (0.05, 10000, add(kron(Z(0), Z(1)), kron(X(2), X(3)))),
+        (0.2, 1000, hamiltonian_factory(4, detuning=Z)),
+        (0.1, 500, kron(Z(0), Z(1)) + CNOT(0, 1)),
     ],
 )
 def test_readout_error_quantum_model(
     error_probability: float,
     n_shots: int,
     block: AbstractBlock,
-    backend: BackendName,
 ) -> None:
-    diff_mode = "ad" if backend == BackendName.PYQTORCH else "gpsr"
-
-    noiseless_samples: list[Counter] = QuantumModel(
+    backend = BackendName.PYQTORCH
+    diff_mode = "ad"
+    model = QuantumModel(
         QuantumCircuit(block.n_qubits, block), backend=backend, diff_mode=diff_mode
-    ).sample(n_shots=n_shots)
+    )
+    noiseless_samples: list[Counter] = model.sample(n_shots=n_shots)
 
-    noisy_samples: list[Counter] = QuantumModel(
-        QuantumCircuit(block.n_qubits, block), backend=backend, diff_mode=diff_mode
-    ).sample(noise=NoiseHandler(protocol=NoiseProtocol.READOUT), n_shots=n_shots)
+    noise_protocol: NoiseHandler = NoiseHandler(
+        protocol=NoiseProtocol.READOUT.INDEPENDENT,
+        options={"error_probability": error_probability},
+    )
+    noisy_samples: list[Counter] = model.sample(noise=noise_protocol, n_shots=n_shots)
 
     for noiseless, noisy in zip(noiseless_samples, noisy_samples):
         assert sum(noiseless.values()) == sum(noisy.values()) == n_shots
@@ -76,6 +93,23 @@ def test_readout_error_quantum_model(
             atol=1e-1,
         )
 
+    rand_confusion = torch.rand(2**block.n_qubits, 2**block.n_qubits)
+    rand_confusion = rand_confusion / rand_confusion.sum(dim=1, keepdim=True)
+    corr_noise_protocol: NoiseHandler = NoiseHandler(
+        protocol=NoiseProtocol.READOUT.CORRELATED,
+        options={"confusion_matrix": rand_confusion},
+    )
+    # assert difference with noiseless samples
+    corr_noisy_samples: list[Counter] = model.sample(noise=corr_noise_protocol, n_shots=n_shots)
+    for noiseless, noisy in zip(noiseless_samples, corr_noisy_samples):
+        assert sum(noiseless.values()) == sum(noisy.values()) == n_shots
+        assert js_divergence(noiseless, noisy) > 0.0
+
+    # assert difference noisy samples
+    for noisy, corr_noisy in zip(noisy_samples, corr_noisy_samples):
+        assert sum(noisy.values()) == sum(corr_noisy.values()) == n_shots
+        assert js_divergence(noisy, corr_noisy) > 0.0
+
 
 @pytest.mark.parametrize("backend", [BackendName.PYQTORCH, BackendName.PULSER])
 def test_readout_error_backends(backend: BackendName) -> None:
@@ -88,7 +122,7 @@ def test_readout_error_backends(backend: BackendName) -> None:
     samples = qd.sample(feature_map, n_shots=1000, values=inputs, backend=backend, noise=None)
     # introduce noise
     options = {"error_probability": error_probability}
-    noise = NoiseHandler(protocol=NoiseProtocol.READOUT, options=options)
+    noise = NoiseHandler(protocol=NoiseProtocol.READOUT.INDEPENDENT, options=options)
     noisy_samples = qd.sample(
         feature_map, n_shots=1000, values=inputs, backend=backend, noise=noise
     )
@@ -120,7 +154,7 @@ def test_readout_error_with_measurements(
     observable = hamiltonian_factory(circuit.n_qubits, detuning=Z)
 
     model = QuantumModel(circuit=circuit, observable=observable, diff_mode=DiffMode.GPSR)
-    noise = NoiseHandler(protocol=NoiseProtocol.READOUT)
+    noise = NoiseHandler(protocol=NoiseProtocol.READOUT.INDEPENDENT)
     measurement = Measurements(protocol=str(measurement_proto), options=options)
 
     noisy = model.expectation(values=inputs, measurement=measurement, noise=noise)
@@ -137,7 +171,16 @@ def test_readout_error_with_measurements(
 
 
 def test_serialization() -> None:
-    noise = NoiseHandler(protocol=NoiseProtocol.READOUT)
+    noise = NoiseHandler(protocol=NoiseProtocol.READOUT.INDEPENDENT)
+    serialized_noise = NoiseHandler._from_dict(noise._to_dict())
+    assert noise == serialized_noise
+
+    rand_confusion = torch.rand(4, 4)
+    rand_confusion = rand_confusion / rand_confusion.sum(dim=1, keepdim=True)
+    noise = NoiseHandler(
+        protocol=NoiseProtocol.READOUT.CORRELATED,
+        options={"seed": 0, "confusion_matrix": rand_confusion},
+    )
     serialized_noise = NoiseHandler._from_dict(noise._to_dict())
     assert noise == serialized_noise
 
@@ -150,10 +193,21 @@ def test_serialization() -> None:
         [NoiseProtocol.DIGITAL.BITFLIP, NoiseProtocol.DIGITAL.PHASEFLIP],
     ],
 )
-def test_append(noise_config: NoiseProtocol | list[NoiseProtocol]) -> None:
-    noise = NoiseHandler(protocol=NoiseProtocol.READOUT)
+@pytest.mark.parametrize(
+    "initial_noise",
+    [
+        NoiseHandler(protocol=NoiseProtocol.READOUT.INDEPENDENT),
+        NoiseHandler(protocol=NoiseProtocol.READOUT.CORRELATED, options=torch.rand((4, 4))),
+    ],
+)
+def test_append(
+    initial_noise: NoiseHandler, noise_config: NoiseProtocol | list[NoiseProtocol]
+) -> None:
     options = {"error_probability": 0.1}
     with pytest.raises(ValueError):
-        noise.append(NoiseHandler(noise_config, options))
+        initial_noise.append(NoiseHandler(noise_config, options))
+    with pytest.raises(ValueError):
+        initial_noise.readout_independent(options)
+
     with pytest.raises(ValueError):
-        noise.readout(options)
+        initial_noise.readout_correlated({"confusion_matrix": torch.rand(4, 4)})