Add exponential fitting to RandomizedBenchMarkResult

cirq-core/cirq/experiments/qubit_characterizations.py

@@ -17,6 +17,7 @@
 
 from typing import Any, cast, Iterator, List, Optional, Sequence, Tuple, TYPE_CHECKING
 import numpy as np
+from scipy.optimize import curve_fit
 
 from matplotlib import pyplot as plt
 
@@ -92,13 +93,44 @@ def plot(self, ax: Optional[plt.Axes] = None, **plot_kwargs: Any) -> plt.Axes:
             fig, ax = plt.subplots(1, 1, figsize=(8, 8))  # pragma: no cover
             ax = cast(plt.Axes, ax)  # pragma: no cover
         ax.set_ylim((0.0, 1.0))  # pragma: no cover
-        ax.plot(self._num_cfds_seq, self._gnd_state_probs, 'ro-', **plot_kwargs)
+        ax.plot(self._num_cfds_seq, self._gnd_state_probs, 'ro', label='data', **plot_kwargs)
+        x = np.linspace(self._num_cfds_seq[0], self._num_cfds_seq[-1], 100)
+        fit = self._fit_exponential()
+        ax.plot(x, fit[0][0] * fit[0][2] ** x + fit[0][1], '--k', label='fit')
+        ax.legend(loc='upper right')
         ax.set_xlabel(r"Number of Cliffords")
         ax.set_ylabel('Ground State Probability')
         if show_plot:
             fig.show()
         return ax
 
+    def pauli_error(self) -> float:
+        """Returns the Pauli error inferred from randomized benchmarking.
+
+        If sequence fidelity (F) decays with number of gates (m) as
+
+        F = A p^m + B
+
+        where 0 < p < 1, then the Pauli error (r_p) is given by
+
+        r_p = (1 - 1/d^2) * (1 - p)
+
+        Where d = 2^num_qubits is the Hilbert space dimension.
+        """
+        fit = self._fit_exponential()
+        p = fit[0][2]
+        return (1.0 - 1.0 / 4.0) * (1.0 - p)
+
+    def _fit_exponential(self) -> np.ndarray:
+        exp_fit = lambda x, A, B, p: A * p**x + B
+        return curve_fit(
+            exp_fit,
+            self._num_cfds_seq,
+            self._gnd_state_probs,
+            [0.5, 0.5, 1.0 - 1e-3],
+            bounds=([0, 0.25, 0], [0.5, 0.75, 1]),
+        )
+
 
 class TomographyResult:
     """Results from a state tomography experiment."""

cirq-core/cirq/experiments/qubit_characterizations_test.py

@@ -85,12 +85,11 @@ def test_single_qubit_randomized_benchmarking():
     # sequences is always unity.
     simulator = sim.Simulator()
     qubit = GridQubit(0, 0)
-    num_cfds = range(5, 20, 5)
-    results = single_qubit_randomized_benchmarking(
-        simulator, qubit, num_clifford_range=num_cfds, repetitions=100
-    )
+    num_cfds = tuple(np.logspace(np.log10(5), 3, 5, dtype=int))
+    results = single_qubit_randomized_benchmarking(simulator, qubit, num_clifford_range=num_cfds)
     g_pops = np.asarray(results.data)[:, 1]
     assert np.isclose(np.mean(g_pops), 1.0)
+    assert np.isclose(results.pauli_error(), 0.0, atol=1e-7)  # warning is expected
 
 
 def test_two_qubit_randomized_benchmarking():