RB experiment

qiskit_experiments/__init__.py

@@ -19,3 +19,5 @@
 
 # Experiment modules
 from . import composite
+from . import analysis
+from . import randomized_benchmarking

qiskit_experiments/randomized_benchmarking/__init__.py

@@ -0,0 +1,16 @@
+# 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.
+
+"""Randomized Benchmarking Experiment Classes."""
+
+from .rb_experiment import RBExperiment
+from .rb_analysis import RBAnalysis

qiskit_experiments/randomized_benchmarking/rb_analysis.py

@@ -0,0 +1,152 @@
+# 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.
+"""
+Standard RB analysis class.
+"""
+
+from typing import Optional, List
+
+import numpy as np
+from qiskit_experiments.base_analysis import BaseAnalysis
+from qiskit_experiments.analysis.curve_fitting import curve_fit
+from qiskit_experiments.analysis.data_processing import (
+    level2_probability,
+    mean_xy_data,
+    filter_data,
+)
+from qiskit_experiments.analysis.plotting import plot_curve_fit, plot_scatter, plot_errorbar
+
+try:
+    from matplotlib import pyplot as plt
+
+    HAS_MATPLOTLIB = True
+except ImportError:
+    HAS_MATPLOTLIB = False
+
+
+class RBAnalysis(BaseAnalysis):
+    """RB Analysis class."""
+
+    # pylint: disable = arguments-differ, invalid-name, attribute-defined-outside-init
+    def _run_analysis(
+        self,
+        experiment_data,
+        p0: Optional[List[float]] = None,
+        plot: bool = True,
+        ax: Optional["AxesSubplot"] = None,
+    ):
+        """Run analysis on circuit data.
+        Args:
+            experiment_data (ExperimentData): the experiment data to analyze.
+            p0: Optional, initial parameter values for curve_fit.
+            plot: If True generate a plot of fitted data.
+            ax: Optional, matplotlib axis to add plot to.
+        Returns:
+            tuple: A pair ``(analysis_result, figures)`` where
+                   ``analysis_results`` may be a single or list of
+                   AnalysisResult objects, and ``figures`` may be
+                   None, a single figure, or a list of figures.
+        """
+        self._num_qubits = len(experiment_data.data[0]["metadata"]["qubits"])
+        xdata, ydata, ydata_sigma = self._extract_data(experiment_data)
+
+        def fit_fun(x, a, alpha, b):
+            return a * alpha ** x + b
+
+        p0 = self._p0(xdata, ydata)
+        analysis_result = curve_fit(
+            fit_fun, xdata, ydata, p0, ydata_sigma, bounds=([0, 0, 0], [1, 1, 1])
+        )
+
+        # Add EPC data
+        popt = analysis_result["popt"]
+        popt_err = analysis_result["popt_err"]
+        scale = (2 ** self._num_qubits - 1) / (2 ** self._num_qubits)
+        analysis_result["EPC"] = scale * (1 - popt[1])
+        analysis_result["EPC_err"] = scale * popt_err[1] / popt[1]
+        analysis_result["plabels"] = ["A", "alpha", "B"]
+
+        if plot:
+            ax = plot_curve_fit(fit_fun, analysis_result, ax=ax)
+            ax = plot_scatter(xdata, ydata, ax=ax)
+            ax = plot_errorbar(xdata, ydata, ydata_sigma, ax=ax)
+            self._format_plot(ax, analysis_result)
+            analysis_result.plt = plt
+        return analysis_result, None
+
+    def _p0(self, xdata, ydata):
+        """Initial guess for the fitting function"""
+        fit_guess = [0.95, 0.99, 1 / 2 ** self._num_qubits]
+        # Use the first two points to guess the decay param
+        dcliff = xdata[1] - xdata[0]
+        dy = (ydata[1] - fit_guess[2]) / (ydata[0] - fit_guess[2])
+        alpha_guess = dy ** (1 / dcliff)
+        if alpha_guess < 1.0:
+            fit_guess[1] = alpha_guess
+
+        if ydata[0] > fit_guess[2]:
+            fit_guess[0] = (ydata[0] - fit_guess[2]) / fit_guess[1] ** xdata[0]
+
+        return fit_guess
+
+    def _extract_data(self, experiment_data, **filters):
+        """Extract the base data for the fitter from the experiment data.
+        Args:
+            data: the experiment data to analyze
+        Returns:
+            tuple: ``(xdata, ydata, ydata_sigma)`` , where
+               ``xdata`` is an array of unique x-values, ``ydata`` is an array of
+               sample mean y-values, and ``ydata_sigma`` is an array of sample standard
+               deviation of y values.
+        """
+        data = filter_data(experiment_data.data, **filters)
+        size = len(data)
+        xdata = np.zeros(size, dtype=int)
+        ydata = np.zeros(size, dtype=float)
+        ydata_var = np.zeros(size, dtype=float)
+        for i, datum in enumerate(data):
+            metadata = datum["metadata"]
+            xdata[i] = metadata["xdata"]
+            ydata[i], ydata_var[i] = level2_probability(datum, metadata["ylabel"])
+
+        ydata_sigma = np.sqrt(ydata_var)
+        xdata, ydata, ydata_sigma = mean_xy_data(xdata, ydata, ydata_sigma)
+        return (xdata, ydata, ydata_sigma)
+
+    @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_err = analysis_result["popt_err"][1]
+            epc = analysis_result["EPC"]
+            epc_err = analysis_result["EPC_err"]
+            box_text = "\u03B1:{:.4f} \u00B1 {:.4f}".format(alpha, alpha_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/rb_experiment.py

@@ -0,0 +1,149 @@
+# 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.
+"""
+Standard RB Experiment class.
+"""
+from typing import Union, Iterable, Optional
+
+import numpy as np
+from numpy.random import Generator, default_rng
+
+from qiskit import QuantumCircuit
+from qiskit.quantum_info import Clifford, random_clifford
+
+from qiskit_experiments.base_experiment import BaseExperiment
+from .rb_analysis import RBAnalysis
+
+
+class RBExperiment(BaseExperiment):
+    """RB Experiment class"""
+
+    # Analysis class for experiment
+    __analysis_class__ = RBAnalysis
+
+    def __init__(
+        self,
+        qubits: Union[int, Iterable[int]],
+        lengths: Iterable[int],
+        num_samples: int = 1,
+        seed: Optional[Union[int, Generator]] = None,
+        full_sampling: bool = False,
+    ):
+        """Standard randomized benchmarking experiment
+        Args:
+            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.
+        """
+        if not isinstance(seed, Generator):
+            self._rng = default_rng(seed=seed)
+        else:
+            self._rng = seed
+        self._lengths = list(lengths)
+        self._num_samples = num_samples
+        self._full_sampling = full_sampling
+        super().__init__(qubits)
+
+    # pylint: disable = arguments-differ
+    def circuits(self, backend=None):
+        """Return a list of RB circuits.
+        Args:
+            backend (Backend): Optional, a backend object.
+        Returns:
+            List[QuantumCircuit]: A list of :class:`QuantumCircuit`s.
+        """
+        circuits = []
+        for _ in range(self._num_samples):
+            circuits += self._sample_circuits(self._lengths, seed=self._rng)
+        return circuits
+
+    def transpiled_circuits(self, backend=None, **kwargs):
+        """Return a list of transpiled RB circuits.
+        Args:
+            backend (Backend): Optional, a backend object to use as the
+                               argument for the :func:`qiskit.transpile`
+                               function.
+            kwargs: kwarg options for the :func:`qiskit.transpile` function.
+        Returns:
+            List[QuantumCircuit]: A list of :class:`QuantumCircuit`s.
+        Raises:
+            QiskitError: if an initial layout is specified in the
+                         kwarg options for transpilation. The initial
+                         layout must be generated from the experiment.
+        """
+        circuits = super().transpiled_circuits(backend=backend, **kwargs)
+        return circuits
+
+    def _sample_circuits(
+        self, lengths: Iterable[int], seed: Optional[Union[int, Generator]] = None
+    ):
+        """Return a list RB circuits for the given lengths.
+        Args:
+            lengths: A list of RB sequences lengths.
+            seed: Seed or generator object for random number
+                  generation. If None default_rng will be used.
+        Returns:
+            List[QuantumCircuit]: A list of :class:`QuantumCircuit`s.
+        """
+        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
+            circuits += self._generate_circuit(elements, element_lengths)
+        return circuits
+
+    def _generate_circuit(self, elements: Iterable[Clifford], lengths: Iterable[int]):
+        """Return the RB circuits constructed from the given element list.
+        Args:
+            elements: A list of Clifford elements
+            lengths: A list of RB sequences lengths.
+        Returns:
+            List[QuantumCircuit]: A list of :class:`QuantumCircuit`s.
+        Additional information:
+            The circuits are constructed iteratively; each circuit is obtained
+            by extending the previous circuit (without the inversion and measurement gates)
+        """
+        qubits = list(range(self.num_qubits))
+        circuits = []
+
+        circ = QuantumCircuit(self.num_qubits)
+        circ.barrier(qubits)
+        circ_op = Clifford(np.eye(2 * self.num_qubits))
+
+        for current_length, group_elt in enumerate(elements):
+            circ_op = circ_op.compose(group_elt)
+            circ.append(group_elt, qubits)
+            circ.barrier(qubits)
+            if current_length + 1 in lengths:
+                # copy circuit and add inverse
+                inv = circ_op.adjoint()
+                rb_circ = circ.copy()
+                rb_circ.append(inv, qubits)
+                rb_circ.barrier(qubits)
+                rb_circ.metadata = {
+                    "experiment_type": self._type,
+                    "xdata": current_length + 1,
+                    "ylabel": self.num_qubits * "0",
+                    "group": "Clifford",
+                    "qubits": self.physical_qubits,
+                }
+                rb_circ.measure_all()
+                circuits.append(rb_circ)
+        return circuits