diff --git a/qiskit_experiments/analysis/plotting.py b/qiskit_experiments/analysis/plotting.py index a64256b712..b9f4150977 100644 --- a/qiskit_experiments/analysis/plotting.py +++ b/qiskit_experiments/analysis/plotting.py @@ -38,7 +38,7 @@ def plot_curve_fit( Wraps ``matplotlib.pyplot.plot``. Args: - func: the fit funcion for curve_fit. + func: the fit function for curve_fit. result: an AnalysisResult from curve_fit. confidence_interval: if True plot the confidence interval from popt_err. ax (matplotlib.axes.Axes): Optional, a matplotlib axes to add the plot to. diff --git a/qiskit_experiments/characterization/__init__.py b/qiskit_experiments/characterization/__init__.py index 3525ea9aa3..7c39b7db9e 100644 --- a/qiskit_experiments/characterization/__init__.py +++ b/qiskit_experiments/characterization/__init__.py @@ -34,3 +34,4 @@ T1Analysis """ from .t1_experiment import T1Experiment, T1Analysis +from .t2star_experiment import T2StarExperiment diff --git a/qiskit_experiments/characterization/t2star_experiment.py b/qiskit_experiments/characterization/t2star_experiment.py new file mode 100644 index 0000000000..7331536c8f --- /dev/null +++ b/qiskit_experiments/characterization/t2star_experiment.py @@ -0,0 +1,267 @@ +# 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. +""" +T2Star Experiment class. +""" + +from typing import List, Optional, Union, Tuple, Dict +import numpy as np + +import qiskit +from qiskit.circuit import QuantumCircuit +from qiskit.utils import apply_prefix +from qiskit_experiments.base_experiment import BaseExperiment +from qiskit_experiments.base_analysis import BaseAnalysis, AnalysisResult +from qiskit_experiments.analysis.curve_fitting import curve_fit, process_curve_data +from qiskit_experiments.analysis.data_processing import level2_probability +from qiskit_experiments.analysis import plotting +from ..experiment_data import ExperimentData + +# pylint: disable = invalid-name +class T2StarAnalysis(BaseAnalysis): + """T2Star Experiment result analysis class.""" + + def __init__( + self, + ): + self._conversion_factor = None + + # pylint: disable=arguments-differ, unused-argument + def _run_analysis( + self, + experiment_data: ExperimentData, + user_p0: Dict[str, float], + user_bounds: Tuple[List[float], List[float]], + plot: bool = True, + ax: Optional["AxesSubplot"] = None, + **kwargs, + ) -> Tuple[AnalysisResult, List["matplotlib.figure.Figure"]]: + r""" + Calculate T2Star experiment + The probability of measuring `+` is assumed to be of the form + .. math:: + f(t) = a\mathrm{e}^{-t / T_2^*}\cos(2\pi freq t + \phi) + b + for unknown parameters :math:`a, b, freq, \phi, T_2^*`. + Args: + experiment_data (ExperimentData): the experiment data to analyze + user_p0: contains initial values given by the user, for the + fit parameters :math:`(a, T_2^*, freq, \phi, b)` + User_bounds: lower and upper bounds on the parameters in p0, + given by the user. + The first tuple is the lower bounds, + The second tuple is the upper bounds. + For both params, the order is :math:`a, T_2^*, freq, \phi, b`. + plot: if True, create the plot, otherwise, do not create the plot. + ax: the plot object + **kwargs: additional parameters + Returns: + The analysis result with the estimated :math:`T_2^*` and 'freq' (frequency) + The graph of the function. + + """ + + def osc_fit_fun(x, a, t2star, freq, phi, c): + """ + Decay cosine fit function + """ + return a * np.exp(-x / t2star) * np.cos(2 * np.pi * freq * x + phi) + c + + def _format_plot(ax, unit): + """Format curve fit plot""" + # Formatting + ax.tick_params(labelsize=10) + ax.set_xlabel("Delay (" + str(unit) + ")", fontsize=12) + ax.set_ylabel("Probability to measure |0>", fontsize=12) + + # implementation of _run_analysis + unit = experiment_data._data[0]["metadata"]["unit"] + self._conversion_factor = experiment_data._data[0]["metadata"].get("dt_factor", None) + if self._conversion_factor is None: + self._conversion_factor = 1 if unit == "s" else apply_prefix(1, unit) + xdata, ydata, sigma = process_curve_data( + experiment_data._data, lambda datum: level2_probability(datum, "0") + ) + + si_xdata = xdata * self._conversion_factor + t2star_estimate = np.mean(si_xdata) + + p0, bounds = self._t2star_default_params(user_p0, user_bounds, t2star_input=t2star_estimate) + fit_result = curve_fit( + osc_fit_fun, si_xdata, ydata, p0=list(p0.values()), sigma=sigma, bounds=bounds + ) + + if plot and plotting.HAS_MATPLOTLIB: + ax = plotting.plot_curve_fit(osc_fit_fun, fit_result, ax=ax) + ax = plotting.plot_scatter(si_xdata, ydata, ax=ax) + ax = plotting.plot_errorbar(si_xdata, ydata, sigma, ax=ax) + _format_plot(ax, unit) + figures = [ax.get_figure()] + else: + figures = None + + # Output unit is 'sec', regardless of the unit used in the input + analysis_result = AnalysisResult( + { + "t2star_value": fit_result["popt"][1], + "frequency_value": fit_result["popt"][2], + "stderr": fit_result["popt_err"][1], + "unit": "s", + "label": "T2*", + "fit": fit_result, + "quality": self._fit_quality( + fit_result["popt"], fit_result["popt_err"], fit_result["reduced_chisq"] + ), + } + ) + + analysis_result["fit"]["circuit_unit"] = unit + if unit == "dt": + analysis_result["fit"]["dt"] = self._conversion_factor + return analysis_result, figures + + def _t2star_default_params( + self, + user_p0, + user_bounds, + t2star_input: float, + ) -> Tuple[List[float], Tuple[List[float]]]: + """ + Default fit parameters for oscillation data + Note that :math:`T_2^*` and 'freq' units are converted to 'sec' and + will be output in 'sec'. + Args: + t2star_input: default for t2star if p0==None + Returns: + Fit guessed parameters: either from the input (if given) or + else assign default values. + """ + if user_p0 is None: + a = 0.5 + t2star = t2star_input * self._conversion_factor + freq = 0.1 + phi = 0.0 + b = 0.5 + else: + a = user_p0["A"] + t2star = user_p0["t2star"] + t2star *= self._conversion_factor + freq = user_p0["f"] + phi = user_p0["phi"] + b = user_p0["B"] + freq /= self._conversion_factor + p0 = {"a_guess": a, "t2star": t2star, "f_guess": freq, "phi_guess": phi, "b_guess": b} + if user_bounds is None: + a_bounds = [-0.5, 1.5] + t2star_bounds = [0, np.inf] + f_bounds = [0.5 * freq, 1.5 * freq] + phi_bounds = [-np.pi, np.pi] + b_bounds = [-0.5, 1.5] + bounds = [ + [a_bounds[i], t2star_bounds[i], f_bounds[i], phi_bounds[i], b_bounds[i]] + for i in range(2) + ] + else: + bounds = user_bounds + return p0, bounds + + @staticmethod + def _fit_quality(fit_out, fit_err, reduced_chisq): + # pylint: disable = too-many-boolean-expressions + if ( + (reduced_chisq < 3) + and (fit_err[0] is None or fit_err[0] < 0.1 * fit_out[0]) + and (fit_err[1] is None or fit_err[1] < 0.1 * fit_out[1]) + and (fit_err[2] is None or fit_err[2] < 0.1 * fit_out[2]) + ): + return "computer_good" + else: + return "computer_bad" + + +class T2StarExperiment(BaseExperiment): + """T2Star experiment class""" + + __analysis_class__ = T2StarAnalysis + + def __init__( + self, + qubit: int, + delays: Union[List[float], np.array], + unit: str = "s", + osc_freq: float = 0.0, + experiment_type: Optional[str] = None, + ): + + """Initialize the T2Star experiment class. + + Args: + qubit: the qubit under test + delays: delay times of the experiments + unit: Optional, time unit of `delays`. + Supported units: 's', 'ms', 'us', 'ns', 'ps', 'dt'. + The unit is used for both T2* and the frequency + osc_freq: the oscillation frequency induced using by the user + experiment_type: String indicating the experiment type. + Can be 'RamseyExperiment' or 'T2StarExperiment'. + """ + + self._qubit = qubit + self._delays = delays + self._unit = unit + self._osc_freq = osc_freq + super().__init__([qubit], experiment_type) + + def circuits( + self, backend: Optional["Backend"] = None, **circuit_options + ) -> List[QuantumCircuit]: + """ + Return a list of experiment circuits + Each circuit consists of a Hadamard gate, followed by a fixed delay, + a phase gate (with a linear phase), and an additional Hadamard gate. + Args: + backend: Optional, a backend object + circuit_options: from base class, empty here + Returns: + The experiment circuits + Raises: + AttributeError: if unit is dt but dt parameter is missing in the backend configuration + """ + if self._unit == "dt": + try: + dt_factor = getattr(backend._configuration, "dt") + except AttributeError as no_dt: + raise AttributeError("Dt parameter is missing in backend configuration") from no_dt + + circuits = [] + for delay in self._delays: + circ = qiskit.QuantumCircuit(1, 1) + circ.h(0) + circ.delay(delay, 0, self._unit) + circ.p(2 * np.pi * self._osc_freq, 0) + circ.barrier(0) + circ.h(0) + circ.barrier(0) + circ.measure(0, 0) + + circ.metadata = { + "experiment_type": self._type, + "qubit": self._qubit, + "osc_freq": self._osc_freq, + "xval": delay, + "unit": self._unit, + } + if self._unit == "dt": + circ.metadata["dt_factor"] = dt_factor + + circuits.append(circ) + + return circuits diff --git a/test/test_t2star.py b/test/test_t2star.py new file mode 100644 index 0000000000..4b23f14d28 --- /dev/null +++ b/test/test_t2star.py @@ -0,0 +1,271 @@ +# (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. +""" +Test T2Star experiment +""" +import unittest +import numpy as np + +from qiskit.utils import apply_prefix +from qiskit.providers import BaseBackend +from qiskit.providers.models import QasmBackendConfiguration +from qiskit.result import Result +from qiskit.test import QiskitTestCase +from qiskit_experiments.composite import ParallelExperiment +from qiskit_experiments.characterization import T2StarExperiment + + +# Fix seed for simulations +SEED = 9000 + + +class T2starBackend(BaseBackend): + """ + A simple and primitive backend, to be run by the T2Star tests + """ + + def __init__( + self, p0=None, initial_prob_plus=None, readout0to1=None, readout1to0=None, dt_factor=1 + ): + """ + Initialize the T2star backend + """ + + configuration = QasmBackendConfiguration( + backend_name="t2star_simulator", + backend_version="0", + n_qubits=int(1e6), + basis_gates=["barrier", "h", "p", "delay", "measure"], + gates=[], + local=True, + simulator=True, + conditional=False, + open_pulse=False, + memory=False, + max_shots=int(1e6), + coupling_map=None, + dt=dt_factor, + ) + + self._t2star = p0["t2star"] + self._a_guess = p0["a_guess"] + self._f_guess = p0["f_guess"] + self._phi_guess = p0["phi_guess"] + self._b_guess = p0["b_guess"] + self._initial_prob_plus = initial_prob_plus + self._readout0to1 = readout0to1 + self._readout1to0 = readout1to0 + self._dt_factor = dt_factor + super().__init__(configuration) + + # pylint: disable = arguments-differ + def run(self, qobj): + """ + Run the T2star backend + """ + shots = qobj.config.shots + result = { + "backend_name": "T2star backend", + "backend_version": "0", + "qobj_id": 0, + "job_id": 0, + "success": True, + "results": [], + } + + for circ in qobj.experiments: + nqubits = circ.config.n_qubits + counts = dict() + if self._readout0to1 is None: + ro01 = np.zeros(nqubits) + else: + ro01 = self._readout0to1 + if self._readout1to0 is None: + ro10 = np.zeros(nqubits) + else: + ro10 = self._readout1to0 + for _ in range(shots): + if self._initial_prob_plus is None: + prob_plus = np.ones(nqubits) + else: + prob_plus = self._initial_prob_plus.copy() + + clbits = np.zeros(circ.config.memory_slots, dtype=int) + for op in circ.instructions: + qubit = op.qubits[0] + + if op.name == "delay": + delay = op.params[0] + t2star = self._t2star[qubit] * self._dt_factor + freq = self._f_guess[qubit] / self._dt_factor + + prob_plus[qubit] = ( + self._a_guess[qubit] + * np.exp(-delay / t2star) + * np.cos(2 * np.pi * freq * delay + self._phi_guess[qubit]) + + self._b_guess[qubit] + ) + + if op.name == "measure": + # we measure in |+> basis which is the same as measuring |0> + meas_res = np.random.binomial( + 1, + (1 - prob_plus[qubit]) * (1 - ro10[qubit]) + + prob_plus[qubit] * ro01[qubit], + ) + clbits[op.memory[0]] = meas_res + + clstr = "" + for clbit in clbits[::-1]: + clstr = clstr + str(clbit) + + if clstr in counts: + counts[clstr] += 1 + else: + counts[clstr] = 1 + result["results"].append( + { + "shots": shots, + "success": True, + "header": {"metadata": circ.header.metadata}, + "data": {"counts": counts}, + } + ) + return Result.from_dict(result) + + +class TestT2Star(QiskitTestCase): + """Test T2Star experiment""" + + def test_t2star_run_end2end(self): + """ + Run the T2 backend on all possible units + """ + # For some reason, 'ps' was not precise enough - need to check this + for unit in ["s", "ms", "us", "ns", "dt"]: + if unit in ("s", "dt"): + dt_factor = 1 + else: + dt_factor = apply_prefix(1, unit) + estimated_t2star = 20 + estimated_freq = 0.1 + # Set up the circuits + qubit = 0 + if unit == "dt": + delays = list(range(1, 46)) + else: + delays = np.append( + (np.linspace(1.0, 15.0, num=15)).astype(float), + (np.linspace(16.0, 45.0, num=59)).astype(float), + ) + + # dummy numbers to avoid exception triggerring + instruction_durations = [ + ("measure", [0], 3, unit), + ("h", [0], 3, unit), + ("p", [0], 3, unit), + ("delay", [0], 3, unit), + ] + + exp = T2StarExperiment(qubit, delays, unit=unit) + + backend = T2starBackend( + p0={ + "a_guess": [0.5], + "t2star": [estimated_t2star], + "f_guess": [estimated_freq], + "phi_guess": [0.0], + "b_guess": [0.5], + }, + initial_prob_plus=[0.0], + readout0to1=[0.02], + readout1to0=[0.02], + dt_factor=dt_factor, + ) + if unit == "dt": + dt_factor = getattr(backend._configuration, "dt") + + # run circuits + result = exp.run( + backend=backend, + user_p0={ + "A": 0.5, + "t2star": estimated_t2star, + "f": estimated_freq, + "phi": 0, + "B": 0.5, + }, + user_bounds=None, + # plot=False, + instruction_durations=instruction_durations, + shots=2000, + ).analysis_result(0) + self.assertAlmostEqual( + result["t2star_value"], + estimated_t2star * dt_factor, + delta=0.08 * result["t2star_value"], + ) + self.assertAlmostEqual( + result["frequency_value"], + estimated_freq / dt_factor, + delta=0.08 * result["frequency_value"], + ) + self.assertEqual( + result["quality"], "computer_good", "Result quality bad for unit " + str(unit) + ) + + def test_t2star_parallel(self): + """ + Test parallel experiments of T2* using a simulator. + """ + + t2star = [30, 25] + estimated_freq = [0.1, 0.12] + delays = [list(range(1, 60)), list(range(1, 50))] + + exp0 = T2StarExperiment(0, delays[0]) + exp2 = T2StarExperiment(2, delays[1]) + par_exp = ParallelExperiment([exp0, exp2]) + + p0 = { + "a_guess": [0.5, None, 0.5], + "t2star": [t2star[0], None, t2star[1]], + "f_guess": [estimated_freq[0], None, estimated_freq[1]], + "phi_guess": [0, None, 0], + "b_guess": [0.5, None, 0.5], + } + backend = T2starBackend(p0) + res = par_exp.run( + backend=backend, + user_p0=None, + user_bounds=None, + # plot=False, + shots=1000, + ) + + for i in range(2): + sub_res = res.component_experiment_data(i).analysis_result(0) + self.assertAlmostEqual( + sub_res["t2star_value"], t2star[i], delta=0.08 * sub_res["t2star_value"] + ) + self.assertAlmostEqual( + sub_res["frequency_value"], + estimated_freq[i], + delta=0.08 * sub_res["frequency_value"], + ) + self.assertEqual( + sub_res["quality"], + "computer_good", + "Result quality bad for experiment on qubit " + str(i), + ) + + +if __name__ == "__main__": + unittest.main()