[Obs] 4.4 - Readout error mitigation (quantumlib#4323)

Add a facility to take readout calibrations by measuring the <Z> observables under the empty circuit with readout symmetrization turned on. This can be plumbed through so other observables measured with `measure_observables` can be corrected. Note that this all only works for symmetric readout error enforced through readout_symmetrization.

Part of quantumlib#3647

cirq-core/cirq/work/__init__.py

@@ -39,6 +39,9 @@
     RepetitionsStoppingCriteria,
     measure_grouped_settings,
 )
+from cirq.work.observable_readout_calibration import (
+    calibrate_readout_error,
+)
 from cirq.work.sampler import (
     Sampler,
 )

cirq-core/cirq/work/observable_measurement.py

@@ -16,11 +16,10 @@
 import dataclasses
 import itertools
 import warnings
-from typing import Iterable, Dict, List, Tuple, TYPE_CHECKING, Set, Sequence
+from typing import Optional, Iterable, Dict, List, Tuple, TYPE_CHECKING, Set, Sequence
 
 import numpy as np
 import sympy
-
 from cirq import circuits, study, ops, value
 from cirq._doc import document
 from cirq.protocols import json_serializable_dataclass
@@ -371,6 +370,7 @@ def measure_grouped_settings(
     *,
     readout_symmetrization: bool = False,
     circuit_sweep: 'cirq.study.sweepable.SweepLike' = None,
+    readout_calibrations: Optional[BitstringAccumulator] = None,
 ) -> List[BitstringAccumulator]:
     """Measure a suite of grouped InitObsSetting settings.
 
@@ -398,7 +398,11 @@ def measure_grouped_settings(
         circuit_sweep: Additional parameter sweeps for parameters contained
             in `circuit`. The total sweep is the product of the circuit sweep
             with parameter settings for the single-qubit basis-change rotations.
+        readout_calibrations: The result of `calibrate_readout_error`.
     """
+    if readout_calibrations is not None and not readout_symmetrization:
+        raise ValueError("Readout calibration only works if `readout_symmetrization` is enabled.")
+
     qubits = sorted({q for ms in grouped_settings.keys() for q in ms.init_state.qubits})
     qubit_to_index = {q: i for i, q in enumerate(qubits)}
 
@@ -427,6 +431,7 @@ def measure_grouped_settings(
             meas_spec=meas_spec,
             simul_settings=grouped_settings[max_setting],
             qubit_to_index=qubit_to_index,
+            readout_calibration=readout_calibrations,
         )
         accumulators[meas_spec] = accumulator
         meas_specs_todo += [meas_spec]

cirq-core/cirq/work/observable_measurement_test.py

@@ -401,3 +401,32 @@ def test_measure_grouped_settings(with_circuit_sweep):
         else:
             (result,) = results  # one group
             assert result.means() == [coef]
+
+
+def _get_some_grouped_settings():
+    qubits = cirq.LineQubit.range(2)
+    q0, q1 = qubits
+    terms = [
+        cirq.X(q0),
+        cirq.Y(q1),
+    ]
+    settings = list(cirq.work.observables_to_settings(terms, qubits))
+    grouped_settings = cirq.work.group_settings_greedy(settings)
+    return grouped_settings, qubits
+
+
+def test_measure_grouped_settings_calibration_validation():
+    dummy_ro_calib = _MockBitstringAccumulator()
+    grouped_settings, qubits = _get_some_grouped_settings()
+
+    with pytest.raises(
+        ValueError, match=r'Readout calibration only works if `readout_symmetrization` is enabled'
+    ):
+        cw.measure_grouped_settings(
+            circuit=cirq.Circuit(cirq.I.on_each(*qubits)),
+            grouped_settings=grouped_settings,
+            sampler=cirq.Simulator(),
+            stopping_criteria=cw.RepetitionsStoppingCriteria(10_000),
+            readout_calibrations=dummy_ro_calib,
+            readout_symmetrization=False,  # no-no!
+        )

cirq-core/cirq/work/observable_readout_calibration.py

@@ -0,0 +1,52 @@
+import dataclasses
+from typing import Union, Iterable, TYPE_CHECKING
+
+from cirq import circuits, study, ops
+from cirq.work.observable_measurement import measure_grouped_settings, StoppingCriteria
+from cirq.work.observable_settings import (
+    InitObsSetting,
+    zeros_state,
+)
+
+if TYPE_CHECKING:
+    import cirq
+
+
+def calibrate_readout_error(
+    qubits: Iterable[ops.Qid],
+    sampler: Union['cirq.Simulator', 'cirq.Sampler'],
+    stopping_criteria: StoppingCriteria,
+):
+    # We know there won't be any fancy sweeps or observables so we can
+    # get away with more repetitions per job
+    stopping_criteria = dataclasses.replace(stopping_criteria, repetitions_per_chunk=100_000)
+
+    # Simultaneous readout characterization:
+    # We can measure all qubits simultaneously (i.e. _max_setting is ZZZ..ZZ
+    # for all qubits). We will extract individual qubit quantities, so there
+    # are `n_qubits` InitObsSetting, each responsible for one <Z>.
+    #
+    # Readout symmetrization means we just need to measure the "identity"
+    # circuit. In reality, this corresponds to measuring I for half the time
+    # and X for the other half.
+    init_state = zeros_state(qubits)
+    max_setting = InitObsSetting(
+        init_state=init_state, observable=ops.PauliString({q: ops.Z for q in qubits})
+    )
+    grouped_settings = {
+        max_setting: [
+            InitObsSetting(init_state=init_state, observable=ops.PauliString({q: ops.Z}))
+            for q in qubits
+        ]
+    }
+
+    results = measure_grouped_settings(
+        circuit=circuits.Circuit(),
+        grouped_settings=grouped_settings,
+        sampler=sampler,
+        stopping_criteria=stopping_criteria,
+        circuit_sweep=study.UnitSweep,
+        readout_symmetrization=True,
+    )
+    (result,) = list(results)
+    return result

cirq-core/cirq/work/observable_readout_calibration_test.py

@@ -0,0 +1,48 @@
+from typing import Sequence
+
+import cirq
+import cirq.work as cw
+import numpy as np
+
+
+class DepolarizingWithDampedReadoutNoiseModel(cirq.NoiseModel):
+    """This simulates asymmetric readout error.
+
+    The noise is structured so the T1 decay is applied, then the readout bitflip, then measurement.
+    If a circuit contains measurements, they must be in moments that don't also contain gates.
+    """
+
+    def __init__(self, depol_prob: float, bitflip_prob: float, decay_prob: float):
+        self.qubit_noise_gate = cirq.DepolarizingChannel(depol_prob)
+        self.readout_noise_gate = cirq.BitFlipChannel(bitflip_prob)
+        self.readout_decay_gate = cirq.AmplitudeDampingChannel(decay_prob)
+
+    def noisy_moment(self, moment: 'cirq.Moment', system_qubits: Sequence['cirq.Qid']):
+        if cirq.devices.noise_model.validate_all_measurements(moment):
+            return [
+                cirq.Moment(self.readout_decay_gate(q) for q in system_qubits),
+                cirq.Moment(self.readout_noise_gate(q) for q in system_qubits),
+                moment,
+            ]
+        else:
+            return [moment, cirq.Moment(self.qubit_noise_gate(q) for q in system_qubits)]
+
+
+def test_calibrate_readout_error():
+    sampler = cirq.DensityMatrixSimulator(
+        noise=DepolarizingWithDampedReadoutNoiseModel(
+            depol_prob=1e-3,
+            bitflip_prob=0.03,
+            decay_prob=0.03,
+        ),
+        seed=10,
+    )
+    readout_calibration = cw.calibrate_readout_error(
+        qubits=cirq.LineQubit.range(2),
+        sampler=sampler,
+        stopping_criteria=cw.RepetitionsStoppingCriteria(100_000),
+    )
+    means = readout_calibration.means()
+    assert len(means) == 2, 'Two qubits'
+    assert np.all(means > 0.89)
+    assert np.all(means < 0.91)