From 2ce4ed124f4f857ae895f6bbffee97fe973c9a5e Mon Sep 17 00:00:00 2001
From: eliottrosenberg <61400172+eliottrosenberg@users.noreply.github.com>
Date: Mon, 18 Dec 2023 15:44:37 -0800
Subject: [PATCH] Add exponential fitting to RandomizedBenchMarkResult (#6385)

---
 .../experiments/qubit_characterizations.py    | 34 ++++++++++++++++++-
 .../qubit_characterizations_test.py           |  7 ++--
 2 files changed, 36 insertions(+), 5 deletions(-)

diff --git a/cirq-core/cirq/experiments/qubit_characterizations.py b/cirq-core/cirq/experiments/qubit_characterizations.py
index ed12b311e22..aa39631829f 100644
--- a/cirq-core/cirq/experiments/qubit_characterizations.py
+++ b/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)
+        opt_params, _ = self._fit_exponential()
+        ax.plot(x, opt_params[0] * opt_params[2] ** x + opt_params[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:
+        r"""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^N_Q$ is the Hilbert space dimension and $N_Q$ is the number of qubits.
+        """
+        opt_params, _ = self._fit_exponential()
+        p = opt_params[2]
+        return (1.0 - 1.0 / 4.0) * (1.0 - p)
+
+    def _fit_exponential(self) -> Tuple[np.ndarray, np.ndarray]:
+        exp_fit = lambda x, A, B, p: A * p**x + B
+        return curve_fit(
+            f=exp_fit,
+            xdata=self._num_cfds_seq,
+            ydata=self._gnd_state_probs,
+            p0=[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."""
diff --git a/cirq-core/cirq/experiments/qubit_characterizations_test.py b/cirq-core/cirq/experiments/qubit_characterizations_test.py
index baa82c02ffc..e8f09694e42 100644
--- a/cirq-core/cirq/experiments/qubit_characterizations_test.py
+++ b/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():