diff --git a/qiskit_experiments/curve_analysis/curve_analysis.py b/qiskit_experiments/curve_analysis/curve_analysis.py index 2e15010783..573db1b7ce 100644 --- a/qiskit_experiments/curve_analysis/curve_analysis.py +++ b/qiskit_experiments/curve_analysis/curve_analysis.py @@ -515,6 +515,17 @@ def _extra_database_entry(self, fit_data: FitData) -> List[AnalysisResultData]: """ return [] + def _post_process_fit_result(self, fit_result: FitData) -> FitData: + """A hook that sub-classes can override to manipulate the result of the fit. + + Args: + fit_result: A result from the fitting. + + Returns: + A fit result that might be post-processed. + """ + return fit_result + # pylint: disable=unused-argument def _evaluate_quality(self, fit_data: FitData) -> Union[str, None]: """Evaluate quality of the fit result. @@ -857,6 +868,7 @@ def _run_analysis( fit_result = None else: fit_result = sorted(fit_results, key=lambda r: r.reduced_chisq)[0] + fit_result = self._post_process_fit_result(fit_result) # # 5. Create database entry diff --git a/qiskit_experiments/library/characterization/analysis/drag_analysis.py b/qiskit_experiments/library/characterization/analysis/drag_analysis.py index e3fb887912..c13a873c3e 100644 --- a/qiskit_experiments/library/characterization/analysis/drag_analysis.py +++ b/qiskit_experiments/library/characterization/analysis/drag_analysis.py @@ -25,48 +25,49 @@ class DragCalAnalysis(curve.CurveAnalysis): # section: fit_model - Analyse a Drag calibration experiment by fitting three series each to a cosine function. - The three functions share the phase parameter (i.e. beta) but each have their own amplitude, - baseline, and frequency parameters (which therefore depend on the number of repetitions of - xp-xm). Several initial guesses are tried if the user does not provide one. + Analyse a Drag calibration experiment by fitting three series each to a cosine + function. The three functions share the phase parameter (i.e. beta), amplitude, and + baseline. The frequencies of the oscillations are related through the number of + repetitions of the Drag gates. Several initial guesses are tried if the user + does not provide one. The fit function is .. math:: - y_i = {\rm amp} \cos\left(2 \pi\cdot {\rm freq}_i\cdot x - - 2 \pi\cdot {\rm freq}_i\cdot \beta\right) + {\rm base} + y_i = {\rm amp} \cos\left(2 \pi\cdot {\rm reps}_i \cdot {\rm freq}\cdot x - + 2 \pi\cdot {\rm reps}_i \cdot {\rm freq}\cdot \beta\right) + {\rm base} - Note that the aim of the Drag calibration is to find the :math:`\beta` that minimizes the - phase shifts. This implies that the optimal :math:`\beta` occurs when all three :math:`y` - curves are minimum, i.e. they produce the ground state. Therefore, + Here, the fit parameter :math:`freq` is the frequency of the oscillation of a + single pair of Drag plus and minus rotations and :math:`{\rm reps}_i` is the number + of times that the Drag plus and minus rotations are repeated in curve :math:`i`. + Note that the aim of the Drag calibration is to find the :math:`\beta` that + minimizes the phase shifts. This implies that the optimal :math:`\beta` occurs when + all three :math:`y` curves are minimum, i.e. they produce the ground state. This + occurs when .. math:: - y_i = 0 \quad \Longrightarrow \quad -{\rm amp} \cos(2 \pi\cdot X_i) = {\rm base} + {\rm reps}_i * {\rm freq} * (x - \beta) = N - Here, we abbreviated :math:`{\rm freq}_i\cdot x - {\rm freq}_i\cdot \beta` by :math:`X_i`. - For a signal between 0 and 1 the :math:`{\rm base}` will typically fit to 0.5. However, the - equation has an ambiguity if the amplitude is not properly bounded. Indeed, - - - if :math:`{\rm amp} < 0` then we require :math:`2 \pi\cdot X_i = 0` mod :math:`2\pi`, and - - if :math:`{\rm amp} > 0` then we require :math:`2 \pi\cdot X_i = \pi` mod :math:`2\pi`. - - This will result in an ambiguity in :math:`\beta` which we avoid by bounding the amplitude - from above by 0. + is satisfied with :math:`N` an integer. Note, however, that this condition + produces a minimum only when the amplitude is negative. To ensure this is + the case, we bound the amplitude to be less than 0. # section: fit_parameters defpar \rm amp: desc: Amplitude of all series. - init_guess: The maximum y value less the minimum y value. 0.5 is also tried. - bounds: [-2, 2] scaled to the maximum signal value. + init_guess: The maximum y value scaled by -1, -0.5, and -0.25. + bounds: [-2, 0] scaled to the maximum signal value. defpar \rm base: desc: Base line of all series. - init_guess: The average of the data. 0.5 is also tried. - bounds: [-1, 1] scaled to the maximum signal value. - - defpar {\rm freq}_i: - desc: Frequency of the :math:`i` th oscillation. - init_guess: The frequency with the highest power spectral density. + init_guess: Half the maximum y-value of the data. + bounds: [-1, 1] scaled to the maximum y-value. + + defpar {\rm freq}: + desc: Frequency of oscillation as a function of :math:`\beta` for a single pair + of DRAG plus and minus pulses. + init_guess: For the curve with the most Drag pulse repetitions, the peak frequency + of the power spectral density is found and then divided by the number of repetitions. bounds: [0, inf]. defpar \beta: @@ -77,37 +78,37 @@ class DragCalAnalysis(curve.CurveAnalysis): __series__ = [ curve.SeriesDef( - fit_func=lambda x, amp, freq0, freq1, freq2, beta, base: cos( - x, amp=amp, freq=freq0, phase=-2 * np.pi * freq0 * beta, baseline=base + fit_func=lambda x, amp, freq, reps0, reps1, reps2, beta, base: cos( + x, amp=amp, freq=reps0 * freq, phase=-2 * np.pi * reps0 * freq * beta, baseline=base ), plot_color="blue", name="series-0", filter_kwargs={"series": 0}, plot_symbol="o", - model_description=r"{\rm amp} \cos\left(2 \pi\cdot {\rm freq}_0\cdot x " - r"- 2 \pi\cdot {\rm freq}_0\cdot \beta\right) + {\rm base}", + model_description=r"{\rm amp} \cos\left(2 \pi\cdot {\rm reps}_0\cdot {\rm freq} [x " + r"- \beta]\right) + {\rm base}", ), curve.SeriesDef( - fit_func=lambda x, amp, freq0, freq1, freq2, beta, base: cos( - x, amp=amp, freq=freq1, phase=-2 * np.pi * freq1 * beta, baseline=base + fit_func=lambda x, amp, freq, reps0, reps1, reps2, beta, base: cos( + x, amp=amp, freq=reps1 * freq, phase=-2 * np.pi * reps1 * freq * beta, baseline=base ), plot_color="green", name="series-1", filter_kwargs={"series": 1}, plot_symbol="^", - model_description=r"{\rm amp} \cos\left(2 \pi\cdot {\rm freq}_1\cdot x " - r"- 2 \pi\cdot {\rm freq}_1\cdot \beta\right) + {\rm base}", + model_description=r"{\rm amp} \cos\left(2 \pi\cdot {\rm reps}_1\cdot {\rm freq} [x " + r"- \beta]\right) + {\rm base}", ), curve.SeriesDef( - fit_func=lambda x, amp, freq0, freq1, freq2, beta, base: cos( - x, amp=amp, freq=freq2, phase=-2 * np.pi * freq2 * beta, baseline=base + fit_func=lambda x, amp, freq, reps0, reps1, reps2, beta, base: cos( + x, amp=amp, freq=reps2 * freq, phase=-2 * np.pi * reps2 * freq * beta, baseline=base ), plot_color="red", name="series-2", filter_kwargs={"series": 2}, plot_symbol="v", - model_description=r"{\rm amp} \cos\left(2 \pi\cdot {\rm freq}_2\cdot x " - r"- 2 \pi\cdot {\rm freq}_2\cdot \beta\right) + {\rm base}", + model_description=r"{\rm amp} \cos\left(2 \pi\cdot {\rm reps}_2\cdot {\rm freq} [x " + r"- \beta]\right) + {\rm base}", ), ] @@ -122,6 +123,8 @@ def _default_options(cls): default_options.result_parameters = ["beta"] default_options.xlabel = "Beta" default_options.ylabel = "Signal (arb. units)" + default_options.fixed_parameters = {"reps0": 1, "reps1": 3, "reps2": 5} + default_options.normalization = True return default_options @@ -140,24 +143,31 @@ def _generate_fit_guesses( x_data = self._data("series-0").x min_beta, max_beta = min(x_data), max(x_data) - freqs_guesses = {} - for i in range(3): - curve_data = self._data(f"series-{i}") - freqs_guesses[f"freq{i}"] = curve.guess.frequency(curve_data.x, curve_data.y) - user_opt.p0.set_if_empty(**freqs_guesses) + # Use the highest-frequency curve to estimate the oscillation frequency. + series_label, reps_label = max( + ("series-0", "reps0"), + ("series-1", "reps1"), + ("series-2", "reps2"), + key=lambda x: self.options.fixed_parameters[x[1]], + ) + curve_data = self._data(series_label) + reps2 = self.options.fixed_parameters[reps_label] + freqs_guess = curve.guess.frequency(curve_data.x, curve_data.y) / reps2 + user_opt.p0.set_if_empty(freq=freqs_guess) - max_abs_y, _ = curve.guess.max_height(self._data().y, absolute=True) - freq_bound = max(10 / user_opt.p0["freq0"], max(x_data)) + avg_x = (max(x_data) + min(x_data)) / 2 + span_x = max(x_data) - min(x_data) + beta_bound = max(5 / user_opt.p0["freq"], span_x) + ptp_y = np.ptp(self._data().y) user_opt.bounds.set_if_empty( - amp=(-2 * max_abs_y, 0), - freq0=(0, np.inf), - freq1=(0, np.inf), - freq2=(0, np.inf), - beta=(-freq_bound, freq_bound), - base=(-max_abs_y, max_abs_y), + amp=(-2 * ptp_y, 0), + freq=(0, np.inf), + beta=(avg_x - beta_bound, avg_x + beta_bound), + base=(min(self._data().y) - ptp_y, max(self._data().y) + ptp_y), ) - user_opt.p0.set_if_empty(base=0.5) + base_guess = (max(self._data().y) - min(self._data().y)) / 2 + user_opt.p0.set_if_empty(base=(user_opt.p0["amp"] or base_guess)) # Drag curves can sometimes be very flat, i.e. averages of y-data # and min-max do not always make good initial guesses. We therefore add @@ -165,14 +175,42 @@ def _generate_fit_guesses( # becomes +1 at zero phase, i.e. optimal beta, in which y data should become zero # in discriminated measurement level. options = [] - for amp_guess in (0.5, -0.5): + for amp_factor in (-1, -0.5, -0.25): for beta_guess in np.linspace(min_beta, max_beta, 20): new_opt = user_opt.copy() - new_opt.p0.set_if_empty(amp=amp_guess, beta=beta_guess) + new_opt.p0.set_if_empty(amp=ptp_y * amp_factor, beta=beta_guess) options.append(new_opt) return options + def _post_process_fit_result(self, fit_result: curve.FitData) -> curve.FitData: + r"""Post-process the fit result from a Drag analysis. + + The Drag analysis should return the beta value that is closest to zero. + Since the oscillating term is of the form + + .. math:: + + \cos(2 \pi\cdot {\rm reps}_i \cdot {\rm freq}\cdot [x - \beta]) + + There is a periodicity in beta. This post processing finds the beta that is + closest to zero by performing the minimization using the modulo function. + + .. math:: + + n_\text{min} = \min_{n}|\beta_\text{fit} + n / {\rm freq}| + + and assigning the new beta value to + + .. math:: + + \beta = \beta_\text{fit} + n_\text{min} / {\rm freq}. + """ + beta = fit_result.popt[2] + freq = fit_result.popt[1] + fit_result.popt[2] = ((beta + 1 / freq / 2) % (1 / freq)) - 1 / freq / 2 + return fit_result + def _evaluate_quality(self, fit_data: curve.FitData) -> Union[str, None]: """Algorithmic criteria for whether the fit is good or bad. @@ -182,11 +220,11 @@ def _evaluate_quality(self, fit_data: curve.FitData) -> Union[str, None]: - an error on the drag beta smaller than the beta. """ fit_beta = fit_data.fitval("beta") - fit_freq0 = fit_data.fitval("freq0") + fit_freq = fit_data.fitval("freq") criteria = [ fit_data.reduced_chisq < 3, - fit_beta.nominal_value < 1 / fit_freq0.nominal_value, + abs(fit_beta.nominal_value) < 1 / fit_freq.nominal_value / 2, curve.is_error_not_significant(fit_beta), ] diff --git a/qiskit_experiments/library/characterization/drag.py b/qiskit_experiments/library/characterization/drag.py index a8294949ee..7a9c2031d0 100644 --- a/qiskit_experiments/library/characterization/drag.py +++ b/qiskit_experiments/library/characterization/drag.py @@ -92,14 +92,6 @@ def _default_experiment_options(cls) -> Options: return options - @classmethod - def _default_analysis_options(cls) -> Options: - """Default analysis options.""" - options = Options() - options.normalization = True - - return options - # pylint: disable=arguments-differ def set_experiment_options(self, reps: Optional[List] = None, **fields): """Raise if reps has a length different from three. @@ -107,19 +99,21 @@ def set_experiment_options(self, reps: Optional[List] = None, **fields): Raises: CalibrationError: if the number of repetitions is different from three. """ - - if reps is None: - reps = [1, 3, 5] - else: + if reps is not None: + if len(reps) != 3: + raise CalibrationError( + f"{self.__class__.__name__} must use exactly three repetition numbers. " + f"Received {reps} with length {len(reps)} != 3." + ) reps = sorted(reps) # ensure reps 1 is the lowest frequency. + super().set_experiment_options(reps=reps) - if len(reps) != 3: - raise CalibrationError( - f"{self.__class__.__name__} must use exactly three repetition numbers. " - f"Received {reps} with length {len(reps)} != 3." - ) + if isinstance(self.analysis, DragCalAnalysis): + self.analysis.set_options( + fixed_parameters={"reps0": reps[0], "reps1": reps[1], "reps2": reps[2]} + ) - super().set_experiment_options(reps=reps, **fields) + super().set_experiment_options(**fields) def __init__( self, @@ -143,7 +137,6 @@ def __init__( """ super().__init__([qubit], analysis=DragCalAnalysis(), backend=backend) - self.analysis.set_options(**self._default_analysis_options.__dict__) if betas is not None: self.set_experiment_options(betas=betas) diff --git a/qiskit_experiments/test/mock_iq_backend.py b/qiskit_experiments/test/mock_iq_backend.py index b9fc60eaff..9e194cc401 100644 --- a/qiskit_experiments/test/mock_iq_backend.py +++ b/qiskit_experiments/test/mock_iq_backend.py @@ -291,25 +291,36 @@ def __init__( self, iq_cluster_centers: Tuple[float, float, float, float] = (1.0, 1.0, -1.0, -1.0), iq_cluster_width: float = 1.0, - error: float = 0.03, + freq: float = 0.02, ideal_beta=2.0, gate_name: str = "Rp", rng_seed: int = 0, + max_prob: float = 1.0, + offset_prob: float = 0.0, ): """Initialize the rabi backend.""" - self._error = error + self._freq = freq self._gate_name = gate_name self.ideal_beta = ideal_beta + if max_prob + offset_prob > 1: + raise ValueError("Probabilities need to be between 0 and 1.") + + self._max_prob = max_prob + self._offset_prob = offset_prob + super().__init__(iq_cluster_centers, iq_cluster_width, rng_seed=rng_seed) def _compute_probability(self, circuit: QuantumCircuit) -> float: """Returns the probability based on the beta, number of gates, and leakage.""" - n_gates = sum(circuit.count_ops().values()) + n_gates = circuit.count_ops()[self._gate_name] beta = next(iter(circuit.calibrations[self._gate_name].keys()))[1][0] - return np.sin(n_gates * self._error * (beta - self.ideal_beta)) ** 2 + prob = np.sin(2 * np.pi * n_gates * self._freq * (beta - self.ideal_beta) / 4) ** 2 + rescaled_prob = self._max_prob * prob + self._offset_prob + + return rescaled_prob class RabiBackend(MockIQBackend): diff --git a/test/calibration/experiments/test_drag.py b/test/calibration/experiments/test_drag.py index f70beb1fa4..3c9e2d187f 100644 --- a/test/calibration/experiments/test_drag.py +++ b/test/calibration/experiments/test_drag.py @@ -14,6 +14,7 @@ from test.base import QiskitExperimentsTestCase import unittest +from ddt import ddt, data, unpack import numpy as np from qiskit.circuit import Parameter @@ -30,6 +31,7 @@ from qiskit_experiments.calibration_management import Calibrations +@ddt class TestDragEndToEnd(QiskitExperimentsTestCase): """Test the drag experiment.""" @@ -43,7 +45,7 @@ def setUp(self): pulse.play(Drag(duration=160, amp=0.208519, sigma=40, beta=beta), DriveChannel(0)) self.x_plus = xp - self.test_tol = 0.05 + self.test_tol = 0.1 def test_reps(self): """Test that setting reps raises and error if reps is not of length three.""" @@ -69,10 +71,9 @@ def test_end_to_end(self): # Small leakage will make the curves very flat, in this case one should # rather increase beta. - backend = DragBackend(error=0.0051, gate_name="Drag(xp)") + backend = DragBackend(freq=0.0044, gate_name="Drag(xp)") drag = RoughDrag(0, self.x_plus) - drag.analysis.set_options(p0={"beta": 1.2}) exp_data = drag.run(backend) self.assertExperimentDone(exp_data) result = exp_data.analysis_results(1) @@ -81,11 +82,11 @@ def test_end_to_end(self): self.assertEqual(result.quality, "good") # Large leakage will make the curves oscillate quickly. - backend = DragBackend(error=0.05, gate_name="Drag(xp)") + backend = DragBackend(freq=0.04, gate_name="Drag(xp)") drag = RoughDrag(1, self.x_plus, betas=np.linspace(-4, 4, 31)) drag.set_run_options(shots=200) - drag.analysis.set_options(p0={"beta": 1.8, "freq0": 0.08, "freq1": 0.16, "freq2": 0.32}) + drag.analysis.set_options(p0={"beta": 1.8, "freq": 0.08}) exp_data = drag.run(backend) self.assertExperimentDone(exp_data) result = exp_data.analysis_results(1) @@ -96,6 +97,28 @@ def test_end_to_end(self): self.assertTrue(abs(result.value.n - backend.ideal_beta) < self.test_tol) self.assertEqual(result.quality, "good") + @data( + (0.0040, 1.0, 0.00, [1, 3, 5], None, 0.1), # partial oscillation. + (0.0020, 0.5, 0.00, [1, 3, 5], None, 0.5), # even slower oscillation with amp < 1 + (0.0040, 0.8, 0.05, [3, 5, 7], None, 0.1), # constant offset, i.e. lower SNR. + (0.0800, 0.9, 0.05, [1, 3, 5], np.linspace(-1, 1, 51), 0.1), # Beta not in range + (0.2000, 0.5, 0.10, [1, 3, 5], np.linspace(-2.5, 2.5, 51), 0.1), # Max closer to zero + ) + @unpack + def test_nasty_data(self, freq, amp, offset, reps, betas, tol): + """A set of tests for non-ideal data.""" + backend = DragBackend(freq=freq, gate_name="Drag(xp)", max_prob=amp, offset_prob=offset) + + drag = RoughDrag(0, self.x_plus, betas=betas) + drag.set_experiment_options(reps=reps) + + exp_data = drag.run(backend) + self.assertExperimentDone(exp_data) + result = exp_data.analysis_results("beta") + + self.assertTrue(abs(result.value.n - backend.ideal_beta) < tol) + self.assertEqual(result.quality, "good") + class TestDragCircuits(QiskitExperimentsTestCase): """Test the circuits of the drag calibration.""" @@ -114,7 +137,7 @@ def setUp(self): def test_default_circuits(self): """Test the default circuit.""" - backend = DragBackend(error=0.005, gate_name="Drag(xp)") + backend = DragBackend(freq=0.005, gate_name="Drag(xp)") drag = RoughDrag(0, self.x_plus) drag.set_experiment_options(reps=[2, 4, 8]) diff --git a/test/calibration/experiments/test_fine_drag.py b/test/calibration/experiments/test_fine_drag.py index 812be102ab..1b90429b2f 100644 --- a/test/calibration/experiments/test_fine_drag.py +++ b/test/calibration/experiments/test_fine_drag.py @@ -34,7 +34,7 @@ def _compute_probability(self, circuit: QuantumCircuit) -> float: """Returns the probability based on the beta, number of gates, and leakage.""" n_gates = circuit.count_ops().get("rz", 0) // 2 - return 0.5 * np.sin(n_gates * self._error) + 0.5 + return 0.5 * np.sin(n_gates * self._freq) + 0.5 class TestFineDrag(QiskitExperimentsTestCase):