diff --git a/docs/tutorials/randomized_benchmarking.rst b/docs/tutorials/randomized_benchmarking.rst index d1c2aefd37..29640a8066 100644 --- a/docs/tutorials/randomized_benchmarking.rst +++ b/docs/tutorials/randomized_benchmarking.rst @@ -4,7 +4,7 @@ Randomized Benchmarking A randomized benchmarking (RB) experiment consists of the generation of random Clifford circuits on the given qubits such that the unitary computed by the circuits is the identity. After running the circuits, -the number of shots resulting in an error (i.e. an output different than +the number of shots resulting in an error (i.e. an output different than the ground state) are counted, and from this data one can infer error estimates for the quantum device, by calculating the Error Per Clifford. See `Qiskit @@ -15,8 +15,7 @@ for an explanation on the RB method, which is based on Ref. [1, 2]. import numpy as np from qiskit_experiments.library import StandardRB, InterleavedRB - from qiskit_experiments.framework import ParallelExperiment - from qiskit_experiments.library.randomized_benchmarking import RBUtils + from qiskit_experiments.framework import ParallelExperiment, BatchExperiment import qiskit.circuit.library as circuits # For simulation @@ -47,7 +46,7 @@ in order to generate the RB circuits and run them on a backend: sequences are constructed by appending additional Clifford samples to shorter sequences. The default is ``False`` -The analysis results of the RB Experiment includes: +The analysis results of the RB Experiment may include: - ``EPC``: The estimated Error Per Clifford @@ -60,6 +59,40 @@ The analysis results of the RB Experiment includes: Running a 1-qubit RB experiment ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +Standard RB experiment will provide you gate errors for every basis gates +constituting averaged Clifford gate. Note that you can only obtain a single EPC value :math:`\cal E` +from a single RB experiment. As such, computing the error values for multiple gates :math:`\{g_i\}` +requires some assumption of contribution of each gate to the total depolarizing error. +This is so called ``gate_error_ratio`` option you can find in analysis options. + +Provided that we have :math:`n_i` gates with independent error :math:`e_i` per Clifford, +the total EPC is estimated by the composition of error from every basis gate, + +.. math:: + + {\cal E} = 1 - \prod_{i} (1 - e_i)^{n_i} \sim \sum_{i} n_i e_i + O(e^2), + +where :math:`e_i \ll 1` and the higher order terms can be ignored. + +We cannot distinguish :math:`e_i` with a single EPC value :math:`\cal E` as explained, +however by defining an error ratio :math:`r_i` with respect to +some standard value :math:`e_0`, we can compute EPG :math:`e_i` for each basis gate. + +.. math:: + + {\cal E} \sim e_0 \sum_{i} n_i r_i + +The EPG of :math:`i` th basis gate will be + +.. math:: + + e_i \sim r_i e_0 = \dfrac{r_i{\cal E}}{\sum_{i} n_i r_i}. + +Because EPGs are computed based on this simple assumption, +this is not necessary representing the true gate error on the hardware. +If you have multiple kinds of basis gates with unclear error ratio :math:`r_i`, +interleaved RB experiment will always give you accurate error value :math:`e_i`. + .. jupyter-execute:: lengths = np.arange(1, 800, 200) @@ -73,6 +106,7 @@ Running a 1-qubit RB experiment results1 = expdata1.analysis_results() # View result data + print("Gate error ratio: %s" % expdata1.experiment.analysis.options.gate_error_ratio) display(expdata1.figure(0)) for result in results1: print(result) @@ -82,55 +116,92 @@ Running a 1-qubit RB experiment Running a 2-qubit RB experiment ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -Running a 1-qubit RB experiment and a 2-qubit RB experiment, in order to -calculate the gate error (EPG) of the ``cx`` gate: +In the same way we can compute EPC for two-qubit RB experiment. +However, the EPC value obtained by the experiment indicates a depolarization +which is a composition of underlying error channels for 2Q gates and 1Q gates in each qubit. +Usually 1Q gate contribution is small enough to ignore, but in case this +contribution is significant comparing to the 2Q gate error, +we can decompose the contribution of 1Q gates [3]. + +.. math:: + + \alpha_{2Q,C} = \frac{1}{5} \left( \alpha_0^{N_1/2} + \alpha_1^{N_1/2} + + 3 \alpha_0^{N_1/2} \alpha_1^{N_1/2} \right) \alpha_{01}^{N_2}, + +where :math:`\alpha_i` is the single qubit depolarizing parameter of channel :math:`i`, +and :math:`\alpha_{01}` is the two qubit depolarizing parameter of interest. +:math:`N_1` and :math:`N_2` are total count of single and two qubit gates, respectively. + +Note that the single qubit gate sequence in the channel :math:`i` may consist of +multiple kinds of basis gates :math:`\{g_{ij}\}` with different EPG :math:`e_{ij}`. +Therefore the :math:`\alpha_i^{N_1/2}` should be computed from EPGs, +rather than directly using the :math:`\alpha_i`, which is usually a composition of +depolarizing maps of every single qubit gate. +As such, EPGs should be measured in the separate single-qubit RBs in advance. + +.. math:: + + \alpha_i^{N_1/2} = \alpha_{i0}^{n_{i0}} \cdot \alpha_{i1}^{n_{i1}} \cdot ..., + +where :math:`\alpha_{ij}^{n_{ij}}` indicates a depolarization due to +a particular basis gate :math:`j` in the channel :math:`i`. +Here we assume EPG :math:`e_{ij}` corresponds to the depolarizing probability +of the map of :math:`g_{ij}`, and thus we can express :math:`\alpha_{ij}` with EPG. + +.. math:: + + e_{ij} = \frac{2^n - 1}{2^n} (1 - \alpha_{ij}) = \frac{1 - \alpha_{ij}}{2}, + +for the single qubit channel :math:`n=1`. Accordingly, + +.. math:: + + \alpha_i^{N_1/2} = \prod_{j} (1 - 2 e_{ij})^{n_{ij}}, + +as a composition of depolarization from every primitive gates per qubit. +This correction will give you two EPC values as a result of the two-qubit RB experiment. +The corrected EPC must be closer to the outcome of of interleaved RB. +The EPGs of two-qubit RB are analyzed with the corrected EPC if available. .. jupyter-execute:: - lengths = np.arange(1, 200, 30) + lengths_2_qubit = np.arange(1, 200, 30) + lengths_1_qubit = np.arange(1, 800, 200) num_samples = 10 seed = 1010 - qubits = (1,4) - + qubits = (1, 4) + # Run a 1-qubit RB expriment on qubits 1, 4 to determine the error-per-gate of 1-qubit gates - expdata_1q = {} - epg_1q = [] - lengths_1_qubit = np.arange(1, 800, 200) - for qubit in qubits: - exp = StandardRB([qubit], lengths_1_qubit, num_samples=num_samples, seed=seed) - expdata = exp.run(backend).block_for_results() - expdata_1q[qubit] = expdata - epg_1q += expdata.analysis_results() + single_exps = BatchExperiment( + [ + StandardRB([qubit], lengths_1_qubit, num_samples=num_samples, seed=seed) + for qubit in qubits + ], + flatten_results=True, + ) + expdata_1q = single_exps.run(backend).block_for_results() + .. jupyter-execute:: # Run an RB experiment on qubits 1, 4 - exp2 = StandardRB(qubits, lengths, num_samples=num_samples, seed=seed) + exp_2q = StandardRB(qubits, lengths_2_qubit, num_samples=num_samples, seed=seed) # Use the EPG data of the 1-qubit runs to ensure correct 2-qubit EPG computation - exp2.analysis.set_options(epg_1_qubit=epg_1q) + exp_2q.analysis.set_options(epg_1_qubit=expdata_1q.analysis_results()) # Run the 2-qubit experiment - expdata2 = exp2.run(backend).block_for_results() - - # View result data - results2 = expdata2.analysis_results() - -.. jupyter-execute:: + expdata_2q = exp_2q.run(backend).block_for_results() # View result data - display(expdata2.figure(0)) - for result in results2: + print("Gate error ratio: %s" % expdata_2q.experiment.analysis.options.gate_error_ratio) + display(expdata_2q.figure(0)) + for result in expdata_2q.analysis_results(): print(result) -.. jupyter-execute:: - # Compare the computed EPG of the cx gate with the backend's recorded cx gate error: - expected_epg = RBUtils.get_error_dict_from_backend(backend, qubits)[(qubits, 'cx')] - exp2_epg = expdata2.analysis_results("EPG_cx").value - - print("Backend's reported EPG of the cx gate:", expected_epg) - print("Experiment computed EPG of the cx gate:", exp2_epg) +Note that ``EPC_corrected`` value is smaller than one of raw ``EPC``, which indicates +contribution of depolarization from single-qubit error channels. Displaying the RB circuits @@ -253,12 +324,6 @@ different qubits (see Ref. [5]) par_exp = ParallelExperiment(exps) par_expdata = par_exp.run(backend).block_for_results() par_results = par_expdata.analysis_results() - - # View result data - for result in par_results: - print(result) - print("\nextra:") - print(result.extra) Viewing sub experiment data diff --git a/qiskit_experiments/curve_analysis/__init__.py b/qiskit_experiments/curve_analysis/__init__.py index 90a2c74555..efcf552de8 100644 --- a/qiskit_experiments/curve_analysis/__init__.py +++ b/qiskit_experiments/curve_analysis/__init__.py @@ -90,6 +90,7 @@ guess.constant_sinusoidal_offset guess.constant_spectral_offset guess.exp_decay + guess.rb_decay guess.full_width_half_max guess.frequency guess.max_height diff --git a/qiskit_experiments/curve_analysis/curve_analysis.py b/qiskit_experiments/curve_analysis/curve_analysis.py index 4100119922..581ab95e05 100644 --- a/qiskit_experiments/curve_analysis/curve_analysis.py +++ b/qiskit_experiments/curve_analysis/curve_analysis.py @@ -21,13 +21,12 @@ import inspect import warnings from abc import ABC -from typing import Any, Dict, List, Tuple, Callable, Union, Optional +from typing import Dict, List, Tuple, Callable, Union, Optional import numpy as np import uncertainties from uncertainties import unumpy as unp -from qiskit.providers import Backend from qiskit.utils import detach_prefix from qiskit_experiments.curve_analysis.curve_data import ( CurveData, @@ -253,14 +252,11 @@ def __init__(self): p: self.options.get(p, None) for p in self.__fixed_parameters__ } - #: Dict[str, Any]: Experiment metadata - self.__experiment_metadata = None - #: List[CurveData]: Processed experiment data set. self.__processed_data_set = list() - #: Backend: backend object used for experimentation - self.__backend = None + #: List[int]: Index of physical qubits + self._physical_qubits = None @classmethod def _fit_params(cls) -> List[str]: @@ -665,92 +661,6 @@ def _is_target_series(datum, **filters): raise AnalysisError(f"Not expected data label {formatted_data.label} != fit_ready.") self.__processed_data_set.append(formatted_data) - @property - def _experiment_type(self) -> str: - """Return type of experiment.""" - try: - return self.__experiment_metadata["experiment_type"] - except (TypeError, KeyError): - # Ignore experiment metadata is not set or key is not found - return None - - @property - def _num_qubits(self) -> int: - """Getter for qubit number.""" - try: - return len(self.__experiment_metadata["physical_qubits"]) - except (TypeError, KeyError): - # Ignore experiment metadata is not set or key is not found - return None - - @property - def _physical_qubits(self) -> List[int]: - """Getter for physical qubit indices.""" - try: - return list(self.__experiment_metadata["physical_qubits"]) - except (TypeError, KeyError): - # Ignore experiment metadata is not set or key is not found - return None - - @property - def _backend(self) -> Backend: - """Getter for backend object.""" - return self.__backend - - def _experiment_options(self, index: int = -1) -> Dict[str, Any]: - """Return the experiment options of given job index. - - Args: - index: Index of job metadata to extract. Default to -1 (latest). - - Returns: - Experiment options. This option is used for circuit generation. - """ - try: - return self.__experiment_metadata["job_metadata"][index]["experiment_options"] - except (TypeError, KeyError, IndexError): - # Ignore experiment metadata or job metadata is not set or key is not found - return None - - def _run_options(self, index: int = -1) -> Dict[str, Any]: - """Returns the run options of given job index. - - Args: - index: Index of job metadata to extract. Default to -1 (latest). - - Returns: - Run options. This option is used for backend execution. - """ - try: - return self.__experiment_metadata["job_metadata"][index]["run_options"] - except (TypeError, KeyError, IndexError): - # Ignore experiment metadata or job metadata is not set or key is not found - return None - - def _transpile_options(self, index: int = -1) -> Dict[str, Any]: - """Returns the transpile options of given job index. - - Args: - index: Index of job metadata to extract. Default to -1 (latest). - - Returns: - Transpile options. This option is used for circuit optimization. - """ - try: - return self.__experiment_metadata["job_metadata"][index]["transpile_options"] - except (TypeError, KeyError, IndexError): - # Ignore experiment metadata or job metadata is not set or key is not found - return None - - def _extra_metadata(self) -> Dict[str, Any]: - """Returns extra metadata. - - Returns: - Extra metadata explicitly added by the experiment subclass. - """ - exclude = ["experiment_type", "num_qubits", "physical_qubits", "job_metadata"] - return {k: v for k, v in self.__experiment_metadata.items() if k not in exclude} - def _data( self, series_name: Optional[str] = None, @@ -793,6 +703,10 @@ def _data( raise AnalysisError(f"Specified series {series_name} is not defined in this analysis.") + @property + def _num_qubits(self) -> int: + return len(self._physical_qubits) + def _run_analysis( self, experiment_data: ExperimentData ) -> Tuple[List[AnalysisResultData], List["pyplot.Figure"]]: @@ -824,15 +738,8 @@ def _run_analysis( # get experiment metadata try: - self.__experiment_metadata = experiment_data.metadata - - except AttributeError: - pass - - # get backend - try: - self.__backend = experiment_data.backend - except AttributeError: + self._physical_qubits = experiment_data.metadata["physical_qubits"] + except KeyError: pass # @@ -850,11 +757,42 @@ def _run_analysis( # Qiskit DataProcessor instance. May need calibration. data_processor.train(data=experiment_data.data()) + # Initialize fit figure canvas + if self.options.plot: + self.drawer.initialize_canvas() + # # 3. Extract curve entries from experiment data # self._extract_curves(experiment_data=experiment_data, data_processor=data_processor) + # TODO remove _data method dependency in follow-up + # self.__processed_data_set will be removed from instance. + + # Draw raw data + if self.options.plot and self.options.plot_raw_data: + for s in self.__series__: + raw_data = self._data(label="raw_data", series_name=s.name) + self.drawer.draw_raw_data( + x_data=raw_data.x, + y_data=raw_data.y, + ax_index=s.canvas, + ) + + # Draw formatted data + if self.options.plot: + for s in self.__series__: + curve_data = self._data(label="fit_ready", series_name=s.name) + self.drawer.draw_formatted_data( + x_data=curve_data.x, + y_data=curve_data.y, + y_err_data=curve_data.y_err, + name=s.name, + ax_index=s.canvas, + color=s.plot_color, + marker=s.plot_symbol, + ) + # # 4. Run fitting # @@ -985,75 +923,51 @@ def _run_analysis( ) analysis_results.append(raw_data_entry) - # - # 6. Create figures - # - if self.options.plot: - # Initialize axis - self.drawer.initialize_canvas() - # Write raw data - if self.options.plot_raw_data: - for s in self.__series__: - raw_data = self._data(label="raw_data", series_name=s.name) - self.drawer.draw_raw_data( - x_data=raw_data.x, - y_data=raw_data.y, - ax_index=s.canvas, - ) - # Write data points + # Draw fit results if fitting succeeded + if self.options.plot and fit_result: for s in self.__series__: - curve_data = self._data(label="fit_ready", series_name=s.name) - self.drawer.draw_formatted_data( - x_data=curve_data.x, - y_data=curve_data.y, - y_err_data=curve_data.y_err, - name=s.name, + interp_x = np.linspace(*fit_result.x_range, 100) + + params = {} + for fitpar in s.signature: + if fitpar in self.options.fixed_parameters: + params[fitpar] = self.options.fixed_parameters[fitpar] + else: + params[fitpar] = fit_result.fitval(fitpar) + + y_data_with_uncertainty = s.fit_func(interp_x, **params) + y_mean = unp.nominal_values(y_data_with_uncertainty) + y_std = unp.std_devs(y_data_with_uncertainty) + # Draw fit line + self.drawer.draw_fit_line( + x_data=interp_x, + y_data=y_mean, ax_index=s.canvas, color=s.plot_color, - marker=s.plot_symbol, ) - # Write fit results if fitting succeeded - if fit_result: - for s in self.__series__: - interp_x = np.linspace(*fit_result.x_range, 100) - - params = {} - for fitpar in s.signature: - if fitpar in self.options.fixed_parameters: - params[fitpar] = self.options.fixed_parameters[fitpar] - else: - params[fitpar] = fit_result.fitval(fitpar) - - y_data_with_uncertainty = s.fit_func(interp_x, **params) - y_mean = unp.nominal_values(y_data_with_uncertainty) - y_std = unp.std_devs(y_data_with_uncertainty) - # Draw fit line - self.drawer.draw_fit_line( - x_data=interp_x, - y_data=y_mean, - ax_index=s.canvas, - color=s.plot_color, - ) - # Draw confidence intervals with different n_sigma - sigmas = unp.std_devs(y_data_with_uncertainty) - if np.isfinite(sigmas).all(): - for n_sigma, alpha in self.drawer.options.plot_sigma: - self.drawer.draw_confidence_interval( - x_data=interp_x, - y_ub=y_mean + n_sigma * y_std, - y_lb=y_mean - n_sigma * y_std, - ax_index=s.canvas, - alpha=alpha, - color=s.plot_color, - ) - - # Write fitting report - report_description = "" - for res in analysis_results: - if isinstance(res.value, (float, uncertainties.UFloat)): - report_description += f"{analysis_result_to_repr(res)}\n" - report_description += r"Fit $\chi^2$ = " + f"{fit_result.reduced_chisq: .4g}" - self.drawer.draw_fit_report(description=report_description) + # Draw confidence intervals with different n_sigma + sigmas = unp.std_devs(y_data_with_uncertainty) + if np.isfinite(sigmas).all(): + for n_sigma, alpha in self.drawer.options.plot_sigma: + self.drawer.draw_confidence_interval( + x_data=interp_x, + y_ub=y_mean + n_sigma * y_std, + y_lb=y_mean - n_sigma * y_std, + ax_index=s.canvas, + alpha=alpha, + color=s.plot_color, + ) + + # Draw fitting report + report_description = "" + for res in analysis_results: + if isinstance(res.value, (float, uncertainties.UFloat)): + report_description += f"{analysis_result_to_repr(res)}\n" + report_description += r"Fit $\chi^2$ = " + f"{fit_result.reduced_chisq: .4g}" + self.drawer.draw_fit_report(description=report_description) + + # Output figure + if self.options.plot: self.drawer.format_canvas() figures = [self.drawer.figure] else: diff --git a/qiskit_experiments/curve_analysis/guess.py b/qiskit_experiments/curve_analysis/guess.py index 6abeae5c78..248ffc2116 100644 --- a/qiskit_experiments/curve_analysis/guess.py +++ b/qiskit_experiments/curve_analysis/guess.py @@ -345,3 +345,67 @@ def constant_sinusoidal_offset(y: np.ndarray) -> float: minv, _ = min_height(y, percentile=5) return 0.5 * (maxv + minv) + + +def rb_decay( + x: np.ndarray, + y: np.ndarray, + b: float = 0.5, + a: Optional[float] = None, +) -> float: + r"""Get base of exponential decay function which is assumed to be close to 1. + + This assumes following model: + + .. math:: + + y(x) = a \alpha^x + b. + + This model can be often seen in randomized benchmarking analysis. + Now we assume :math:`\alpha \simeq 1` and replace it with :math:`\alpha = 1 - p` + where :math:`p \ll 1`. Then, above function can be approximated to + + .. math:: + + y(x) = a \left(1 - px + \frac{x(x-1)}{2} p^2 + ...\right) + b, + + here we assume :math:`px \ll 1` to ignore higher order terms of :math:`p`. + Likely parameter :math:`p` of interest can be solved as + + .. math:: + + p(x) = \frac{a + b - y(x)}{ax}. + + The :math:`p(x)` is averaged with the weight of :math:`1/x`. + When estimated :math:`\alpha = 1-p < 0.9`, this calcualtes following value instead. + + .. math:: + + \alpha \simeq dy ^ {1 / dx}, + + where :math:`dx = x[1] - x[0]` and :math:`dy = (y[1] - b) / (y[0] - b)`. + + Args: + x: Array of x values. + y: Array of y values. + a: Coefficient of decay function. If not provided this is defaults to :math:`1-b`. + :math:`y(0) = 1.0`. + b: Asymptote of decay function. + + Returns: + Base of decay function. + """ + if a is None: + a = 1 - b + + a_guess_tmp = np.average(1 - (a + b - y) / (a * x), weights=1 / x) + + if a_guess_tmp > 0.9: + return a_guess_tmp + + # Violate assumption, then switch to the conventional method + # Use the first two points to guess the decay param + dx = x[1] - x[0] + dy = (y[1] - b) / (y[0] - b) + + return dy ** (1 / dx) diff --git a/qiskit_experiments/framework/base_experiment.py b/qiskit_experiments/framework/base_experiment.py index 8f08fd3481..e0656eb986 100644 --- a/qiskit_experiments/framework/base_experiment.py +++ b/qiskit_experiments/framework/base_experiment.py @@ -261,12 +261,12 @@ def run( # Finalize experiment before executions experiment._finalize() - # Initialize result container - experiment_data = experiment._initialize_experiment_data() - # Generate and transpile circuits transpiled_circuits = experiment._transpiled_circuits() + # Initialize result container + experiment_data = experiment._initialize_experiment_data() + # Run options run_opts = experiment.run_options.__dict__ diff --git a/qiskit_experiments/library/randomized_benchmarking/interleaved_rb_analysis.py b/qiskit_experiments/library/randomized_benchmarking/interleaved_rb_analysis.py index 8e9d91e016..e7e5fd7555 100644 --- a/qiskit_experiments/library/randomized_benchmarking/interleaved_rb_analysis.py +++ b/qiskit_experiments/library/randomized_benchmarking/interleaved_rb_analysis.py @@ -18,10 +18,8 @@ import qiskit_experiments.curve_analysis as curve from qiskit_experiments.framework import AnalysisResultData -from .rb_analysis import RBAnalysis - -class InterleavedRBAnalysis(RBAnalysis): +class InterleavedRBAnalysis(curve.CurveAnalysis): r"""A class to analyze interleaved randomized benchmarking experiment. # section: overview @@ -69,22 +67,20 @@ class InterleavedRBAnalysis(RBAnalysis): # section: fit_parameters defpar a: desc: Height of decay curve. - init_guess: Determined by the average :math:`a` of the standard and interleaved RB. + init_guess: Determined by :math:`1 - b`. bounds: [0, 1] defpar b: desc: Base line. - init_guess: Determined by the average :math:`b` of the standard and interleaved RB. - Usually equivalent to :math:`(1/2)^n` where :math:`n` is number of qubit. + init_guess: Determined by :math:`(1/2)^n` where :math:`n` is number of qubit. bounds: [0, 1] defpar \alpha: desc: Depolarizing parameter. - init_guess: Determined by the slope of :math:`(y - b)^{-x}` of the first and the - second data point of the standard RB. + init_guess: Determined by :func:`~rb_decay` with standard RB curve. bounds: [0, 1] defpar \alpha_c: desc: Ratio of the depolarizing parameter of interleaved RB to standard RB curve. - init_guess: Estimate :math:`\alpha' = \alpha_c \alpha` from the - interleaved RB curve, then divide this by the initial guess of :math:`\alpha`. + init_guess: Determined by alpha of interleaved RB curve devided by one of + standard RB curve. Both alpha values are estimated by :func:`~rb_decay`. bounds: [0, 1] # section: reference @@ -140,27 +136,60 @@ def _generate_fit_guesses( b=(0, 1), ) + b_guess = 1 / 2**self._num_qubits + a_guess = 1 - b_guess + # for standard RB curve std_curve = self._data(series_name="Standard") - opt_std = user_opt.copy() - opt_std = self._initial_guess(opt_std, std_curve.x, std_curve.y, self._num_qubits) + alpha_std = curve.guess.rb_decay(std_curve.x, std_curve.y, a=a_guess, b=b_guess) # for interleaved RB curve int_curve = self._data(series_name="Interleaved") - opt_int = user_opt.copy() - if opt_int.p0["alpha_c"] is not None: - opt_int.p0["alpha"] = opt_std.p0["alpha"] * opt_int.p0["alpha_c"] - opt_int = self._initial_guess(opt_int, int_curve.x, int_curve.y, self._num_qubits) + alpha_int = curve.guess.rb_decay(int_curve.x, int_curve.y, a=a_guess, b=b_guess) + + alpha_c = min(alpha_int / alpha_std, 1.0) user_opt.p0.set_if_empty( - a=np.mean([opt_std.p0["a"], opt_int.p0["a"]]), - alpha=opt_std.p0["alpha"], - alpha_c=min(opt_int.p0["alpha"] / opt_std.p0["alpha"], 1), - b=np.mean([opt_std.p0["b"], opt_int.p0["b"]]), + b=b_guess, + a=a_guess, + alpha=alpha_std, + alpha_c=alpha_c, ) return user_opt + def _format_data(self, data: curve.CurveData) -> curve.CurveData: + """Data format with averaging with sampling strategy.""" + # TODO Eventually move this to data processor, then create RB data processor. + + # take average over the same x value by regenerating sigma from variance of y values + series, xdata, ydata, sigma, shots = curve.data_processing.multi_mean_xy_data( + series=data.data_index, + xdata=data.x, + ydata=data.y, + sigma=data.y_err, + shots=data.shots, + method="sample", + ) + + # sort by x value in ascending order + series, xdata, ydata, sigma, shots = curve.data_processing.data_sort( + series=series, + xdata=xdata, + ydata=ydata, + sigma=sigma, + shots=shots, + ) + + return curve.CurveData( + label="fit_ready", + x=xdata, + y=ydata, + y_err=sigma, + shots=shots, + data_index=series, + ) + def _extra_database_entry(self, fit_data: curve.FitData) -> List[AnalysisResultData]: """Calculate EPC.""" nrb = 2**self._num_qubits diff --git a/qiskit_experiments/library/randomized_benchmarking/rb_analysis.py b/qiskit_experiments/library/randomized_benchmarking/rb_analysis.py index d3253fbe65..43e0b66bdf 100644 --- a/qiskit_experiments/library/randomized_benchmarking/rb_analysis.py +++ b/qiskit_experiments/library/randomized_benchmarking/rb_analysis.py @@ -12,16 +12,22 @@ """ Standard RB analysis class. """ +import warnings +from collections import defaultdict +from typing import Dict, List, Sequence, Tuple, Union, Optional, TYPE_CHECKING -from typing import List, Union +from qiskit.exceptions import QiskitError -import numpy as np import qiskit_experiments.curve_analysis as curve -from qiskit_experiments.curve_analysis.data_processing import multi_mean_xy_data, data_sort -from qiskit_experiments.database_service.device_component import Qubit -from qiskit_experiments.framework import AnalysisResultData +from qiskit_experiments.exceptions import AnalysisError +from qiskit_experiments.framework import AnalysisResultData, ExperimentData +from qiskit_experiments.database_service import DbAnalysisResultV1 -from .rb_utils import RBUtils +if TYPE_CHECKING: + from uncertainties import UFloat + +# A dictionary key of qubit aware quantum instruciton; type alias for better readability +QubitGateTuple = Tuple[Tuple[int, ...], str] class RBAnalysis(curve.CurveAnalysis): @@ -32,6 +38,12 @@ class RBAnalysis(curve.CurveAnalysis): This series is fit by the exponential decay function. From the fit :math:`\alpha` value this analysis estimates the error per Clifford (EPC). + When analysis option ``gate_error_ratio`` is provided, this analysis also estimates + errors of individual gates assembling a Clifford gate. + In computation of two-qubit EPC, this analysis can also decompose + the contribution from the underlying single qubit depolarizing channels when + ``epg_1_qubit`` analysis option is provided [1]. + # section: fit_model .. math:: @@ -40,19 +52,20 @@ class RBAnalysis(curve.CurveAnalysis): # section: fit_parameters defpar a: desc: Height of decay curve. - init_guess: Determined by :math:`(y - b) / \alpha^x`. + init_guess: Determined by :math:`1 - b`. bounds: [0, 1] defpar b: desc: Base line. - init_guess: Determined by the average :math:`b` of the standard and interleaved RB. - Usually equivalent to :math:`(1/2)^n` where :math:`n` is number of qubit. + init_guess: Determined by :math:`(1/2)^n` where :math:`n` is number of qubit. bounds: [0, 1] defpar \alpha: desc: Depolarizing parameter. - init_guess: Determined by the slope of :math:`(y - b)^{-x}` of the first and the - second data point. + init_guess: Determined by :func:`~rb_decay`. bounds: [0, 1] + # section: reference + .. ref_arxiv:: 1 1712.06550 + """ __series__ = [ @@ -65,20 +78,23 @@ class RBAnalysis(curve.CurveAnalysis): ) ] + def __init__(self): + super().__init__() + self._gate_counts_per_clifford = None + @classmethod def _default_options(cls): """Default analysis options. Analysis Options: - error_dict (Dict[Tuple[Iterable[int], str], float]): Optional. - Error estimates for gates from the backend properties. - epg_1_qubit (Dict[int, Dict[str, float]]) : Optional. - EPG data for the 1-qubit gate involved, - assumed to have been obtained from previous experiments. - This is used to estimate the 2-qubit EPG. - gate_error_ratio (Dict[str, float]): An estimate for the ratios - between errors on different gates. - + gate_error_ratio (Optional[Dict[str, float]]): A dictionary with gate name keys + and error ratio values used when calculating EPG from the estimated EPC. + The default value will use standard gate error ratios. + If you don't know accurate error ratio between your basis gates, + you can skip analysis of EPGs by setting this options to ``None``. + epg_1_qubit (List[DbAnalysisResultV1]): Analysis results from previous RB experiments + for individual single qubit gates. If this is provided, EPC of + 2Q RB is corected to exclude the deporalization of underlying 1Q channels. """ default_options = super()._default_options() default_options.curve_plotter.set_options( @@ -87,12 +103,20 @@ def _default_options(cls): ) default_options.plot_raw_data = True default_options.result_parameters = ["alpha"] - default_options.error_dict = None + default_options.gate_error_ratio = "default" default_options.epg_1_qubit = None - default_options.gate_error_ratio = None return default_options + def set_options(self, **fields): + if "error_dict" in fields: + warnings.warn( + "Option 'error_dict' has been removed and merged into 'gate_error_ratio'.", + DeprecationWarning, + ) + fields["gate_error_ratio"] = fields.pop("error_dict") + super().set_options(**fields) + def _generate_fit_guesses( self, user_opt: curve.FitOptions ) -> Union[curve.FitOptions, List[curve.FitOptions]]: @@ -112,33 +136,24 @@ def _generate_fit_guesses( b=(0, 1), ) - return self._initial_guess(user_opt, curve_data.x, curve_data.y, self._num_qubits) + b_guess = 1 / 2**self._num_qubits + a_guess = 1 - b_guess + alpha_guess = curve.guess.rb_decay(curve_data.x, curve_data.y, a=a_guess, b=b_guess) - @staticmethod - def _initial_guess( - opt: curve.FitOptions, x_values: np.ndarray, y_values: np.ndarray, num_qubits: int - ) -> curve.FitOptions: - """Create initial guess with experiment data.""" - opt.p0.set_if_empty(b=1 / 2**num_qubits) - - # Use the first two points to guess the decay param - dcliff = x_values[1] - x_values[0] - dy = (y_values[1] - opt.p0["b"]) / (y_values[0] - opt.p0["b"]) - alpha_guess = dy ** (1 / dcliff) - - opt.p0.set_if_empty(alpha=alpha_guess if alpha_guess < 1.0 else 0.99) - - if y_values[0] > opt.p0["b"]: - opt.p0.set_if_empty(a=(y_values[0] - opt.p0["b"]) / (opt.p0["alpha"] ** x_values[0])) - else: - opt.p0.set_if_empty(a=0.95) + user_opt.p0.set_if_empty( + b=b_guess, + a=a_guess, + alpha=alpha_guess, + ) - return opt + return user_opt def _format_data(self, data: curve.CurveData) -> curve.CurveData: """Data format with averaging with sampling strategy.""" + # TODO Eventually move this to data processor, then create RB data processor. + # take average over the same x value by regenerating sigma from variance of y values - series, xdata, ydata, sigma, shots = multi_mean_xy_data( + series, xdata, ydata, sigma, shots = curve.data_processing.multi_mean_xy_data( series=data.data_index, xdata=data.x, ydata=data.y, @@ -148,7 +163,7 @@ def _format_data(self, data: curve.CurveData) -> curve.CurveData: ) # sort by x value in ascending order - series, xdata, ydata, sigma, shots = data_sort( + series, xdata, ydata, sigma, shots = curve.data_processing.data_sort( series=series, xdata=xdata, ydata=ydata, @@ -168,6 +183,8 @@ def _format_data(self, data: curve.CurveData) -> curve.CurveData: def _extra_database_entry(self, fit_data: curve.FitData) -> List[AnalysisResultData]: """Calculate EPC.""" extra_entries = [] + # pylint: disable=assignment-from-none + quality = self._evaluate_quality(fit_data) # Calculate EPC alpha = fit_data.fitval("alpha") @@ -179,60 +196,241 @@ def _extra_database_entry(self, fit_data: curve.FitData) -> List[AnalysisResultD name="EPC", value=epc, chisq=fit_data.reduced_chisq, - quality=self._evaluate_quality(fit_data), + quality=quality, ) ) - # Calculate EPG - if not self.options.gate_error_ratio: - # we attempt to get the ratio from the backend properties - if not self.options.error_dict: - gate_error_ratio = RBUtils.get_error_dict_from_backend( - backend=self._backend, qubits=self._physical_qubits - ) - else: - gate_error_ratio = self.options.error_dict - else: - gate_error_ratio = self.options.gate_error_ratio - - count_ops = [] - for meta in self._data(label="raw_data").metadata: - count_ops += meta.get("count_ops", []) - - if len(count_ops) > 0 and gate_error_ratio is not None: - gates_per_clifford = RBUtils.gates_per_clifford(count_ops) - num_qubits = len(self._physical_qubits) - - if num_qubits == 1: - epg_dict = RBUtils.calculate_1q_epg( - epc, - self._physical_qubits, - gate_error_ratio, - gates_per_clifford, - ) - elif num_qubits == 2: - epg_1_qubit = self.options.epg_1_qubit - epg_dict = RBUtils.calculate_2q_epg( - epc, - self._physical_qubits, - gate_error_ratio, - gates_per_clifford, - epg_1_qubit=epg_1_qubit, + # Correction for 1Q depolarizing channel if EPGs are provided + if self.options.epg_1_qubit and self._num_qubits == 2: + epc = _exclude_1q_error( + epc=epc, + qubits=self._physical_qubits, + gate_counts_per_clifford=self._gate_counts_per_clifford, + extra_analyses=self.options.epg_1_qubit, + ) + extra_entries.append( + AnalysisResultData( + name="EPC_corrected", + value=epc, + chisq=fit_data.reduced_chisq, + quality=quality, ) - else: - # EPG calculation is not supported for more than 3 qubits RB - epg_dict = None + ) + # Calculate EPG + if self._gate_counts_per_clifford is not None and self.options.gate_error_ratio: + epg_dict = _calculate_epg( + epc=epc, + qubits=self._physical_qubits, + gate_error_ratio=self.options.gate_error_ratio, + gate_counts_per_clifford=self._gate_counts_per_clifford, + ) if epg_dict: - for qubits, gate_dict in epg_dict.items(): - for gate, value in gate_dict.items(): - extra_entries.append( - AnalysisResultData( - f"EPG_{gate}", - value, - chisq=fit_data.reduced_chisq, - quality=self._evaluate_quality(fit_data), - device_components=[Qubit(i) for i in qubits], - ) + for gate, epg_val in epg_dict.items(): + extra_entries.append( + AnalysisResultData( + name=f"EPG_{gate}", + value=epg_val, + chisq=fit_data.reduced_chisq, + quality=quality, ) + ) + return extra_entries + + def _run_analysis( + self, experiment_data: ExperimentData + ) -> Tuple[List[AnalysisResultData], List["pyplot.Figure"]]: + + if self.options.gate_error_ratio is not None: + # If gate error ratio is not False, EPG analysis is enabled. + # Here analysis prepares gate error ratio and gate counts for EPC to EPG conversion. + + # If gate count dictionary is not set it will compute counts from circuit metadata. + avg_gpc = defaultdict(float) + n_circs = len(experiment_data.data()) + for circ_result in experiment_data.data(): + try: + count_ops = circ_result["metadata"]["count_ops"] + except KeyError as ex: + raise AnalysisError( + "'count_ops' key is not found in the circuit metadata. " + "This analysis cannot compute error per gates. " + "Please disable this with 'gate_error_ratio=False'." + ) from ex + nclif = circ_result["metadata"]["xval"] + for (qinds, gate), count in count_ops: + formatted_key = tuple(sorted(qinds)), gate + avg_gpc[formatted_key] += count / nclif / n_circs + self._gate_counts_per_clifford = dict(avg_gpc) + + if self.options.gate_error_ratio == "default": + # Gate error dict is computed for gates appearing in counts dictionary + # Error ratio among gates is determined based on the predefined lookup table. + # This is not always accurate for every quantum backends. + gate_error_ratio = {} + for qinds, gate in self._gate_counts_per_clifford.keys(): + if set(qinds) != set(experiment_data.metadata["physical_qubits"]): + continue + gate_error_ratio[gate] = _lookup_epg_ratio(gate, len(qinds)) + self.set_options(gate_error_ratio=gate_error_ratio) + + return super()._run_analysis(experiment_data) + + +def _lookup_epg_ratio(gate: str, n_qubits: int) -> Union[None, int]: + """A helper method to look-up preset gate error ratio for given basis gate name. + + In the table the error ratio is defined based on the count of + typical assembly gate in the gate decomposition. + For example, "u3" gate can be decomposed into two "sx" gates. + In this case, the ratio of "u3" gate error becomes 2. + + .. note:: + + This table is not aware of the actual waveform played on the hardware, + and the returned error ratio is just a guess. + To be precise, user can always set "gate_error_ratio" option of the experiment. + + Args: + gate: Name of the gate. + n_qubits: Number of qubits measured in the RB experiments. + + Returns: + Corresponding error ratio. + + Raises: + QiskitError: When number of qubit is more than three. + """ + + # Gate count in (X, SX)-based decomposition. VZ gate contribution is ignored. + # Amplitude or duration modulated pulse implementation is not considered. + standard_1q_ratio = { + "u1": 0.0, + "u2": 1.0, + "u3": 2.0, + "u": 2.0, + "p": 0.0, + "x": 1.0, + "y": 1.0, + "z": 0.0, + "t": 0.0, + "tdg": 0.0, + "s": 0.0, + "sdg": 0.0, + "sx": 1.0, + "sxdg": 1.0, + "rx": 2.0, + "ry": 2.0, + "rz": 0.0, + "id": 0.0, + "h": 1.0, + } + + # Gate count in (CX, CSX)-based decomposition, 1q gate contribution is ignored. + # Amplitude or duration modulated pulse implementation is not considered. + standard_2q_ratio = { + "swap": 3.0, + "rxx": 2.0, + "rzz": 2.0, + "cx": 1.0, + "cy": 1.0, + "cz": 1.0, + "ch": 1.0, + "crx": 2.0, + "cry": 2.0, + "crz": 2.0, + "csx": 1.0, + "cu1": 2.0, + "cp": 2.0, + "cu": 2.0, + "cu3": 2.0, + } + + if n_qubits == 1: + return standard_1q_ratio.get(gate, None) + + if n_qubits == 2: + return standard_2q_ratio.get(gate, None) + + raise QiskitError( + f"Standard gate error ratio for {n_qubits} qubit RB is not provided. " + "Please explicitly set 'gate_error_ratio' option of the experiment." + ) + + +def _calculate_epg( + epc: Union[float, "UFloat"], + qubits: Sequence[int], + gate_error_ratio: Dict[str, float], + gate_counts_per_clifford: Dict[QubitGateTuple, float], +) -> Dict[str, Union[float, "UFloat"]]: + """A helper mehtod to compute EPGs of basis gates from fit EPC value. + + Args: + epc: Error per Clifford. + qubits: List of qubits used in the experiment. + gate_error_ratio: A dictionary of assumed ratio of errors among basis gates. + gate_counts_per_clifford: Basis gate counts per Clifford gate. + + Returns: + A dictionary of gate errors keyed on the gate name. + """ + norm = 0 + for gate, r_epg in gate_error_ratio.items(): + formatted_key = tuple(sorted(qubits)), gate + norm += r_epg * gate_counts_per_clifford.get(formatted_key, 0.0) + + epgs = {} + for gate, r_epg in gate_error_ratio.items(): + epgs[gate] = r_epg * epc / norm + return epgs + + +def _exclude_1q_error( + epc: Union[float, "UFloat"], + qubits: Tuple[int, int], + gate_counts_per_clifford: Dict[QubitGateTuple, float], + extra_analyses: Optional[List[DbAnalysisResultV1]], +) -> Union[float, "UFloat"]: + """A helper method to exclude contribution of single qubit gates from 2Q EPC. + + Args: + epc: EPC from 2Q RB experiment. + qubits: Index of two qubits used for 2Q RB experiment. + gate_counts_per_clifford: Basis gate counts per 2Q Clifford gate. + extra_analyses: Analysis results containing depolarizing parameters of 1Q RB experiments. + + Returns: + Corrected 2Q EPC. + """ + # Extract EPC of non-measured qubits from previous experiments + epg_1qs = {} + for analyis_data in extra_analyses: + if ( + not analyis_data.name.startswith("EPG_") + or len(analyis_data.device_components) > 1 + or not str(analyis_data.device_components[0]).startswith("Q") + ): + continue + qind = analyis_data.device_components[0]._index + gate = analyis_data.name[4:] + formatted_key = (qind,), gate + epg_1qs[formatted_key] = analyis_data.value + + if not epg_1qs: + return epc + + # Convert 2Q EPC into depolarizing parameter alpha + alpha_c_2q = 1 - 4 / 3 * epc + + # Estimate composite alpha of 1Q channels + alpha_i = [1.0, 1.0] + for q_gate_tup, epg in epg_1qs.items(): + n_gate = gate_counts_per_clifford.get(q_gate_tup, 0.0) + aind = qubits.index(q_gate_tup[0][0]) + alpha_i[aind] *= (1 - 2 * epg) ** n_gate + alpha_c_1q = 1 / 5 * (alpha_i[0] + alpha_i[1] + 3 * alpha_i[0] * alpha_i[1]) + + # Corrected 2Q channel EPC + return 3 / 4 * (1 - (alpha_c_2q / alpha_c_1q)) diff --git a/qiskit_experiments/library/randomized_benchmarking/rb_experiment.py b/qiskit_experiments/library/randomized_benchmarking/rb_experiment.py index ca3db9749f..265f0dc23a 100644 --- a/qiskit_experiments/library/randomized_benchmarking/rb_experiment.py +++ b/qiskit_experiments/library/randomized_benchmarking/rb_experiment.py @@ -12,6 +12,8 @@ """ Standard RB Experiment class. """ +import logging +from collections import defaultdict from typing import Union, Iterable, Optional, List, Sequence import numpy as np @@ -20,14 +22,14 @@ from qiskit import QuantumCircuit, QiskitError from qiskit.quantum_info import Clifford -from qiskit.circuit import Gate from qiskit.providers.backend import Backend -from qiskit_experiments.framework import BaseExperiment, ParallelExperiment, Options +from qiskit_experiments.framework import BaseExperiment, Options from qiskit_experiments.framework.restless_mixin import RestlessMixin from .rb_analysis import RBAnalysis from .clifford_utils import CliffordUtils -from .rb_utils import RBUtils + +LOG = logging.getLogger(__name__) class StandardRB(BaseExperiment, RestlessMixin): @@ -46,17 +48,12 @@ class StandardRB(BaseExperiment, RestlessMixin): the ground state, fits an exponentially decaying curve, and estimates the Error Per Clifford (EPC), as described in Refs. [1, 2]. - See :class:`RBUtils` documentation for additional information - on estimating the Error Per Gate (EPG) for 1-qubit and 2-qubit gates, - from 1-qubit and 2-qubit standard RB experiments, by Ref. [3]. - # section: analysis_ref :py:class:`RBAnalysis` # section: reference .. ref_arxiv:: 1 1009.3639 .. ref_arxiv:: 2 1109.6887 - .. ref_arxiv:: 3 1712.06550 """ @@ -77,13 +74,11 @@ def __init__( backend: The backend to run the experiment on. num_samples: Number of samples to generate for each sequence length. seed: Optional, seed used to initialize ``numpy.random.default_rng``. - when generating circuits. The ``default_rng`` will be initialized - with this seed value everytime :meth:`circuits` is called. + when generating circuits. The ``default_rng`` will be initialized + with this seed value everytime :meth:`circuits` is called. 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. - The default is False. + all lengths. If False for sample of lengths longer sequences are constructed by + appending additional Clifford samples to shorter sequences. The default is ``False``. """ # Initialize base experiment super().__init__(qubits, analysis=RBAnalysis(), backend=backend) @@ -192,7 +187,7 @@ def _generate_circuit( group_elt_gate = group_elt_circ group_elt_op = Clifford(group_elt_circ) - if not isinstance(group_elt_gate, Gate): + if hasattr(group_elt_gate, "to_gate"): group_elt_gate = group_elt_gate.to_gate() circ_op = circ_op.compose(group_elt_op) for circ in circs: @@ -214,25 +209,31 @@ def _generate_circuit( circuits.append(rb_circ) return circuits - def _get_circuit_metadata(self, circuit): - if circuit.metadata["experiment_type"] == self._type: - return circuit.metadata - if circuit.metadata["experiment_type"] == ParallelExperiment.__name__: - for meta in circuit.metadata["composite_metadata"]: - if meta["physical_qubits"] == self.physical_qubits: - return meta - return None - def _transpiled_circuits(self) -> List[QuantumCircuit]: """Return a list of experiment circuits, transpiled.""" transpiled = super()._transpiled_circuits() - for c in transpiled: - meta = self._get_circuit_metadata(c) - if meta is not None: - c_count_ops = RBUtils.count_ops(c, self.physical_qubits) - circuit_length = meta["xval"] - count_ops = [(key, (value, circuit_length)) for key, value in c_count_ops.items()] - meta.update({"count_ops": count_ops}) + + if self.analysis.options.get("gate_error_ratio", None) is None: + # Gate errors are not computed, then counting ops is not necessary. + return transpiled + + # Compute average basis gate numbers per Clifford operation + # This is probably main source of performance regression. + # This should be integrated into transpile pass in future. + for circ in transpiled: + count_ops_result = defaultdict(int) + # This is physical circuits, i.e. qargs is physical index + for inst, qargs, _ in circ.data: + if inst.name in ("measure", "reset", "delay", "barrier", "snapshot"): + continue + qinds = [circ.find_bit(q).index for q in qargs] + if not set(self.physical_qubits).issuperset(qinds): + continue + # Not aware of multi-qubit gate direction + formatted_key = tuple(sorted(qinds)), inst.name + count_ops_result[formatted_key] += 1 + circ.metadata["count_ops"] = tuple(count_ops_result.items()) + return transpiled def _metadata(self): @@ -242,4 +243,5 @@ def _metadata(self): for run_opt in ["meas_level", "meas_return"]: if hasattr(self.run_options, run_opt): metadata[run_opt] = getattr(self.run_options, run_opt) + return metadata diff --git a/qiskit_experiments/library/randomized_benchmarking/rb_utils.py b/qiskit_experiments/library/randomized_benchmarking/rb_utils.py index 488f6c511a..db52b4a554 100644 --- a/qiskit_experiments/library/randomized_benchmarking/rb_utils.py +++ b/qiskit_experiments/library/randomized_benchmarking/rb_utils.py @@ -13,7 +13,6 @@ """ RB Helper functions """ - from typing import Tuple, Dict, Optional, List, Union, Sequence import numpy as np @@ -23,6 +22,7 @@ from qiskit_experiments.database_service.device_component import Qubit from qiskit_experiments.framework import DbAnalysisResultV1, AnalysisResultData +from qiskit_experiments.warnings import deprecated_function class RBUtils: @@ -30,6 +30,13 @@ class RBUtils: from randomized benchmarking experiments""" @staticmethod + @deprecated_function( + last_version="0.4", + msg=( + "This method may return errorneous error ratio. " + "Please directly provide known gate error ratio to the analysis option." + ), + ) def get_error_dict_from_backend( backend: Backend, qubits: Sequence[int] ) -> Dict[Tuple[Sequence[int], str], float]: @@ -62,6 +69,13 @@ def get_error_dict_from_backend( return error_dict @staticmethod + @deprecated_function( + last_version="0.4", + msg=( + "Now this method is integarated into 'StandardRB._transpiled_circuits' method. " + "You don't need to explicitly call this method." + ), + ) def count_ops( circuit: QuantumCircuit, qubits: Optional[Sequence[int]] = None ) -> Dict[Tuple[Sequence[int], str], int]: @@ -94,6 +108,7 @@ def count_ops( return count_ops_result @staticmethod + @deprecated_function(last_version="0.4") def gates_per_clifford( ops_count: List, ) -> Dict[Tuple[Sequence[int], str], float]: @@ -185,6 +200,10 @@ def coherence_limit(nQ=2, T1_list=None, T2_list=None, gatelen=0.1): return coherence_limit_err @staticmethod + @deprecated_function( + last_version="0.4", + msg="Please use calculate_epg function instead. This works regardless of qubit number.", + ) def calculate_1q_epg( epc_1_qubit: Union[float, uncertainties.UFloat], qubits: Sequence[int], @@ -218,6 +237,10 @@ def calculate_1q_epg( return out @staticmethod + @deprecated_function( + last_version="0.4", + msg="Please use calculate_epg function instead. This works regardless of qubit number.", + ) def calculate_2q_epg( epc_2_qubit: Union[uncertainties.UFloat, float], qubits: Sequence[int], diff --git a/releasenotes/notes/cleanup-rb-experiment-f17b6e674ae4e473.yaml b/releasenotes/notes/cleanup-rb-experiment-f17b6e674ae4e473.yaml new file mode 100644 index 0000000000..845c5575f4 --- /dev/null +++ b/releasenotes/notes/cleanup-rb-experiment-f17b6e674ae4e473.yaml @@ -0,0 +1,26 @@ +--- +features: + - | + Curve fit parameter guess function :func:`~rb_decay` has been added. + This improves inital parameter estimation of randomized benchmark experiments. +upgrade: + - | + Computation of error per gates (EPGs) from EPC in :class:`RBAnalysis` has been upgraded. + To compute these values from a single EPC value obtained by the experiment, + we should provide a guess of contribution per basis gate to the depolarization. + This ratio has been extrated from backend properties + with :meth:`RBUtils.get_error_dict_from_backend`, + however, this approach may result in unreproducible EPG outcomes under certain circumstances. + See qiskit-experiments/#762 for more details. + Not this error ratio is provided from hard-coded lookup table, + and user can still provide custom values with analysis option ``gate_error_ratio``. + One can skip computation of EPGs by setting ``False`` to the option. + - | + :class:`RBAnalysis` has been upgraded to compute corrected EPC for 2Q RB. + When the analysis option ``epg_1_qubit`` is provided, + it returns two EPG analysis results, with and without correction for + underlying single qubit depolarization channels. + New result is added under the name ``EPC_corrected``. +deprecations: + - | + Calling :class:`RBUtils` methods have been deprecated and will be removed after 0.4. diff --git a/requirements.txt b/requirements.txt index 4ccdc3d6d2..77a6a7a598 100644 --- a/requirements.txt +++ b/requirements.txt @@ -4,4 +4,3 @@ qiskit-terra>=0.19.0 qiskit-ibmq-provider>=0.16.0 matplotlib>=3.4 uncertainties - diff --git a/test/curve_analysis/test_guess.py b/test/curve_analysis/test_guess.py index 0f7b496817..a93243398b 100644 --- a/test/curve_analysis/test_guess.py +++ b/test/curve_analysis/test_guess.py @@ -186,3 +186,24 @@ def test_baseline_sinusoidal(self, b0, x0, a, freq): b0_guess = guess.constant_sinusoidal_offset(y) self.assertAlmostEqual(b0, b0_guess, delta=0.1) + + @data( + # typical 1Q + [0.5, 0.5, 0.99], + # typical 2Q + [0.25, 0.75, 0.97], + # alpha around equation switching + [0.48, 0.46, 0.85], + # bad limit + [0.20, 0.36, 0.72], + [0.55, 0.40, 0.65], + ) + @unpack + def test_rb_decay(self, a, b, alpha): + """Test of rb decay basis guess.""" + x = np.arange(1, 100, 5) + y = a * alpha**x + b + + alpha_guess = guess.rb_decay(x, y, a=a, b=b) + + self.assertAlmostEqual(alpha, alpha_guess, delta=alpha * 0.1) diff --git a/test/randomized_benchmarking/test_randomized_benchmarking.py b/test/randomized_benchmarking/test_randomized_benchmarking.py index 69e92cd8eb..37abad28cb 100644 --- a/test/randomized_benchmarking/test_randomized_benchmarking.py +++ b/test/randomized_benchmarking/test_randomized_benchmarking.py @@ -23,6 +23,7 @@ from qiskit.quantum_info import Clifford from qiskit_experiments.library import randomized_benchmarking as rb +from qiskit_experiments.database_service.exceptions import DbExperimentEntryNotFound class RBTestCase(QiskitExperimentsTestCase): @@ -86,22 +87,13 @@ def test_single_qubit(self): seed=123, backend=self.backend, ) + exp.analysis.set_options(gate_error_ratio=None) exp.set_transpile_options(**self.transpiler_options) - exp.analysis.set_options(gate_error_ratio={((0,), "sx"): 1.0, ((0,), "rz"): 0.0}) self.assertAllIdentity(exp.circuits()) expdata = exp.run() self.assertExperimentDone(expdata) - # Since this is 1Q RB, the gate error is factored by 0.5 := (2^1 - 1) / 2^1 - # We know the depolarizing parameter since we defined the noise model - sx_epg = expdata.analysis_results("EPG_sx") - rz_epg = expdata.analysis_results("EPG_rz") - - # Allow for 10 percent tolerance - self.assertAlmostEqual(sx_epg.value.n, self.p1q / 2, delta=0.1 * self.p1q) - self.assertAlmostEqual(rz_epg.value.n, self.pvz / 2, delta=0.1 * self.p1q) - # Given we have gate number per Clifford n_gpc, we can compute EPC as # EPC = 1 - (1 - r)^n_gpc # where r is gate error of SX gate, i.e. dep-parameter divided by 2. @@ -122,6 +114,7 @@ def test_two_qubit(self): seed=123, backend=self.backend, ) + exp.analysis.set_options(gate_error_ratio=None) exp.set_transpile_options(**self.transpiler_options) self.assertAllIdentity(exp.circuits()) @@ -147,6 +140,7 @@ def test_add_more_circuit_yields_lower_variance(self): backend=self.backend, num_samples=3, ) + exp1.analysis.set_options(gate_error_ratio=None) exp1.set_transpile_options(**self.transpiler_options) expdata1 = exp1.run() self.assertExperimentDone(expdata1) @@ -158,6 +152,7 @@ def test_add_more_circuit_yields_lower_variance(self): backend=self.backend, num_samples=5, ) + exp2.analysis.set_options(gate_error_ratio=None) exp2.set_transpile_options(**self.transpiler_options) expdata2 = exp2.run() self.assertExperimentDone(expdata2) @@ -387,3 +382,158 @@ def test_expdata_serialization(self): self.assertExperimentDone(expdata) self.assertRoundTripSerializable(expdata, check_func=self.experiment_data_equiv) self.assertRoundTripPickle(expdata, check_func=self.experiment_data_equiv) + + +class TestEPGAnalysis(QiskitExperimentsTestCase): + """Test case for EPG colculation from EPC. + + EPG and depplarizing probability p are assumed to have following relationship + + EPG = (2^n - 1) / 2^n · p + + This p is provided to the Aer noise model, thus we verify EPG computation + by comparing the value with the depolarizing probability. + """ + + def setUp(self): + """Setup the tests.""" + super().setUp() + + # Setup noise model, including more gate for complicated EPG computation + # Note that 1Q channel error is amplified to check 1q channel correction mechanism + x_error = depolarizing_error(0.04, 1) + h_error = depolarizing_error(0.02, 1) + s_error = depolarizing_error(0.00, 1) + cx_error = depolarizing_error(0.08, 2) + + noise_model = NoiseModel() + noise_model.add_all_qubit_quantum_error(x_error, "x") + noise_model.add_all_qubit_quantum_error(h_error, "h") + noise_model.add_all_qubit_quantum_error(s_error, "s") + noise_model.add_all_qubit_quantum_error(cx_error, "cx") + + # Need level1 for consecutive gate cancellation for reference EPC value calculation + transpiler_options = { + "basis_gates": ["x", "h", "s", "cx"], + "optimization_level": 1, + } + + # Aer simulator + backend = AerSimulator(noise_model=noise_model, seed_simulator=123) + + # Prepare experiment data and cache for analysis + exp_1qrb_q0 = rb.StandardRB( + qubits=(0,), + lengths=[1, 10, 30, 50, 80, 120, 150, 200], + seed=123, + backend=backend, + ) + exp_1qrb_q0.set_transpile_options(**transpiler_options) + expdata_1qrb_q0 = exp_1qrb_q0.run(analysis=None).block_for_results(timeout=300) + + exp_1qrb_q1 = rb.StandardRB( + qubits=(1,), + lengths=[1, 10, 30, 50, 80, 120, 150, 200], + seed=123, + backend=backend, + ) + exp_1qrb_q1.set_transpile_options(**transpiler_options) + expdata_1qrb_q1 = exp_1qrb_q1.run(analysis=None).block_for_results(timeout=300) + + exp_2qrb = rb.StandardRB( + qubits=(0, 1), + lengths=[1, 3, 5, 10, 15, 20, 30, 50], + seed=123, + backend=backend, + ) + exp_2qrb.set_transpile_options(**transpiler_options) + expdata_2qrb = exp_2qrb.run(analysis=None).block_for_results(timeout=300) + + self.expdata_1qrb_q0 = expdata_1qrb_q0 + self.expdata_1qrb_q1 = expdata_1qrb_q1 + self.expdata_2qrb = expdata_2qrb + + def test_default_epg_ratio(self): + """Calculate EPG with default ratio dictionary. H and X have the same ratio.""" + analysis = rb.RBAnalysis() + analysis.set_options(outcome="0") + result = analysis.run(self.expdata_1qrb_q0, replace_results=False) + self.assertExperimentDone(result) + + s_epg = result.analysis_results("EPG_s") + h_epg = result.analysis_results("EPG_h") + x_epg = result.analysis_results("EPG_x") + + self.assertEqual(s_epg.value.n, 0.0) + + # H and X gate EPG are assumed to be the same, so this underestimate X and overestimate H + self.assertEqual(h_epg.value.n, x_epg.value.n) + self.assertLess(x_epg.value.n, 0.04 * 0.5) + self.assertGreater(h_epg.value.n, 0.02 * 0.5) + + def test_no_epg(self): + """Calculate no EPGs.""" + analysis = rb.RBAnalysis() + analysis.set_options(outcome="0", gate_error_ratio=None) + result = analysis.run(self.expdata_1qrb_q0, replace_results=False) + self.assertExperimentDone(result) + + with self.assertRaises(DbExperimentEntryNotFound): + result.analysis_results("EPG_s") + + with self.assertRaises(DbExperimentEntryNotFound): + result.analysis_results("EPG_h") + + with self.assertRaises(DbExperimentEntryNotFound): + result.analysis_results("EPG_x") + + def test_with_custom_epg_ratio(self): + """Calculate no EPGs with custom EPG ratio dictionary.""" + analysis = rb.RBAnalysis() + analysis.set_options(outcome="0", gate_error_ratio={"x": 2, "h": 1, "s": 0}) + result = analysis.run(self.expdata_1qrb_q0, replace_results=False) + self.assertExperimentDone(result) + + h_epg = result.analysis_results("EPG_h") + x_epg = result.analysis_results("EPG_x") + + self.assertAlmostEqual(x_epg.value.n, 0.04 * 0.5, delta=0.005) + self.assertAlmostEqual(h_epg.value.n, 0.02 * 0.5, delta=0.005) + + def test_2q_epg(self): + """Compute 2Q EPG without correction. + + Since 1Q gates are designed to have comparable EPG with CX gate, + this will overestimate the error of CX gate. + """ + analysis = rb.RBAnalysis() + analysis.set_options(outcome="00") + result = analysis.run(self.expdata_2qrb, replace_results=False) + self.assertExperimentDone(result) + + cx_epg = result.analysis_results("EPG_cx") + + self.assertGreater(cx_epg.value.n, 0.08 * 0.75) + + def test_correct_1q_depolarization(self): + """Compute 2Q EPG with 1Q depolarization correction.""" + analysis_1qrb_q0 = rb.RBAnalysis() + analysis_1qrb_q0.set_options(outcome="0", gate_error_ratio={"x": 2, "h": 1, "s": 0}) + result_q0 = analysis_1qrb_q0.run(self.expdata_1qrb_q0, replace_results=False) + self.assertExperimentDone(result_q0) + + analysis_1qrb_q1 = rb.RBAnalysis() + analysis_1qrb_q1.set_options(outcome="0", gate_error_ratio={"x": 2, "h": 1, "s": 0}) + result_q1 = analysis_1qrb_q1.run(self.expdata_1qrb_q1, replace_results=False) + self.assertExperimentDone(result_q1) + + analysis_2qrb = rb.RBAnalysis() + analysis_2qrb.set_options( + outcome="00", + epg_1_qubit=result_q0.analysis_results() + result_q1.analysis_results(), + ) + result_2qrb = analysis_2qrb.run(self.expdata_2qrb) + self.assertExperimentDone(result_2qrb) + + cx_epg = result_2qrb.analysis_results("EPG_cx") + self.assertAlmostEqual(cx_epg.value.n, 0.08 * 0.75, delta=0.006) diff --git a/test/randomized_benchmarking/test_rb_utils.py b/test/randomized_benchmarking/test_rb_utils.py index d1d85832ca..c9163e287d 100644 --- a/test/randomized_benchmarking/test_rb_utils.py +++ b/test/randomized_benchmarking/test_rb_utils.py @@ -74,7 +74,8 @@ def test_count_ops(self, num_qubits, expected_counts): rng.shuffle(gates_to_add) for qubits, gate in gates_to_add: circuit.append(self.instructions[gate], qubits) - counts = rb.RBUtils.count_ops(circuit) + with self.assertWarns(DeprecationWarning): + counts = rb.RBUtils.count_ops(circuit) self.assertDictEqual(expected_counts, counts) def test_calculate_1q_epg(self): @@ -87,7 +88,10 @@ def test_calculate_1q_epg(self): qubits = [0] gate_error_ratio = {((0,), "id"): 1, ((0,), "rz"): 0, ((0,), "sx"): 1, ((0,), "x"): 1} gates_per_clifford = {((0,), "rz"): 10.5, ((0,), "sx"): 8.15, ((0,), "x"): 0.25} - epg = rb.RBUtils.calculate_1q_epg(epc_1_qubit, qubits, gate_error_ratio, gates_per_clifford) + with self.assertWarns(DeprecationWarning): + epg = rb.RBUtils.calculate_1q_epg( + epc_1_qubit, qubits, gate_error_ratio, gates_per_clifford + ) error_dict = { ((0,), "rz"): ufloat(0, 0), ((0,), "sx"): ufloat(0.0004432101747785104, 0), @@ -143,9 +147,10 @@ def test_calculate_2q_epg(self): AnalysisResultData("EPG_x", 0.00036207066403884814, device_components=[1]), AnalysisResultData("EPG_x", 0.0005429962529239195, device_components=[4]), ] - epg = rb.RBUtils.calculate_2q_epg( - epc_2_qubit, qubits, gate_error_ratio, gates_per_clifford, epg_1_qubit - ) + with self.assertWarns(DeprecationWarning): + epg = rb.RBUtils.calculate_2q_epg( + epc_2_qubit, qubits, gate_error_ratio, gates_per_clifford, epg_1_qubit + ) error_dict = { ((1, 4), "cx"): ufloat(0.012438847900902494, 0),