Drag analysis

qiskit_experiments/library/characterization/analysis/drag_analysis.py

@@ -25,19 +25,23 @@ 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.
 
         .. 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}
 
-        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, :math:`{\rm freq}_i` is given by the fit parameter :math:`freq` multiplied by
+        :math:`{\rm rep}_i` which is the number of times that the Drag plus and minus
+        rotations are repeated. 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,
 
         .. math::
 
@@ -56,17 +60,18 @@ class DragCalAnalysis(curve.CurveAnalysis):
     # 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.
+            init_guess: Half the maximum y-value of the data.
+            bounds: [-1, 1] scaled to the maximum y-value.
 
-        defpar {\rm freq}_i:
-            desc: Frequency of the :math:`i` th oscillation.
-            init_guess: The frequency with the highest power spectral density.
+        defpar {\rm freq}:
+            desc: Frequency of the :math:`i` th oscillation divided by the number of repetitions.
+            init_guess: The frequency with the highest power spectral density of the curve with
+                the highest number of Drag plus and minus rotation repetitions.
             bounds: [0, inf].
 
         defpar \beta:
@@ -77,8 +82,8 @@ 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, reps1, reps2, reps3, beta, base: cos(
+                x, amp=amp, freq=reps1 * freq, phase=-2 * np.pi * reps1 * freq * beta, baseline=base
             ),
             plot_color="blue",
             name="series-0",
@@ -88,8 +93,8 @@ class DragCalAnalysis(curve.CurveAnalysis):
             r"- 2 \pi\cdot {\rm freq}_0\cdot \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, reps1, reps2, reps3, beta, base: cos(
+                x, amp=amp, freq=reps2 * freq, phase=-2 * np.pi * reps2 * freq * beta, baseline=base
             ),
             plot_color="green",
             name="series-1",
@@ -99,8 +104,8 @@ class DragCalAnalysis(curve.CurveAnalysis):
             r"- 2 \pi\cdot {\rm freq}_1\cdot \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, reps1, reps2, reps3, beta, base: cos(
+                x, amp=amp, freq=reps3 * freq, phase=-2 * np.pi * reps3 * freq * beta, baseline=base
             ),
             plot_color="red",
             name="series-2",
@@ -111,6 +116,8 @@ class DragCalAnalysis(curve.CurveAnalysis):
         ),
     ]
 
+    __fixed_parameters__ = ["reps1", "reps2", "reps3"]
+
     @classmethod
     def _default_options(cls):
         """Return the default analysis options.
@@ -122,6 +129,9 @@ def _default_options(cls):
         default_options.result_parameters = ["beta"]
         default_options.xlabel = "Beta"
         default_options.ylabel = "Signal (arb. units)"
+        default_options.rep1 = 1
+        default_options.rep2 = 3
+        default_options.rep3 = 5
 
         return default_options
 
@@ -140,35 +150,32 @@ 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.
+        curve_data = self._data(f"series-2")
+        freqs_guess = curve.guess.frequency(curve_data.x, curve_data.y) / self.options.reps3
+        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))
+        freq_bound = max(10 / user_opt.p0["freq"], max(x_data))
 
         user_opt.bounds.set_if_empty(
             amp=(-2 * max_abs_y, 0),
-            freq0=(0, np.inf),
-            freq1=(0, np.inf),
-            freq2=(0, np.inf),
+            freq=(0, np.inf),
             beta=(-freq_bound, freq_bound),
             base=(-max_abs_y, max_abs_y),
         )
-        user_opt.p0.set_if_empty(base=0.5)
+        user_opt.p0.set_if_empty(base=(user_opt.p0["amp"] or max_abs_y) / 2)
 
         # 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
         # 0.5 to the initial guesses. Note that we also set amp=-0.5 because the cosine function
         # 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=max_abs_y * amp_factor, beta=beta_guess)
                 options.append(new_opt)
 
         return options
@@ -182,11 +189,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,
+            fit_beta.nominal_value < 1 / fit_freq.nominal_value,
             curve.is_error_not_significant(fit_beta),
         ]
 

qiskit_experiments/library/characterization/drag.py

@@ -121,6 +121,11 @@ def set_experiment_options(self, reps: Optional[List] = None, **fields):
 
         super().set_experiment_options(reps=reps, **fields)
 
+        if isinstance(self.analysis, DragCalAnalysis):
+            self.analysis.set_options(reps1=reps[0])
+            self.analysis.set_options(reps2=reps[1])
+            self.analysis.set_options(reps3=reps[2])
+
     def __init__(
         self,
         qubit: int,

qiskit_experiments/test/mock_iq_backend.py

@@ -154,21 +154,32 @@ def __init__(
         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._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(n_gates * self._error * (beta - self.ideal_beta)) ** 2
+        rescaled_prob = self._max_prob * prob + self._offset_prob
+
+        return rescaled_prob
 
 
 class RabiBackend(MockIQBackend):

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."""
@@ -72,7 +74,6 @@ def test_end_to_end(self):
         backend = DragBackend(error=0.0051, 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)
@@ -96,6 +97,27 @@ 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.0050, 1.0, 0.00, [1, 3, 5], None),  # partial oscillation.
+        (0.0025, 0.5, 0.00, [1, 3, 5], None),  # even slower oscillation with amp < 1
+        (0.0050, 0.8, 0.05, [3, 5, 7], None),  # constant offset, i.e. lower SNR.
+        (0.1000, 0.5, 0.10, [1, 3, 5], np.linspace(-1, 1, 51)),  # Beta not in range
+    )
+    @unpack
+    def test_nasty_data(self, error, amp, offset, reps, betas):
+        """A set of tests for non-ideal data."""
+        backend = DragBackend(error=error, 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(1)
+
+        self.assertTrue(abs(result.value.n - backend.ideal_beta) < self.test_tol)
+        self.assertEqual(result.quality, "good")
+
 
 class TestDragCircuits(QiskitExperimentsTestCase):
     """Test the circuits of the drag calibration."""