Interleaved RB

qiskit_experiments/analysis/curve_fitting.py

@@ -293,7 +293,7 @@ def process_multi_curve_data(
     ydata = np.zeros(size, dtype=float)
     ydata_var = np.zeros(size, dtype=float)
 
-    for i, datum in enumerate(filter_data):
+    for i, datum in enumerate(filtered_data):
         metadata = datum["metadata"]
         series[i] = metadata[series_key]
         xdata[i] = metadata[x_key]

qiskit_experiments/analysis/data_processing.py

@@ -120,6 +120,45 @@ def mean_xy_data(
     raise QiskitError(f"Unsupported method {method}")
 
 
+def multi_mean_xy_data(
+    series: np.ndarray,
+    xdata: np.ndarray,
+    ydata: np.ndarray,
+    sigma: Optional[np.ndarray] = None,
+    method: str = "sample",
+):
+    r"""Return (series, x, y_mean, sigma) data.
+    Performs `mean_xy_data` for each series
+    and returns the concatenated results
+    """
+    series_vals = np.unique(series)
+
+    series_means = []
+    xdata_means = []
+    ydata_means = []
+    sigma_means = []
+
+    # Get x, y, sigma data for series and process mean data
+    for i in series_vals:
+        idxs = series == series_vals[i]
+        sigma_i = sigma[idxs] if sigma is not None else None
+        x_mean, y_mean, sigma_mean = mean_xy_data(
+            xdata[idxs], ydata[idxs], sigma=sigma_i, method=method
+        )
+        series_means.append(i * np.ones(x_mean.size, dtype=int))
+        xdata_means.append(x_mean)
+        ydata_means.append(y_mean)
+        sigma_means.append(sigma_mean)
+
+    # Concatenate lists
+    return (
+        np.concatenate(series_means),
+        np.concatenate(xdata_means),
+        np.concatenate(ydata_means),
+        np.concatenate(sigma_means),
+    )
+
+
 def level2_probability(data: Dict[str, any], outcome: str) -> Tuple[float]:
     """Return the outcome probability mean and variance.
 

qiskit_experiments/randomized_benchmarking/__init__.py

@@ -13,4 +13,6 @@
 """Randomized Benchmarking Experiment Classes."""
 
 from .rb_experiment import RBExperiment
+from .interleaved_rb_experiment import InterleavedRBExperiment
 from .rb_analysis import RBAnalysis
+from .interleaved_rb_analysis import InterleavedRBAnalysis

qiskit_experiments/randomized_benchmarking/interleaved_rb_analysis.py

@@ -0,0 +1,176 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2021.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+"""
+Interleaved RB analysis class.
+"""
+from typing import Optional, List
+import numpy as np
+from qiskit_experiments.analysis.curve_fitting import (
+    process_multi_curve_data,
+    multi_curve_fit,
+)
+
+from qiskit_experiments.analysis.plotting import plot_curve_fit, plot_scatter, plot_errorbar
+from qiskit_experiments.analysis.data_processing import (
+    level2_probability,
+    multi_mean_xy_data,
+)
+from .rb_analysis import RBAnalysis
+
+try:
+    from matplotlib import pyplot as plt
+
+    HAS_MATPLOTLIB = True
+except ImportError:
+    HAS_MATPLOTLIB = False
+
+
+class InterleavedRBAnalysis(RBAnalysis):
+    r"""Interleaved RB Analysis class.
+    According to the paper: "Efficient measurement of quantum gate
+    error by interleaved randomized benchmarking" (arXiv:1203.4550)
+
+    The epc estimate is obtained using the equation
+    :math:`r_{\mathcal{C}}^{\text{est}}=
+    \frac{\left(d-1\right)\left(1-p_{\overline{\mathcal{C}}}/p\right)}{d}`
+
+    The error bounds are given by
+    :math:`E=\min\left\{ \begin{array}{c}
+    \frac{\left(d-1\right)\left[\left|p-p_{\overline{\mathcal{C}}}/p\right|+\left(1-p\right)\right]}{d}\\
+    \frac{2\left(d^{2}-1\right)\left(1-p\right)}{pd^{2}}+\frac{4\sqrt{1-p}\sqrt{d^{2}-1}}{p}
+    \end{array}\right.`
+    """
+
+    # pylint: disable=invalid-name
+    def _run_analysis(
+        self,
+        experiment_data,
+        p0: Optional[List[float]] = None,
+        plot: bool = True,
+        ax: Optional["AxesSubplot"] = None,
+    ):
+        def data_processor(datum):
+            return level2_probability(datum, datum["metadata"]["ylabel"])
+
+        num_qubits = len(experiment_data.data[0]["metadata"]["qubits"])
+        series, x, y, sigma = process_multi_curve_data(experiment_data.data, data_processor)
+        series, xdata, ydata, ydata_sigma = multi_mean_xy_data(series, x, y, sigma)
+
+        def fit_fun_standard(x, a, alpha_std, _, b):
+            return a * alpha_std ** x + b
+
+        def fit_fun_interleaved(x, a, _, alpha_int, b):
+            return a * alpha_int ** x + b
+
+        std_idx = series == 0
+        std_xdata = xdata[std_idx]
+        std_ydata = ydata[std_idx]
+        std_ydata_sigma = ydata_sigma[std_idx]
+        p0_std = self._p0(std_xdata, std_ydata, num_qubits)
+
+        int_idx = series == 1
+        int_xdata = xdata[int_idx]
+        int_ydata = ydata[int_idx]
+        int_ydata_sigma = ydata_sigma[int_idx]
+        p0_int = self._p0(int_xdata, int_ydata, num_qubits)
+
+        p0 = (
+            np.mean([p0_std[0], p0_int[0]]),
+            p0_std[1],
+            p0_int[1],
+            np.mean([p0_std[2], p0_int[2]]),
+        )
+
+        analysis_result = multi_curve_fit(
+            [fit_fun_standard, fit_fun_interleaved],
+            series,
+            xdata,
+            ydata,
+            p0,
+            ydata_sigma,
+            bounds=([0, 0, 0, 0], [1, 1, 1, 1]),
+        )
+
+        # Add EPC data
+        nrb = 2 ** num_qubits
+        scale = (nrb - 1) / (2 ** nrb)
+        _, alpha, alpha_c, _ = analysis_result["popt"]
+        _, alpha_err, alpha_c_err, _ = analysis_result["popt_err"]
+
+        # Calculate epc_est (=r_c^est) - Eq. (4):
+        epc_est = scale * (1 - alpha_c / alpha)
+        # Calculate the systematic error bounds - Eq. (5):
+        systematic_err_1 = scale * (abs(alpha - alpha_c / alpha) + (1 - alpha))
+        systematic_err_2 = (
+            2 * (nrb * nrb - 1) * (1 - alpha) / (alpha * nrb * nrb)
+            + 4 * (np.sqrt(1 - alpha)) * (np.sqrt(nrb * nrb - 1)) / alpha
+        )
+        systematic_err = min(systematic_err_1, systematic_err_2)
+        systematic_err_l = epc_est - systematic_err
+        systematic_err_r = epc_est + systematic_err
+
+        alpha_err_sq = (alpha_err / alpha) ** 2
+        alpha_c_err_sq = (alpha_c_err / alpha_c) ** 2
+        epc_est_err = (
+            ((nrb - 1) / nrb) * (alpha_c / alpha) * (np.sqrt(alpha_err_sq + alpha_c_err_sq))
+        )
+
+        analysis_result["EPC"] = epc_est
+        analysis_result["EPC_err"] = epc_est_err
+        analysis_result["systematic_err"] = systematic_err
+        analysis_result["systematic_err_L"] = systematic_err_l
+        analysis_result["systematic_err_R"] = systematic_err_r
+        analysis_result["plabels"] = ["A", "alpha", "alpha_c", "B"]
+
+        if plot:
+            ax = plot_curve_fit(fit_fun_standard, analysis_result, ax=ax)
+            ax = plot_curve_fit(fit_fun_interleaved, analysis_result, ax=ax)
+            ax = plot_scatter(std_xdata, std_ydata, ax=ax)
+            ax = plot_scatter(int_xdata, int_ydata, ax=ax)
+            ax = plot_errorbar(std_xdata, std_ydata, std_ydata_sigma, ax=ax)
+            ax = plot_errorbar(int_xdata, int_ydata, int_ydata_sigma, ax=ax)
+            self._format_plot(ax, analysis_result)
+            analysis_result.plt = plt
+
+        return analysis_result, None
+
+    @classmethod
+    def _format_plot(cls, ax, analysis_result, add_label=True):
+        """Format curve fit plot"""
+        # Formatting
+        ax.tick_params(labelsize=14)
+        ax.set_xlabel("Clifford Length", fontsize=16)
+        ax.set_ylabel("Ground State Population", fontsize=16)
+        ax.grid(True)
+
+        if add_label:
+            alpha = analysis_result["popt"][1]
+            alpha_c = analysis_result["popt"][2]
+            alpha_err = analysis_result["popt_err"][1]
+            alpha_c_err = analysis_result["popt_err"][2]
+            epc = analysis_result["EPC"]
+            epc_err = analysis_result["EPC_err"]
+            box_text = "\u03B1:{:.4f} \u00B1 {:.4f}".format(alpha, alpha_err)
+            box_text += "\n\u03B1_c:{:.4f} \u00B1 {:.4f}".format(alpha_c, alpha_c_err)
+            box_text += "\nEPC: {:.4f} \u00B1 {:.4f}".format(epc, epc_err)
+            bbox_props = dict(boxstyle="square,pad=0.3", fc="white", ec="black", lw=1)
+            ax.text(
+                0.6,
+                0.9,
+                box_text,
+                ha="center",
+                va="center",
+                size=14,
+                bbox=bbox_props,
+                transform=ax.transAxes,
+            )
+        return ax

qiskit_experiments/randomized_benchmarking/interleaved_rb_experiment.py

@@ -0,0 +1,93 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2021.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+"""
+Interleaved RB Experiment class.
+"""
+from typing import Union, Iterable, Optional, List
+
+from numpy.random import Generator
+
+from qiskit import QuantumCircuit
+from qiskit.circuit import Instruction
+from qiskit.quantum_info import Clifford, random_clifford
+
+from .rb_experiment import RBExperiment
+from .interleaved_rb_analysis import InterleavedRBAnalysis
+
+
+class InterleavedRBExperiment(RBExperiment):
+    """Interleaved RB Experiment class"""
+
+    # Analysis class for experiment
+    __analysis_class__ = InterleavedRBAnalysis
+
+    def __init__(
+        self,
+        interleaved_element: Union[QuantumCircuit, Instruction, Clifford],
+        qubits: Union[int, Iterable[int]],
+        lengths: Iterable[int],
+        num_samples: int = 1,
+        seed: Optional[Union[int, Generator]] = None,
+        full_sampling: bool = False,
+    ):
+        """Interleaved randomized benchmarking experiment
+        Args:
+            interleaved_element: the element to interleave,
+                    given either as a group element or as an instruction/circuit
+            qubits: the number of qubits or list of
+                    physical qubits for the experiment.
+            lengths: A list of RB sequences lengths.
+            num_samples: number of samples to generate for each
+                         sequence length
+            seed: Seed or generator object for random number
+                  generation. If None default_rng will be used.
+            full_sampling: If True all Cliffords are independently sampled for
+                           all lengths. If False for sample of lengths longer
+                           sequences are constructed by appending additional
+                           Clifford samples to shorter sequences.
+        """
+        self._interleaved_element = interleaved_element
+        super().__init__(qubits, lengths, num_samples, seed, full_sampling)
+
+    def _sample_circuits(self, lengths, seed=None):
+        circuits = []
+        for length in lengths if self._full_sampling else [lengths[-1]]:
+            elements = [random_clifford(self.num_qubits, seed=seed) for _ in range(length)]
+            element_lengths = [len(elements)] if self._full_sampling else lengths
+            std_circuits = self._generate_circuit(elements, element_lengths)
+            for circuit in std_circuits:
+                circuit.metadata["series"] = 0
+                circuit.metadata["series_name"] = "standard"
+            circuits += std_circuits
+
+            int_elements = self._interleave(elements)
+            int_elements_lengths = [length * 2 for length in element_lengths]
+            int_circuits = self._generate_circuit(int_elements, int_elements_lengths)
+            for circuit in int_circuits:
+                circuit.metadata["series"] = 1
+                circuit.metadata["series_name"] = "interleaved"
+                circuit.metadata["xval"] = circuit.metadata["xval"] // 2
+            circuits += int_circuits
+        return circuits
+
+    def _interleave(self, element_list: List) -> List:
+        """Interleaving the interleaved element inside the element list
+        Args:
+            element_list: The list of elements we add the interleaved element to
+        Returns:
+            The new list with the element interleaved
+        """
+        new_element_list = []
+        for element in element_list:
+            new_element_list.append(element)
+            new_element_list.append(self._interleaved_element)
+        return new_element_list

qiskit_experiments/randomized_benchmarking/rb_experiment.py

@@ -139,7 +139,7 @@ def _generate_circuit(self, elements: Iterable[Clifford], lengths: Iterable[int]
                 rb_circ.barrier(qubits)
                 rb_circ.metadata = {
                     "experiment_type": self._type,
-                    "xdata": current_length + 1,
+                    "xval": current_length + 1,
                     "ylabel": self.num_qubits * "0",
                     "group": "Clifford",
                     "qubits": self.physical_qubits,

test/test_curve_fitting.py

@@ -17,7 +17,11 @@
 from qiskit import QuantumCircuit, execute
 from qiskit.providers.basicaer import QasmSimulatorPy
 from qiskit_experiments.analysis.curve_fitting import curve_fit, multi_curve_fit, process_curve_data
-from qiskit_experiments.analysis.data_processing import level2_probability
+from qiskit_experiments.analysis.data_processing import (
+    level2_probability,
+    mean_xy_data,
+    multi_mean_xy_data,
+)
 
 
 class TestCurveFitting(QiskitTestCase):
@@ -108,3 +112,38 @@ def test_multi_curve_fit(self):
             [self.objective0, self.objective1], series, xdata, ydata, p0, sigma=sigma, bounds=bounds
         )
         self.assertTrue(abs(sol["popt"][0] - 0.5) < 0.05)
+
+    def test_mean_xy_data(self):
+        """Test mean_xy_data function"""
+        # pylint: disable=unbalanced-tuple-unpacking
+        x = np.array([1, 1, 1, 2, 2, 2, 2, 3, 3, 4, 5, 5, 5, 5])
+        y = np.array([1, 2, 3, 8, 10, 50, 60, 10, 11, 17, 10, 10, 10, 10])
+        x_mean, y_mean, y_sigma = mean_xy_data(x, y, method="sample")
+
+        expected_x_mean = np.array([1, 2, 3, 4, 5])
+        expected_y_mean = np.array([2, 32, 10.5, 17, 10])
+        expected_y_sigma = np.sqrt(np.array([2 / 3, 542, 1 / 4, 0, 0]))
+        self.assertTrue(np.allclose(expected_x_mean, x_mean))
+        self.assertTrue(np.allclose(expected_y_mean, y_mean))
+        self.assertTrue(np.allclose(expected_y_sigma, y_sigma))
+
+        sigma = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14])
+        x_mean, y_mean, y_sigma = mean_xy_data(x, y, sigma, method="iwv")
+        expected_y_mean = np.array([1.34693878, 23.31590234, 10.44137931, 17.0, 10.0])
+        expected_y_sigma = np.array([0.85714286, 2.57610543, 5.97927455, 10.0, 6.17470935])
+        self.assertTrue(np.allclose(expected_x_mean, x_mean))
+        self.assertTrue(np.allclose(expected_y_mean, y_mean))
+        self.assertTrue(np.allclose(expected_y_sigma, y_sigma))
+
+        x = np.array([1, 1, 1, 1, 2, 2, 2, 2])
+        y = np.array([2, 6, 100, 200, 17, 50, 60, 70])
+        series = np.array([0, 0, 1, 1, 0, 1, 1, 1])
+        series, x_mean, y_mean, y_sigma = multi_mean_xy_data(series, x, y, method="sample")
+        expected_x_mean = np.array([1, 2, 1, 2])
+        expected_y_mean = np.array([4, 17, 150, 60])
+        expected_y_sigma = np.sqrt(np.array([4.0, 0.0, 2500.0, 66.66666667]))
+        expected_series = np.array([0, 0, 1, 1])
+        self.assertTrue(np.allclose(expected_x_mean, x_mean))
+        self.assertTrue(np.allclose(expected_y_mean, y_mean))
+        self.assertTrue(np.allclose(expected_y_sigma, y_sigma))
+        self.assertTrue(np.allclose(expected_series, series))