Merge pull request #605 from qiboteam/avoided_crossing

Avoided crossing

src/qibocal/auto/operation.py

@@ -9,7 +9,7 @@
 import numpy as np
 import numpy.typing as npt
 from qibolab.platform import Platform
-from qibolab.qubits import Qubit, QubitId, QubitPairId
+from qibolab.qubits import Qubit, QubitId, QubitPair, QubitPairId
 
 from qibocal.config import log
 
@@ -21,7 +21,7 @@
 """Valid value for a routine and runcard parameter."""
 Qubits = dict[QubitId, Qubit]
 """Convenient way of passing qubit pairs in the routines."""
-QubitsPairs = dict[tuple[QubitId, QubitId], Qubit]
+QubitsPairs = dict[tuple[QubitId, QubitId], QubitPair]
 
 
 DATAFILE = "data.npz"

src/qibocal/protocols/characterization/__init__.py

@@ -20,6 +20,7 @@
 from .dispersive_shift_qutrit import dispersive_shift_qutrit
 from .fast_reset.fast_reset import fast_reset
 from .flipping import flipping
+from .flux_dependence.avoided_crossing import avoided_crossing
 from .flux_dependence.qubit_crosstalk import qubit_crosstalk
 from .flux_dependence.qubit_flux_dependence import qubit_flux
 from .flux_dependence.qubit_flux_tracking import qubit_flux_tracking
@@ -105,6 +106,7 @@ class Operation(Enum):
     qubit_spectroscopy_ef = qubit_spectroscopy_ef
     qutrit_classification = qutrit_classification
     resonator_amplitude = resonator_amplitude
+    avoided_crossing = avoided_crossing
     dispersive_shift_qutrit = dispersive_shift_qutrit
     coupler_resonator_spectroscopy = coupler_resonator_spectroscopy
     coupler_qubit_spectroscopy = coupler_qubit_spectroscopy

src/qibocal/protocols/characterization/flux_dependence/avoided_crossing.py

@@ -0,0 +1,326 @@
+from copy import deepcopy
+from dataclasses import dataclass, field
+
+import numpy as np
+import numpy.typing as npt
+import plotly.graph_objects as go
+from plotly.subplots import make_subplots
+from qibolab.platform import Platform
+from qibolab.qubits import QubitId
+
+from qibocal.auto.operation import Data, QubitsPairs, Results, Routine
+from qibocal.protocols.characterization.two_qubit_interaction.utils import order_pair
+from qibocal.protocols.characterization.utils import HZ_TO_GHZ, table_dict, table_html
+
+from .qubit_flux_dependence import QubitFluxParameters, QubitFluxType
+from .qubit_flux_dependence import _acquisition as flux_acquisition
+
+STEP = 60
+POINT_SIZE = 10
+
+
+@dataclass
+class AvoidedCrossingParameters(QubitFluxParameters):
+    """Avoided Crossing Parameters"""
+
+
+@dataclass
+class AvoidedCrossingResults(Results):
+    """Avoided crossing outputs"""
+
+    parabolas: dict[tuple, list]
+    """Extracted parabolas"""
+    fits: dict[tuple, list]
+    """Fits parameters"""
+    cz: dict[tuple, list]
+    """CZ intersection points """
+    iswap: dict[tuple, list]
+    """iSwap intersection points"""
+
+
+@dataclass
+class AvoidedCrossingData(Data):
+    """Avoided crossing acquisition outputs"""
+
+    qubit_pairs: list
+    """list of qubit pairs ordered following the drive frequency"""
+    drive_frequency_low: dict = field(default_factory=dict)
+    """Lowest drive frequency in each qubit pair"""
+    data: dict[tuple[QubitId, str], npt.NDArray[QubitFluxType]] = field(
+        default_factory=dict
+    )
+    """Raw data acquired."""
+
+
+def _acquisition(
+    params: AvoidedCrossingParameters,
+    platform: Platform,
+    qubits: QubitsPairs,  # qubit pairs
+) -> AvoidedCrossingData:
+    """
+    Data acquisition for avoided crossing.
+    This routine performs the qubit flux dependency for the "01" and "02" transition
+    on the qubit pair. It returns the bias and frequency values to perform a CZ
+    and a iSwap gate.
+
+    Args:
+        params (AvoidedCrossingParameters): experiment's parameters.
+        platform (Platform): Qibolab platform object.
+        qubits (dict): list of targets qubit pairs to perform the action.
+    """
+    qubit_pairs = list(qubits.keys())
+    order_pairs = np.array([order_pair(pair, platform.qubits) for pair in qubit_pairs])
+    data = AvoidedCrossingData(qubit_pairs=order_pairs.tolist())
+    # Extract the qubits in the qubits pairs and evaluate their flux dep
+    unique_qubits = np.unique(
+        order_pairs[:, 1]
+    )  # select qubits with high freq in each couple
+    new_qubits = {key: platform.qubits[key] for key in unique_qubits}
+
+    for transition in ["01", "02"]:
+        params.transition = transition
+        data_transition = flux_acquisition(
+            params=params,
+            platform=platform,
+            qubits=new_qubits,
+        )
+        for qubit in unique_qubits:
+            qubit_data = data_transition.data[qubit]
+            freq = qubit_data["freq"]
+            bias = qubit_data["bias"]
+            msr = qubit_data["msr"]
+            phase = qubit_data["phase"]
+            data.register_qubit(
+                QubitFluxType,
+                (float(qubit), transition),
+                dict(
+                    freq=freq.tolist(),
+                    bias=bias.tolist(),
+                    msr=msr.tolist(),
+                    phase=phase.tolist(),
+                ),
+            )
+
+    unique_low_qubits = np.unique(order_pairs[:, 0])
+    data.drive_frequency_low = {
+        str(qubit): float(platform.qubits[qubit].drive_frequency)
+        for qubit in unique_low_qubits
+    }
+    return data
+
+
+def _fit(data: AvoidedCrossingData) -> AvoidedCrossingResults:
+    """
+    Post-Processing for avoided crossing.
+    """
+    qubit_data = data.data
+    fits = {}
+    cz = {}
+    iswap = {}
+    curves = {key: find_parabola(val) for key, val in qubit_data.items()}
+    for qubit_pair in data.qubit_pairs:
+        qubit_pair = tuple(qubit_pair)
+        low = qubit_pair[0]
+        high = qubit_pair[1]
+        # Fit the 02*2 curve
+        curve_02 = np.array(curves[high, "02"]) * 2
+        x_02 = np.unique(qubit_data[high, "02"]["bias"])
+        fit_02 = np.polyfit(x_02, curve_02, 2)
+        fits[qubit_pair, "02"] = fit_02.tolist()
+
+        # Fit the 01+10 curve
+        curve_01 = np.array(curves[high, "01"])
+        x_01 = np.unique(qubit_data[high, "01"]["bias"])
+        fit_01_10 = np.polyfit(x_01, curve_01 + data.drive_frequency_low[str(low)], 2)
+        fits[qubit_pair, "01+10"] = fit_01_10.tolist()
+        # find the intersection of the two parabolas
+        delta_fit = fit_02 - fit_01_10
+        x1, x2 = solve_eq(delta_fit)
+        cz[qubit_pair] = [
+            [x1, np.polyval(fit_02, x1)],
+            [x2, np.polyval(fit_02, x2)],
+        ]
+        # find the intersection of the 01 parabola and the 10 line
+        fit_01 = np.polyfit(x_01, curve_01, 2)
+        fits[qubit_pair, "01"] = fit_01.tolist()
+        fit_pars = deepcopy(fit_01)
+        line_val = data.drive_frequency_low[str(low)]
+        fit_pars[2] -= line_val
+        x1, x2 = solve_eq(fit_pars)
+        iswap[qubit_pair] = [[x1, line_val], [x2, line_val]]
+    return AvoidedCrossingResults(curves, fits, cz, iswap)
+
+
+def _plot(data: AvoidedCrossingData, fit: AvoidedCrossingResults, qubit):
+    """Plotting function for avoided crossing"""
+    fitting_report = ""
+    figures = []
+    order_pair = tuple(index(data.qubit_pairs, qubit))
+    heatmaps = make_subplots(
+        rows=1,
+        cols=2,
+        subplot_titles=[f"{i} transition qubit {qubit[0]}" for i in ["01", "02"]],
+    )
+    parabolas = make_subplots(rows=1, cols=1, subplot_titles=["Parabolas"])
+    for i, transition in enumerate(["01", "02"]):
+        data_high = data.data[order_pair[1], transition]
+        bias_unique = np.unique(data_high.bias)
+        min_bias = min(bias_unique)
+        max_bias = max(bias_unique)
+        heatmaps.add_trace(
+            go.Heatmap(
+                x=data_high.freq * HZ_TO_GHZ,
+                y=data_high.bias,
+                z=data_high.msr,
+                coloraxis="coloraxis",
+            ),
+            row=1,
+            col=i + 1,
+        )
+
+        # the fit of the parabola in 02 transition was done doubling the frequencies
+        heatmaps.add_trace(
+            go.Scatter(
+                x=np.polyval(fit.fits[order_pair, transition], bias_unique)
+                / (i + 1)
+                * HZ_TO_GHZ,
+                y=bias_unique,
+                mode="markers",
+                marker_color="lime",
+                showlegend=True,
+                marker=dict(size=POINT_SIZE),
+                name=f"Curve estimation {transition}",
+            ),
+            row=1,
+            col=i + 1,
+        )
+        heatmaps.add_trace(
+            go.Scatter(
+                x=np.array(fit.parabolas[order_pair[1], transition]) * HZ_TO_GHZ,
+                y=bias_unique,
+                mode="markers",
+                marker_color="black",
+                showlegend=True,
+                marker=dict(symbol="cross", size=POINT_SIZE),
+                name=f"Parabola {transition}",
+            ),
+            row=1,
+            col=i + 1,
+        )
+    cz = np.array(fit.cz[order_pair])
+    iswap = np.array(fit.iswap[order_pair])
+    min_bias = min(min_bias, *cz[:, 0], *iswap[:, 0])
+    max_bias = max(max_bias, *cz[:, 0], *iswap[:, 0])
+    bias_range = np.linspace(min_bias, max_bias, STEP)
+    for transition in ["01", "02", "01+10"]:
+        parabolas.add_trace(
+            go.Scatter(
+                x=bias_range,
+                y=np.polyval(fit.fits[order_pair, transition], bias_range) * HZ_TO_GHZ,
+                showlegend=True,
+                name=transition,
+            )
+        )
+    parabolas.add_trace(
+        go.Scatter(
+            x=bias_range,
+            y=np.array([data.drive_frequency_low[str(order_pair[0])]] * STEP)
+            * HZ_TO_GHZ,
+            showlegend=True,
+            name="10",
+        )
+    )
+    parabolas.add_trace(
+        go.Scatter(
+            x=cz[:, 0],
+            y=cz[:, 1] * HZ_TO_GHZ,
+            showlegend=True,
+            name="CZ",
+            marker_color="black",
+            mode="markers",
+            marker=dict(symbol="cross", size=POINT_SIZE),
+        )
+    )
+    parabolas.add_trace(
+        go.Scatter(
+            x=iswap[:, 0],
+            y=iswap[:, 1] * HZ_TO_GHZ,
+            showlegend=True,
+            name="iswap",
+            marker_color="blue",
+            mode="markers",
+            marker=dict(symbol="cross", size=10),
+        )
+    )
+    parabolas.update_layout(
+        xaxis_title="Bias[V]",
+        yaxis_title="Frequency[GHz]",
+    )
+    heatmaps.update_layout(
+        coloraxis_colorbar=dict(
+            yanchor="top",
+            y=1,
+            x=-0.08,
+            ticks="outside",
+        ),
+        xaxis_title="Frequency[GHz]",
+        yaxis_title="Bias[V]",
+        xaxis2_title="Frequency[GHz]",
+        yaxis2_title="Bias[V]",
+    )
+    figures.append(heatmaps)
+    figures.append(parabolas)
+    fitting_report = table_html(
+        table_dict(
+            qubit,
+            ["CZ bias", "iSwap bias"],
+            [np.round(cz[:, 0], 3), np.round(iswap[:, 0], 3)],
+        )
+    )
+    return figures, fitting_report
+
+
+avoided_crossing = Routine(_acquisition, _fit, _plot)
+
+
+def find_parabola(data: dict) -> list:
+    """
+    Finds the parabola in `data`
+    """
+    freqs = data["freq"]
+    currs = data["bias"]
+    biass = sorted(np.unique(currs))
+    frequencies = []
+    for bias in biass:
+        data_bias = data[currs == bias]
+        index = data_bias["msr"].argmax()
+        average_msr = np.average(data_bias["msr"])
+        st_dev_msr = np.std(data_bias["msr"])
+        if (
+            data_bias["msr"][index] > average_msr + st_dev_msr
+            or data_bias["msr"][index] > average_msr - st_dev_msr
+        ):
+            frequencies.append(freqs[index])
+    return frequencies
+
+
+def solve_eq(pars: list) -> tuple:
+    """
+    Solver of the quadratic equation
+    .. math::
+        a x^2 + b x + c = 0
+    `pars` is the list [a, b, c].
+    """
+    first_term = -1 * pars[1]
+    second_term = np.sqrt(pars[1] ** 2 - 4 * pars[0] * pars[2])
+    x1 = (first_term + second_term) / pars[0] / 2
+    x2 = (first_term - second_term) / pars[0] / 2
+    return x1, x2
+
+
+def index(pairs: list, item: list) -> list:
+    """Find the ordered pair"""
+    for pair in pairs:
+        if set(pair) == set(item):
+            return pair
+    raise ValueError(f"{item} not in pairs")

src/qibocal/protocols/characterization/two_qubit_interaction/cz_virtualz.py

@@ -14,7 +14,7 @@
 from scipy.optimize import curve_fit
 
 from qibocal import update
-from qibocal.auto.operation import Data, Parameters, Qubits, Results, Routine
+from qibocal.auto.operation import Data, Parameters, QubitsPairs, Results, Routine
 from qibocal.config import log
 from qibocal.protocols.characterization.two_qubit_interaction.chevron import order_pair
 from qibocal.protocols.characterization.utils import table_dict, table_html
@@ -147,7 +147,7 @@ def create_sequence(
 def _acquisition(
     params: CZVirtualZParameters,
     platform: Platform,
-    qubits: Qubits,
+    qubits: QubitsPairs,
 ) -> CZVirtualZData:
     r"""
     Acquisition for CZVirtualZ.
@@ -160,14 +160,13 @@ def _acquisition(
     is undone in the high frequency qubit and a theta90 pulse is applied to the low
     frequency qubit before measurement. That is, a pi-half pulse around the relative phase
     parametereized by the angle theta.
-    Measurements on the low frequency qubit yield the the 2Q-phase of the gate and the
+    Measurements on the low frequency qubit yield the 2Q-phase of the gate and the
     remnant single qubit Z phase aquired during the execution to be corrected.
     Population of the high frequency qubit yield the leakage to the non-computational states
     during the execution of the flux pulse.
     """
 
     theta_absolute = np.arange(params.theta_start, params.theta_end, params.theta_step)
-
     data = CZVirtualZData(thetas=theta_absolute.tolist())
     for pair in qubits:
         # order the qubits so that the low frequency one is the first