Merge pull request #651 from qiboteam/fix_fit_ramsey

Fix fit Rabi

src/qibocal/auto/task.py

@@ -156,6 +156,7 @@ def run(
                 ] = platform.settings.relaxation_time
             operation: Routine = self.operation
             parameters = self.parameters
+
         except (RuntimeError, AttributeError):
             operation = dummy_operation
             parameters = DummyPars()

src/qibocal/protocols/characterization/rabi/amplitude.py

@@ -153,7 +153,7 @@ def _fit(data: RabiAmplitudeData) -> RabiAmplitudeResults:
         pguess = [0.5, 0.5, 1 / f, np.pi / 2]
         try:
             popt, perr = curve_fit(
-                utils.rabi_amplitude_fit,
+                utils.rabi_amplitude_function,
                 x,
                 y,
                 p0=pguess,
@@ -179,7 +179,7 @@ def _fit(data: RabiAmplitudeData) -> RabiAmplitudeResults:
         chi2[qubit] = (
             chi2_reduced(
                 y,
-                utils.rabi_amplitude_fit(x, *popt),
+                utils.rabi_amplitude_function(x, *popt),
                 qubit_data.error,
             ),
             np.sqrt(2 / len(y)),

src/qibocal/protocols/characterization/rabi/amplitude_signal.py

@@ -139,7 +139,7 @@ def _fit(data: RabiAmplitudeVoltData) -> RabiAmplitudeVoltResults:
         pguess = [0.5, 1, 1 / f, np.pi / 2]
         try:
             popt, _ = curve_fit(
-                utils.rabi_amplitude_fit,
+                utils.rabi_amplitude_function,
                 x,
                 y,
                 p0=pguess,

src/qibocal/protocols/characterization/rabi/length.py

@@ -154,10 +154,10 @@ def _fit(data: RabiLengthData) -> RabiLengthResults:
         index = local_maxima[0] if len(local_maxima) > 0 else None
         # 0.5 hardcoded guess for less than one oscillation
         f = x[index] / (x[1] - x[0]) if index is not None else 0.5
-        pguess = [0.5, 0.5, np.max(x) / f, np.pi / 2, 0]
+        pguess = [0.5, 0.5, np.max(x) / f, 0, 0]
         try:
             popt, perr = curve_fit(
-                utils.rabi_length_fit,
+                utils.rabi_length_function,
                 x,
                 y,
                 p0=pguess,
@@ -180,7 +180,7 @@ def _fit(data: RabiLengthData) -> RabiLengthResults:
         chi2[qubit] = (
             chi2_reduced(
                 y,
-                utils.rabi_length_fit(x, *popt),
+                utils.rabi_length_function(x, *popt),
                 qubit_data.error,
             ),
             np.sqrt(2 / len(y)),

src/qibocal/protocols/characterization/rabi/length_signal.py

@@ -131,7 +131,7 @@ def _fit(data: RabiLengthVoltData) -> RabiLengthVoltResults:
         x_min = np.min(rabi_parameter)
         x_max = np.max(rabi_parameter)
         x = (rabi_parameter - x_min) / (x_max - x_min)
-        y = (voltages - y_min) / (y_max - y_min)
+        y = (voltages - y_min) / (y_max - y_min) - 1 / 2
 
         # Guessing period using fourier transform
         ft = np.fft.rfft(y)
@@ -141,21 +141,21 @@ def _fit(data: RabiLengthVoltData) -> RabiLengthVoltResults:
         # 0.5 hardcoded guess for less than one oscillation
         f = x[index] / (x[1] - x[0]) if index is not None else 0.5
 
-        pguess = [0.5, 0.5, 1 / f, np.pi / 2, 0]
+        pguess = [0, np.sign(y[0]) * 0.5, 1 / f, 0, 0]
         try:
             popt, _ = curve_fit(
-                utils.rabi_length_fit,
+                utils.rabi_amplitude_function,
                 x,
                 y,
                 p0=pguess,
                 maxfev=100000,
                 bounds=(
-                    [0, 0, 0, -np.pi, 0],
+                    [0, -1, 0, -np.pi, 0],
                     [1, 1, np.inf, np.pi, np.inf],
                 ),
             )
             translated_popt = [  # change it according to the fit function
-                (y_max - y_min) * popt[0] + y_min,
+                (y_max - y_min) * (popt[0] + 1 / 2) + y_min,
                 (y_max - y_min) * popt[1] * np.exp(x_min * popt[4] / (x_max - x_min)),
                 popt[2] * (x_max - x_min),
                 popt[3] - 2 * np.pi * x_min / popt[2] / (x_max - x_min),
@@ -166,7 +166,6 @@ def _fit(data: RabiLengthVoltData) -> RabiLengthVoltResults:
             log.warning("rabi_fit: the fitting was not succesful")
             pi_pulse_parameter = 0
             translated_popt = [0, 0, 1, 0, 0]
-
         durations[qubit] = pi_pulse_parameter
         fitted_parameters[qubit] = translated_popt
 

src/qibocal/protocols/characterization/rabi/utils.py

@@ -5,36 +5,30 @@
 from ..utils import COLORBAND, COLORBAND_LINE, table_dict, table_html
 
 
-def rabi_amplitude_fit(x, p0, p1, p2, p3):
-    # A fit to Superconducting Qubit Rabi Oscillation
-    #   Offset                       : p[0]
-    #   Oscillation amplitude        : p[1]
-    #   Period    T                  : 1/p[2]
-    #   Phase                        : p[3]
-    #   Arbitrary parameter T_2      : 1/p[4]
-    return p0 + p1 * np.sin(2 * np.pi * x / p2 + p3)
-
-
-def rabi_length_fit(x, p0, p1, p2, p3, p4):
-    # A fit to Superconducting Qubit Rabi Oscillation
-    #   Offset                       : p[0]
-    #   Oscillation amplitude        : p[1]
-    #   Period    T                  : 1/p[2]
-    #   Phase                        : p[3]
-    #   Arbitrary parameter T_2      : 1/p[4]
-    return p0 + p1 * np.sin(2 * np.pi * x / p2 + p3) * np.exp(-x * p4)
+def rabi_amplitude_function(x, offset, amplitude, period, phase):
+    """
+    Fit function of Rabi amplitude signal experiment.
 
+    Args:
+        x: Input data.
+    """
+    return offset + amplitude * np.sin(2 * np.pi * x / period + phase)
+
+
+def rabi_length_function(x, offset, amplitude, period, phase, t2_inv):
+    """
+    Fit function of Rabi length signal experiment.
+
+    Args:
+        x: Input data.
+    """
+    return offset + amplitude * np.cos(2 * np.pi * x / period + phase) * np.exp(
+        -x * t2_inv
+    )
 
-def plot(data, qubit, fit):
-    if "RabiAmplitude" in data.__class__.__name__:
-        quantity = "amp"
-        title = "Amplitude [a.u.]"
-        fitting = rabi_amplitude_fit
-    elif data.__class__.__name__ == "RabiLengthVoltData":
-        quantity = "length"
-        title = "Time [ns]"
-        fitting = rabi_length_fit
 
+def plot(data, qubit, fit):
+    quantity, title, fitting = extract_rabi(data)
     figures = []
     fitting_report = ""
 
@@ -118,15 +112,7 @@ def plot(data, qubit, fit):
 
 
 def plot_probabilities(data, qubit, fit):
-    if data.__class__.__name__ == "RabiAmplitudeData":
-        quantity = "amp"
-        title = "Amplitude [a.u.]"
-        fitting = rabi_amplitude_fit
-    elif data.__class__.__name__ == "RabiLengthData":
-        quantity = "length"
-        title = "Time [ns]"
-        fitting = rabi_length_fit
-
+    quantity, title, fitting = extract_rabi(data)
     figures = []
     fitting_report = ""
 
@@ -193,3 +179,14 @@ def plot_probabilities(data, qubit, fit):
     figures.append(fig)
 
     return figures, fitting_report
+
+
+def extract_rabi(data):
+    """
+    Extract Rabi fit info.
+    """
+    if "RabiAmplitude" in data.__class__.__name__:
+        return "amp", "Amplitude [dimensionless]", rabi_amplitude_function
+    if "RabiLength" in data.__class__.__name__:
+        return "length", "Time [ns]", rabi_length_function
+    raise RuntimeError("Data has to be a data structure of the Rabi routines.")

tests/test_protocols.py

@@ -7,6 +7,14 @@
 from qibolab import create_platform
 
 from qibocal.cli._base import command
+from qibocal.protocols.characterization.rabi.amplitude import RabiAmplitudeData
+from qibocal.protocols.characterization.rabi.ef import RabiAmplitudeEFData
+from qibocal.protocols.characterization.rabi.length import RabiLengthData
+from qibocal.protocols.characterization.rabi.utils import (
+    extract_rabi,
+    rabi_amplitude_function,
+    rabi_length_function,
+)
 
 PATH_TO_RUNCARD = pathlib.Path(__file__).parent / "runcards/protocols.yml"
 PLATFORM = create_platform("dummy")
@@ -119,4 +127,19 @@ def test_fit_command(runcard, tmp_path):
     assert results_plot.exit_code == 0
 
 
+def test_extract_rabi():
+    assert extract_rabi(RabiAmplitudeData()) == (
+        "amp",
+        "Amplitude [dimensionless]",
+        rabi_amplitude_function,
+    )
+    assert extract_rabi(RabiLengthData()) == (
+        "length",
+        "Time [ns]",
+        rabi_length_function,
+    )
+    with pytest.raises(RuntimeError):
+        extract_rabi(RabiAmplitudeEFData)
+
+
 # TODO: compare report by calling qq report