Simultaneous readout (#4307)

Adds an experiment for simultaneous (parallel) readout.

This experiment (estimate_parallel_single_qubit_readout_errors)  initializes a number of bitstrings(trials) and initializes the qubits to those bitstrings and immediately measures them It then estimates zero state and one state errors averaged over all of the trials where the given qubit was initialized to |0> or |1> respectively.  

This change also collapses the isolated single qubit readout experiment to use this function as well.

cirq-core/cirq/__init__.py

@@ -94,6 +94,7 @@
 )
 
 from cirq.experiments import (
+    estimate_parallel_single_qubit_readout_errors,
     estimate_single_qubit_readout_errors,
     hog_score_xeb_fidelity_from_probabilities,
     least_squares_xeb_fidelity_from_expectations,

cirq-core/cirq/experiments/__init__.py

@@ -72,6 +72,7 @@
 )
 
 from cirq.experiments.single_qubit_readout_calibration import (
+    estimate_parallel_single_qubit_readout_errors,
     estimate_single_qubit_readout_errors,
     SingleQubitReadoutCalibrationResult,
 )

cirq-core/cirq/experiments/single_qubit_readout_calibration.py

@@ -11,21 +11,20 @@
 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 # See the License for the specific language governing permissions and
 # limitations under the License.
-
-from typing import Any, Dict, Iterable, TYPE_CHECKING
-
+"""Single qubit readout experiments using parallel or isolated statistics."""
 import dataclasses
 import time
+from typing import Any, Dict, Iterable, List, Optional, TYPE_CHECKING
 
+import sympy
 import numpy as np
-
-from cirq import circuits, ops
+from cirq import circuits, ops, study
 
 if TYPE_CHECKING:
     import cirq
 
 
-@dataclasses.dataclass(frozen=True)
+@dataclasses.dataclass
 class SingleQubitReadoutCalibrationResult:
     """Result of estimating single qubit readout error.
 
@@ -96,19 +95,136 @@ def estimate_single_qubit_readout_errors(
         the probabilities. Also stores a timestamp indicating the time when
         data was finished being collected from the sampler.
     """
+    num_qubits = len(list(qubits))
+    return estimate_parallel_single_qubit_readout_errors(
+        sampler=sampler,
+        qubits=qubits,
+        repetitions=repetitions,
+        trials=2,
+        bit_strings=np.array([[0] * num_qubits, [1] * num_qubits]),
+    )
+
+
+def estimate_parallel_single_qubit_readout_errors(
+    sampler: 'cirq.Sampler',
+    *,
+    qubits: Iterable['cirq.Qid'],
+    trials: int = 20,
+    repetitions: int = 1000,
+    trials_per_batch: Optional[int] = None,
+    bit_strings: np.ndarray = None,
+) -> SingleQubitReadoutCalibrationResult:
+    """Estimate single qubit readout error using parallel operations.
+
+    For each trial, prepare and then measure a random computational basis
+    bitstring on qubits using gates in parallel.
+    Returns a SingleQubitReadoutCalibrationResult which can be used to
+    compute readout errors for each qubit.
+
+    Args:
+        sampler: The `cirq.Sampler` used to run the circuits.
+        qubits: The qubits being tested.
+        repetitions: The number of measurement repetitions to perform for
+            each trial.
+        trials: The number of bitstrings to prepare.
+        trials_per_batch:  If provided, split the experiment into batches
+            with this number of trials in each batch.
+        bit_strings: Optional numpy array of shape (trials, qubits) where the
+            first dimension is the number of the trial and the second
+            dimension is the qubit (ordered by the qubit order from
+            the qubits parameter).  Each value should be a 0 or 1 which
+            specifies which state the qubit should be prepared into during
+            that trial.  If not provided, the function will generate random
+            bit strings for you.
+
+    Returns:
+        A SingleQubitReadoutCalibrationResult storing the readout error
+        probabilities as well as the number of repetitions used to estimate
+        the probabilities. Also stores a timestamp indicating the time when
+        data was finished being collected from the sampler.  Note that,
+        if there did not exist a trial where a given qubit was set to |0〉,
+        the zero-state error will be set to `nan` (not a number).  Likewise
+        for qubits with no |1〉trial and one-state error.
+    """
     qubits = list(qubits)
 
-    zeros_circuit = circuits.Circuit(ops.measure_each(*qubits, key_func=repr))
-    ones_circuit = circuits.Circuit(
-        ops.X.on_each(*qubits), ops.measure_each(*qubits, key_func=repr)
+    if trials <= 0:
+        raise ValueError("Must provide non-zero trials for readout calibration.")
+    if repetitions <= 0:
+        raise ValueError("Must provide non-zero repetition for readout calibration.")
+    if bit_strings is None:
+        bit_strings = np.random.randint(0, 2, size=(trials, len(qubits)))
+    else:
+        if not hasattr(bit_strings, 'shape') or bit_strings.shape != (trials, len(qubits)):
+            raise ValueError(
+                'bit_strings must be numpy array '
+                f'of shape (trials, qubits) ({trials}, {len(qubits)}) '
+                f"but was {bit_strings.shape if hasattr(bit_strings, 'shape') else None}"
+            )
+        if not np.all((bit_strings == 0) | (bit_strings == 1)):
+            raise ValueError('bit_strings values must be all 0 or 1')
+    if trials_per_batch is None:
+        trials_per_batch = trials
+    if trials_per_batch <= 0:
+        raise ValueError("Must provide non-zero trials_per_batch for readout calibration.")
+
+    all_sweeps: List[study.Sweepable] = []
+    num_batches = (trials + trials_per_batch - 1) // trials_per_batch
+
+    # Initialize circuits
+    flip_symbols = sympy.symbols(f'flip_0:{len(qubits)}')
+    flip_circuit = circuits.Circuit(
+        [ops.X(q) ** s for q, s in zip(qubits, flip_symbols)],
+        [ops.measure_each(*qubits, key_func=repr)],
     )
-
-    zeros_result = sampler.run(zeros_circuit, repetitions=repetitions)
-    ones_result = sampler.run(ones_circuit, repetitions=repetitions)
+    all_circuits = [flip_circuit] * num_batches
+
+    # Initialize sweeps
+    for batch in range(num_batches):
+        single_sweeps = []
+        for qubit_idx in range(len(qubits)):
+            trial_range = range(
+                batch * trials_per_batch, min((batch + 1) * trials_per_batch, trials)
+            )
+            single_sweeps.append(
+                study.Points(
+                    key=f'flip_{qubit_idx}',
+                    points=[bit_strings[bit][qubit_idx] for bit in trial_range],
+                )
+            )
+        total_sweeps = study.Zip(*single_sweeps)
+        all_sweeps.append(total_sweeps)
+
+    # Execute circuits
+    results = sampler.run_batch(all_circuits, all_sweeps, repetitions=repetitions)
     timestamp = time.time()
 
-    zero_state_errors = {q: np.mean(zeros_result.measurements[repr(q)]) for q in qubits}
-    one_state_errors = {q: 1 - np.mean(ones_result.measurements[repr(q)]) for q in qubits}
+    # Analyze results
+    zero_state_trials = np.zeros((1, len(qubits)))
+    one_state_trials = np.zeros((1, len(qubits)))
+    zero_state_totals = np.zeros((1, len(qubits)))
+    one_state_totals = np.zeros((1, len(qubits)))
+    for batch_result in results:
+        for trial_idx, trial_result in enumerate(batch_result):
+            all_measurements = trial_result.data[[repr(x) for x in qubits]].to_numpy()
+            sample_counts = np.einsum('ij->j', all_measurements)
+            zero_state_trials += sample_counts * (1 - bit_strings[trial_idx])
+            zero_state_totals += repetitions * (1 - bit_strings[trial_idx])
+            one_state_trials += (repetitions - sample_counts) * bit_strings[trial_idx]
+            one_state_totals += repetitions * bit_strings[trial_idx]
+
+    zero_state_errors = {
+        q: zero_state_trials[0][qubit_idx] / zero_state_totals[0][qubit_idx]
+        if zero_state_totals[0][qubit_idx] > 0
+        else np.nan
+        for qubit_idx, q in enumerate(qubits)
+    }
+    one_state_errors = {
+        q: one_state_trials[0][qubit_idx] / one_state_totals[0][qubit_idx]
+        if one_state_totals[0][qubit_idx] > 0
+        else np.nan
+        for qubit_idx, q in enumerate(qubits)
+    }
 
     return SingleQubitReadoutCalibrationResult(
         zero_state_errors=zero_state_errors,

cirq-core/cirq/experiments/single_qubit_readout_calibration_test.py

@@ -11,14 +11,24 @@
 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 # See the License for the specific language governing permissions and
 # limitations under the License.
-
 from typing import List
+import pytest
 
 import numpy as np
 
 import cirq
 
 
+def test_single_qubit_readout_result_repr():
+    result = cirq.experiments.SingleQubitReadoutCalibrationResult(
+        zero_state_errors={cirq.LineQubit(0): 0.1},
+        one_state_errors={cirq.LineQubit(0): 0.2},
+        repetitions=1000,
+        timestamp=0.3,
+    )
+    cirq.testing.assert_equivalent_repr(result)
+
+
 class NoisySingleQubitReadoutSampler(cirq.Sampler):
     def __init__(self, p0: float, p1: float, seed: 'cirq.RANDOM_STATE_OR_SEED_LIKE' = None):
         """Sampler that flips some bits upon readout.
@@ -79,11 +89,129 @@ def test_estimate_single_qubit_readout_errors_with_noise():
     assert isinstance(result.timestamp, float)
 
 
-def test_single_qubit_readout_calibration_result_repr():
-    result = cirq.experiments.SingleQubitReadoutCalibrationResult(
-        zero_state_errors={cirq.LineQubit(0): 0.1},
-        one_state_errors={cirq.LineQubit(0): 0.2},
-        repetitions=1000,
-        timestamp=0.3,
+def test_estimate_parallel_readout_errors_no_noise():
+    qubits = cirq.LineQubit.range(10)
+    sampler = cirq.Simulator()
+    repetitions = 1000
+    result = cirq.estimate_parallel_single_qubit_readout_errors(
+        sampler, qubits=qubits, repetitions=repetitions
     )
-    cirq.testing.assert_equivalent_repr(result)
+    assert result.zero_state_errors == {q: 0 for q in qubits}
+    assert result.one_state_errors == {q: 0 for q in qubits}
+    assert result.repetitions == repetitions
+    assert isinstance(result.timestamp, float)
+
+
+def test_estimate_parallel_readout_errors_all_zeros():
+    qubits = cirq.LineQubit.range(10)
+    sampler = cirq.ZerosSampler()
+    repetitions = 1000
+    result = cirq.estimate_parallel_single_qubit_readout_errors(
+        sampler, qubits=qubits, repetitions=repetitions
+    )
+    assert result.zero_state_errors == {q: 0 for q in qubits}
+    assert result.one_state_errors == {q: 1 for q in qubits}
+    assert result.repetitions == repetitions
+    assert isinstance(result.timestamp, float)
+
+
+def test_estimate_parallel_readout_errors_bad_bit_string():
+    qubits = cirq.LineQubit.range(4)
+    with pytest.raises(ValueError, match='but was None'):
+        _ = cirq.estimate_parallel_single_qubit_readout_errors(
+            cirq.ZerosSampler(),
+            qubits=qubits,
+            repetitions=1000,
+            trials=35,
+            trials_per_batch=10,
+            bit_strings=[[1] * 4],
+        )
+    with pytest.raises(ValueError, match='0 or 1'):
+        _ = cirq.estimate_parallel_single_qubit_readout_errors(
+            cirq.ZerosSampler(),
+            qubits=qubits,
+            repetitions=1000,
+            trials=2,
+            bit_strings=np.array([[12, 47, 2, -4], [0.1, 7, 0, 0]]),
+        )
+
+
+def test_estimate_parallel_readout_errors_zero_reps():
+    qubits = cirq.LineQubit.range(10)
+    with pytest.raises(ValueError, match='non-zero repetition'):
+        _ = cirq.estimate_parallel_single_qubit_readout_errors(
+            cirq.ZerosSampler(),
+            qubits=qubits,
+            repetitions=0,
+            trials=35,
+            trials_per_batch=10,
+        )
+
+
+def test_estimate_parallel_readout_errors_zero_trials():
+    qubits = cirq.LineQubit.range(10)
+    with pytest.raises(ValueError, match='non-zero trials'):
+        _ = cirq.estimate_parallel_single_qubit_readout_errors(
+            cirq.ZerosSampler(),
+            qubits=qubits,
+            repetitions=1000,
+            trials=0,
+            trials_per_batch=10,
+        )
+
+
+def test_estimate_parallel_readout_errors_zero_batch():
+    qubits = cirq.LineQubit.range(10)
+    with pytest.raises(ValueError, match='non-zero trials_per_batch'):
+        _ = cirq.estimate_parallel_single_qubit_readout_errors(
+            cirq.ZerosSampler(),
+            qubits=qubits,
+            repetitions=1000,
+            trials=10,
+            trials_per_batch=0,
+        )
+
+
+def test_estimate_parallel_readout_errors_batching():
+    qubits = cirq.LineQubit.range(10)
+    sampler = cirq.ZerosSampler()
+    repetitions = 1000
+    result = cirq.estimate_parallel_single_qubit_readout_errors(
+        sampler, qubits=qubits, repetitions=repetitions, trials=45, trials_per_batch=10
+    )
+    assert result.zero_state_errors == {q: 0.0 for q in qubits}
+    assert result.one_state_errors == {q: 1.0 for q in qubits}
+    assert result.repetitions == repetitions
+    assert isinstance(result.timestamp, float)
+
+
+def test_estimate_parallel_readout_errors_with_noise():
+    qubits = cirq.LineQubit.range(5)
+    sampler = NoisySingleQubitReadoutSampler(p0=0.1, p1=0.2, seed=1234)
+    repetitions = 1000
+    result = cirq.estimate_parallel_single_qubit_readout_errors(
+        sampler, qubits=qubits, repetitions=repetitions, trials=40
+    )
+    for error in result.one_state_errors.values():
+        assert 0.17 < error < 0.23
+    for error in result.zero_state_errors.values():
+        assert 0.07 < error < 0.13
+    assert result.repetitions == repetitions
+    assert isinstance(result.timestamp, float)
+
+
+def test_estimate_parallel_readout_errors_missing_qubits():
+    qubits = cirq.LineQubit.range(4)
+
+    result = cirq.estimate_parallel_single_qubit_readout_errors(
+        cirq.ZerosSampler(),
+        qubits=qubits,
+        repetitions=2000,
+        trials=1,
+        bit_strings=np.array([[0] * 4]),
+    )
+    assert result.zero_state_errors == {q: 0 for q in qubits}
+    # Trial did not include a one-state
+    assert all(np.isnan(result.one_state_errors[q]) for q in qubits)
+    assert result.repetitions == 2000
+    assert isinstance(result.timestamp, float)

cirq-core/cirq/protocols/json_test_data/SingleQubitReadoutCalibrationResult.repr

@@ -1 +1 @@
-cirq.experiments.SingleQubitReadoutCalibrationResult(zero_state_errors={cirq.LineQubit(0): 0.1}, one_state_errors={cirq.LineQubit(0): 0.2}, repetitions=1000, timestamp=0.3)
+cirq.experiments.SingleQubitReadoutCalibrationResult(zero_state_errors={cirq.LineQubit(0): 0.1}, one_state_errors={cirq.LineQubit(0): 0.2}, repetitions=1000, timestamp=0.3)