diff --git a/qualibration_graphs/superconducting/calibration_utils/qubit_spectroscopy/__init__.py b/qualibration_graphs/superconducting/calibration_utils/qubit_spectroscopy/__init__.py index 4c8a5edf09..7b7e8728cb 100644 --- a/qualibration_graphs/superconducting/calibration_utils/qubit_spectroscopy/__init__.py +++ b/qualibration_graphs/superconducting/calibration_utils/qubit_spectroscopy/__init__.py @@ -1,5 +1,14 @@ +"""Qubit spectroscopy calibration utilities (v2).""" + from .parameters import Parameters -from .analysis import process_raw_dataset, fit_raw_data, log_fitted_results +from .analysis import ( + process_raw_dataset, + fit_raw_data, + log_fitted_results, + lorentzian_peak_linbg, + fit_qubit_peak, + FitParameters, +) from .plotting import plot_raw_data_with_fit __all__ = [ @@ -7,5 +16,8 @@ "process_raw_dataset", "fit_raw_data", "log_fitted_results", + "lorentzian_peak_linbg", + "fit_qubit_peak", + "FitParameters", "plot_raw_data_with_fit", ] diff --git a/qualibration_graphs/superconducting/calibration_utils/qubit_spectroscopy/analysis.py b/qualibration_graphs/superconducting/calibration_utils/qubit_spectroscopy/analysis.py index 3b5dabce28..210834b9be 100644 --- a/qualibration_graphs/superconducting/calibration_utils/qubit_spectroscopy/analysis.py +++ b/qualibration_graphs/superconducting/calibration_utils/qubit_spectroscopy/analysis.py @@ -1,154 +1,527 @@ +"""Analysis utilities for qubit spectroscopy calibration (v2). + +Replaces the legacy single-pass `peaks_dips` path with an explicit Lorentzian +peak + linear-background `curve_fit` mirroring the resonator-single fitter. +Adds smart initial guess (Savgol smoothing + `find_peaks` with a noise-floor +prominence threshold), narrow-window refit, and R²/FWHM/contrast quality gates +so weak fits stop silently overwriting QUAM state. + +The IQ rotation step (per-qubit `arctan2`-derived `iw_angle` aligning the +signal onto the rotated I-quadrature) is preserved from the legacy pipeline — +that's the qubit-spec-specific piece, not part of the resonator analysis. +""" + import logging +import warnings from dataclasses import dataclass -from typing import Tuple, Dict +from typing import Dict, Tuple + import numpy as np import xarray as xr +from scipy.optimize import curve_fit +from scipy.signal import find_peaks, savgol_filter from qualibrate import QualibrationNode from qualibration_libs.data import add_amplitude_and_phase, convert_IQ_to_V -from qualibration_libs.analysis import peaks_dips from quam_config.instrument_limits import instrument_limits +# --------------------------------------------------------------------------- +# Fit model +# --------------------------------------------------------------------------- + + +def lorentzian_peak_linbg(f, f0, fwhm, amp, bg0, bg1): + """Lorentzian peak (positive amp) with linear background. + + R(f) = [bg0 + bg1*(f - fc)] + amp / [1 + ((f - f0)/(fwhm/2))^2] + + `fc` is the mean of `f` so `bg0` is the baseline level near the peak, + not a value that depends on the choice of absolute frequency offset. + """ + fc = float(np.mean(f)) + return (bg0 + bg1 * (f - fc)) + amp / (1.0 + ((f - f0) / (fwhm / 2.0)) ** 2) + + +# --------------------------------------------------------------------------- +# Peak finder +# --------------------------------------------------------------------------- + + +def find_best_peak(smoothed, raw_residual_std, edge_fraction=0.04, prominence_factor=3.0): + """Return (peak_idx, is_edge_only) for the highest inner local peak. + + `prominence_factor * raw_residual_std` sets the minimum prominence so we + ignore wiggles that are not statistically above the per-point noise. If + no inner peak is found we fall back to the global argmax with the edge + flag set, so the caller can still reject the fit. + """ + N = len(smoothed) + edge = int(N * edge_fraction) + prominence = max(prominence_factor * raw_residual_std, 1e-12) + peaks, _ = find_peaks(smoothed, prominence=prominence) + inner = [p for p in peaks if edge <= p <= N - 1 - edge] + if inner: + return int(inner[np.argmax(smoothed[inner])]), False + return int(np.argmax(smoothed)), True + + +# --------------------------------------------------------------------------- +# Main fitter +# --------------------------------------------------------------------------- + + +def fit_qubit_peak( + freqs, + signal, + *, + window_fwhm_factor: float = 4.0, + min_window_mhz: float = 5.0, + refit_window_fwhm_factor: float = 8.0, + min_refit_window_mhz: float = 20.0, + detrend_window_mhz: float = 10.0, + max_fwhm_mhz: float = 30.0, + fwhm_floor_hz: float = 5e4, + r2_threshold: float = 0.75, + min_contrast: float = 0.05, + edge_fraction: float = 0.04, + smooth_window: int = 11, +): + """Fit a Lorentzian peak with linear background to `signal(freqs)`. + + Pipeline: Savgol smooth → noise-floor-scaled `find_peaks` → smart f0/FWHM + init via local detrend half-max → first `curve_fit` over ±max(2×FWHM, 2.5 MHz) + → **wider refit** over ±max(4×FWHM, 10 MHz) using the first-pass values → + quality gates (R² inside refit window, FWHM bound, contrast = amp / |baseline|). + Frequencies are in Hz; signal is the rotated I-quadrature in V. + + The refit window is intentionally wider than the first-pass window: it + constrains the linear-bg slope on a representative slice of flat baseline + (a handful of FWHMs on each wing) and makes R² a meaningful number. + + `fwhm_floor_hz` is a hard lower bound for the fitted FWHM. Without it, + `curve_fit` can collapse the Lorentzian onto a single noise spike when the + data has no real peak, producing fake-narrow fits with absurdly high + contrast. Default 50 kHz is well below any real qubit linewidth. + + Returns dict with: f0, fwhm, amp, bg0, bg1, r2, contrast, popt(5), success, + edge_peak, peak_idx, fit_window_half_hz (the refit half-width — useful for + plotting so the overlay only spans the region the fit was actually done on). + """ + nan5 = np.full(5, np.nan) + result = dict( + f0=np.nan, + fwhm=np.nan, + amp=np.nan, + bg0=np.nan, + bg1=np.nan, + r2=np.nan, + contrast=np.nan, + popt=nan5.copy(), + success=False, + edge_peak=False, + peak_idx=-1, + fit_window_half_hz=np.nan, + ) + + freqs = np.asarray(freqs, dtype=float) + signal = np.asarray(signal, dtype=float) + N = len(freqs) + if N < 8: + return result + + span_hz = freqs[-1] - freqs[0] + edge_s = int(N * edge_fraction) + step_hz = float(freqs[1] - freqs[0]) if N > 1 else 1e5 + fwhm_min = max(2.0 * step_hz, fwhm_floor_hz) + + # 1. Smooth + win = min(smooth_window, N // 3 * 2 - 1) + win = win if win % 2 == 1 else win - 1 + win = max(win, 5) + try: + smoothed = savgol_filter(signal, win, 3) + except Exception: + smoothed = signal.copy() + residual_std = float(np.std(signal - smoothed)) + + # 2. Find best inner peak above the noise floor + peak_idx, is_edge = find_best_peak( + smoothed, + residual_std, + edge_fraction=edge_fraction, + ) + result["peak_idx"] = peak_idx + result["edge_peak"] = is_edge + if is_edge or peak_idx < edge_s or peak_idx > N - 1 - edge_s: + return result + + f0_init = freqs[peak_idx] + + # 3. Smart FWHM init via local linear detrend + half-max crossings + detr_half = detrend_window_mhz * 1e6 / 2.0 + detr_mask = (freqs >= f0_init - detr_half) & (freqs <= f0_init + detr_half) + if detr_mask.sum() < 8: + detr_mask = np.ones(N, dtype=bool) + f_d = freqs[detr_mask] + s_d = signal[detr_mask] + b0d = (s_d[0] + s_d[-1]) / 2.0 + b1d = (s_d[-1] - s_d[0]) / (f_d[-1] - f_d[0]) if len(f_d) > 1 else 0.0 + s_detr = s_d - (b0d + b1d * (f_d - f_d.mean())) + pi2 = int(np.argmax(s_detr)) + half_d = s_detr[pi2] / 2.0 + lc = np.where(s_detr[: pi2 + 1] <= half_d)[0] + rc = np.where(s_detr[pi2:] <= half_d)[0] + if len(lc) and len(rc): + fwhm_init = f_d[pi2 + rc[0]] - f_d[lc[-1]] + else: + # Fallback: scan-width tenth + fwhm_init = max(span_hz * 0.1, 1e5) + fwhm_init = max(fwhm_init, 1e5) # floor at 100 kHz + + # 4. First fit window (±N × fwhm_init, clamped) + half_win = max(fwhm_init * window_fwhm_factor / 2.0, min_window_mhz * 1e6 / 2.0) + fit_mask = (freqs >= f0_init - half_win) & (freqs <= f0_init + half_win) + if fit_mask.sum() < 8: + fit_mask = np.ones(N, dtype=bool) + f_win = freqs[fit_mask] + s_win = signal[fit_mask] + + bg0_init = float(np.median(s_win)) + bg1_init = 0.0 + amp_init = float(s_win.max() - bg0_init) + if amp_init <= 0: + amp_init = float(np.std(s_win)) + p0 = [f0_init, max(fwhm_init, fwhm_min), amp_init, bg0_init, bg1_init] + bounds = ( + [f_win.min(), fwhm_min, 0.0, -np.inf, -np.inf], + [f_win.max(), span_hz, max(amp_init * 5, 1e-6) * 10, np.inf, np.inf], + ) + try: + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + popt1, _ = curve_fit( + lorentzian_peak_linbg, + f_win, + s_win, + p0=p0, + bounds=bounds, + maxfev=10000, + ) + except Exception: + return result + + # 5. Wider refit: ±max(refit_window_fwhm_factor/2 × FWHM, min_refit_window/2). + # Wider than the first pass so the linear-bg slope is constrained on a + # representative slice of flat baseline and R² is a meaningful number. + f0_a, fwhm_a, amp_a, bg0_a, bg1_a = popt1 + refit_half = max( + fwhm_a * refit_window_fwhm_factor / 2.0, + min_refit_window_mhz * 1e6 / 2.0, + ) + refit_mask = (freqs >= f0_a - refit_half) & (freqs <= f0_a + refit_half) + if refit_mask.sum() < 8: + refit_mask = fit_mask + f_win = freqs[refit_mask] + s_win = signal[refit_mask] + try: + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + popt, _ = curve_fit( + lorentzian_peak_linbg, + f_win, + s_win, + p0=list(popt1), + bounds=( + [f_win.min(), fwhm_min, 0.0, -np.inf, -np.inf], + [f_win.max(), span_hz, max(amp_a * 5, 1e-6) * 10, np.inf, np.inf], + ), + maxfev=10000, + ) + except Exception: + popt = popt1 # accept first pass if refit fails + # Record the half-width actually used so the plot can clip its overlay + # to the region the fit covered. + result["fit_window_half_hz"] = float(refit_half) + + f0_fit, fwhm_fit, amp_fit, bg0_fit, bg1_fit = popt + + # 6. R² inside the refit window (matches what the plot shows) + y_pred = lorentzian_peak_linbg(f_win, *popt) + ss_res = float(np.sum((s_win - y_pred) ** 2)) + ss_tot = float(np.sum((s_win - s_win.mean()) ** 2)) + r2 = 1.0 - ss_res / ss_tot if ss_tot > 0 else 0.0 + + # 7. Contrast: fitted peak height / fitted baseline level at f0 + bg_at_f0 = bg0_fit + bg1_fit * (f0_fit - f_win.mean()) + contrast = amp_fit / abs(bg_at_f0) if abs(bg_at_f0) > 1e-12 else float("inf") + + # 8. Quality gates + success = ( + r2 >= r2_threshold + and 0 < fwhm_fit <= max_fwhm_mhz * 1e6 + and fwhm_fit / span_hz <= 0.5 + and contrast >= min_contrast + and f_win.min() <= f0_fit <= f_win.max() + ) + + result.update( + f0=float(f0_fit), + fwhm=float(fwhm_fit), + amp=float(amp_fit), + bg0=float(bg0_fit), + bg1=float(bg1_fit), + r2=float(r2), + contrast=float(contrast), + popt=np.array(popt, dtype=float), + success=bool(success), + ) + return result + + +# --------------------------------------------------------------------------- +# FitParameters dataclass +# --------------------------------------------------------------------------- + + @dataclass class FitParameters: - """Stores the relevant qubit spectroscopy experiment fit parameters for a single qubit""" + """Stores the fitted qubit-spectroscopy parameters for a single qubit.""" - frequency: float - relative_freq: float + frequency: float # absolute RF frequency of the qubit + relative_freq: float # peak position in the detuning frame fwhm: float - iw_angle: float + r2: float + contrast: float + iw_angle: float # final integration weight angle = prev + delta (rad) saturation_amp: float x180_amp: float success: bool -def log_fitted_results(fit_results: Dict, log_callable=None): - """ - Logs the node-specific fitted results for all qubits from the fit results +# --------------------------------------------------------------------------- +# Logging +# --------------------------------------------------------------------------- - Parameters: - ----------- - fit_results : dict - Dictionary containing the fitted results for all qubits. - logger : logging.Logger, optional - Logger for logging the fitted results. If None, a default logger is used. - """ +def log_fitted_results(fit_results: Dict, log_callable=None): + """Log the fitted results for all qubits.""" if log_callable is None: log_callable = logging.getLogger(__name__).info - for q in fit_results.keys(): - s_qubit = f"Results for qubit {q}: " - s_freq = f"\tQubit frequency: {1e-9 * fit_results[q]['frequency']:.3f} GHz | " - s_fwhm = f"FWHM: {1e-3 * fit_results[q]['fwhm']:.1f} kHz | " - s_angle = f"The integration weight angle: {fit_results[q]['iw_angle']:.3f} rad\n " - s_saturation = f"To get the desired FWHM, the saturation amplitude is updated to: {1e3 * fit_results[q]['saturation_amp']:.1f} mV | " - s_x180 = f"To get the desired x180 gate, the x180 amplitude is updated to: {1e3 * fit_results[q]['x180_amp']:.1f} mV\n " - if fit_results[q]["success"]: - s_qubit += " SUCCESS!\n" - else: - s_qubit += " FAIL!\n" - log_callable(s_qubit + s_freq + s_fwhm + s_freq + s_angle + s_saturation + s_x180) + for q, res in fit_results.items(): + status = "SUCCESS!" if res["success"] else "FAIL!" + log_callable( + f"Results for qubit {q}: {status}\n" + f"\tQubit frequency: {1e-9 * res['frequency']:.4f} GHz | " + f"FWHM: {1e-3 * res['fwhm']:.1f} kHz | " + f"R²: {res['r2']:.3f} | " + f"contrast: {res['contrast']:.3f}\n" + f"\tIntegration-weight angle: {res['iw_angle']:.3f} rad | " + f"saturation amp: {1e3 * res['saturation_amp']:.2f} mV | " + f"x180 amp: {1e3 * res['x180_amp']:.2f} mV" + ) + + +# --------------------------------------------------------------------------- +# Dataset processing +# --------------------------------------------------------------------------- def process_raw_dataset(ds: xr.Dataset, node: QualibrationNode): + """Convert IQ to V, add amplitude/phase, attach full-frequency coordinate.""" ds = convert_IQ_to_V(ds, node.namespace["qubits"]) ds = add_amplitude_and_phase(ds, "detuning", subtract_slope_flag=True) - full_freq = np.array([ds.detuning + q.xy.RF_frequency for q in node.namespace["qubits"]]) + full_freq = np.array([ds.detuning.values + q.xy.RF_frequency for q in node.namespace["qubits"]]) ds = ds.assign_coords(full_freq=(["qubit", "detuning"], full_freq)) ds.full_freq.attrs = {"long_name": "RF frequency", "units": "Hz"} return ds -def fit_raw_data(ds: xr.Dataset, node: QualibrationNode) -> Tuple[xr.Dataset, dict[str, FitParameters]]: - """ - Fit the qubit frequency and FWHM for each qubit in the dataset. - - Parameters: - ----------- - ds : xr.Dataset - Dataset containing the raw data. - node_parameters : Parameters - Parameters related to the node, including whether state discrimination is used. - - Returns: - -------- - xr.Dataset - Dataset containing the fit results. +# --------------------------------------------------------------------------- +# fit_raw_data +# --------------------------------------------------------------------------- + + +def fit_raw_data(ds: xr.Dataset, node: QualibrationNode) -> Tuple[xr.Dataset, Dict[str, FitParameters]]: + """Fit a Lorentzian peak per qubit on the rotated I-quadrature. + + Per qubit: + 1. arctan2-rotate IQ onto the new I axis (re-used from the legacy node). + 2. Smart-init + first `curve_fit` + narrow-window refit on `I_rot(detuning)`. + 3. Compute R² inside the refit window, plus contrast = amp / |baseline|. + 4. `success = R² > threshold && FWHM in (0, max] && contrast > min`. + 5. Derive the absolute integration-weight angle, saturation amplitude + and x180 amplitude using the existing legacy formulas. """ - ds_fit = ds - # search for frequency for which the amplitude the farthest from the mean to indicate the approximate location of the peak - shifts = np.abs((ds_fit.IQ_abs - ds_fit.IQ_abs.mean(dim="detuning"))).idxmax(dim="detuning") - # Find the rotation angle to align the separation along the 'I' axis - angle = np.arctan2( - ds_fit.sel(detuning=shifts).Q - ds_fit.Q.mean(dim="detuning"), - ds_fit.sel(detuning=shifts).I - ds_fit.I.mean(dim="detuning"), + qubits = node.namespace["qubits"] + limits = [instrument_limits(q.xy) for q in qubits] + + # --- Per-qubit IQ rotation (per legacy analysis.py:91-99) --- + shifts = np.abs((ds.IQ_abs - ds.IQ_abs.mean(dim="detuning"))).idxmax(dim="detuning") + rotation_angle = np.arctan2( + ds.sel(detuning=shifts).Q - ds.Q.mean(dim="detuning"), + ds.sel(detuning=shifts).I - ds.I.mean(dim="detuning"), ) - ds_fit = ds_fit.assign({"iw_angle": angle}) - # rotate the data to the new I axis - ds_fit = ds_fit.assign({"I_rot": ds_fit.I * np.cos(ds_fit.iw_angle) + ds_fit.Q * np.sin(ds_fit.iw_angle)}) - # Find the peak with minimal prominence as defined, if no such peak found, returns nan - fit_vals = peaks_dips(ds_fit.I_rot, dim="detuning", prominence_factor=5) - ds_fit = xr.merge([ds_fit, fit_vals]) - # Extract the relevant fitted parameters - fit_data, fit_results = _extract_relevant_fit_parameters(ds_fit, node) - return fit_data, fit_results - - -def _extract_relevant_fit_parameters(fit: xr.Dataset, node: QualibrationNode): - """Add metadata to the dataset and fit results.""" - limits = [instrument_limits(q.xy) for q in node.namespace["qubits"]] - # Add metadata to fit results - fit.attrs = {"long_name": "frequency", "units": "Hz"} - # Get the fitted resonator frequency - full_freq = np.array([q.xy.RF_frequency for q in node.namespace["qubits"]]) - res_freq = fit.position + full_freq - rel_freq = fit.position - fit = fit.assign({"res_freq": ("qubit", res_freq.data)}) - fit = fit.assign({"relative_freq": ("qubit", rel_freq.data)}) - fit.res_freq.attrs = {"long_name": "qubit xy frequency", "units": "Hz"} - # Get the fitted FWHM - fwhm = np.abs(fit.width) - fit = fit.assign({"fwhm": fwhm}) - fit.fwhm.attrs = {"long_name": "qubit fwhm", "units": "Hz"} - # Get optimum iw angle - prev_angles = np.array( - [q.resonator.operations["readout"].integration_weights_angle for q in node.namespace["qubits"]] + ds_fit = ds.assign({"iw_angle_delta": rotation_angle}) + ds_fit = ds_fit.assign({"I_rot": ds_fit.I * np.cos(rotation_angle) + ds_fit.Q * np.sin(rotation_angle)}) + + # --- Per-qubit Lorentzian peak fit on I_rot vs detuning --- + qubit_names = [q.name for q in qubits] + f0_list, fwhm_list, amp_list, bg0_list, bg1_list = [], [], [], [], [] + r2_list, contrast_list, success_list, popt_list, fit_half_list = ( + [], + [], + [], + [], + [], ) - fit = fit.assign({"iw_angle": (prev_angles + fit.iw_angle) % (2 * np.pi)}) - fit.iw_angle.attrs = {"long_name": "integration weight angle", "units": "rad"} - # Get saturation amplitude - x180_length = np.array([q.xy.operations["x180"].length * 1e-9 for q in node.namespace["qubits"]]) + for q in qubits: + det_q = ds_fit.detuning.values + sig_q = ds_fit.sel(qubit=q.name).I_rot.values + res = fit_qubit_peak( + det_q, + sig_q, + max_fwhm_mhz=node.parameters.max_fwhm_mhz, + r2_threshold=node.parameters.r2_threshold, + min_contrast=node.parameters.min_contrast, + ) + f0_list.append(res["f0"]) + fwhm_list.append(res["fwhm"]) + amp_list.append(res["amp"]) + bg0_list.append(res["bg0"]) + bg1_list.append(res["bg1"]) + r2_list.append(0.0 if np.isnan(res["r2"]) else res["r2"]) + contrast_list.append(0.0 if np.isnan(res["contrast"]) else res["contrast"]) + success_list.append(res["success"]) + popt_list.append(res["popt"]) + fit_half_list.append(res["fit_window_half_hz"]) + + popt_arr = np.stack(popt_list, axis=0) # (n_qubits, 5) + f0_arr = np.array(f0_list) + fwhm_arr = np.array(fwhm_list) + + # --- Derive QUAM-state values (mirrors legacy _extract_relevant_fit_parameters) --- + full_freq = np.array([q.xy.RF_frequency for q in qubits]) + res_freq = f0_arr + full_freq + + prev_angles = np.array([q.resonator.operations["readout"].integration_weights_angle for q in qubits]) + iw_angle_final = (prev_angles + rotation_angle.values) % (2 * np.pi) + used_amp = np.array( - [ - q.xy.operations["saturation"].amplitude * node.parameters.operation_amplitude_factor - for q in node.namespace["qubits"] - ] + [q.xy.operations["saturation"].amplitude * node.parameters.operation_amplitude_factor for q in qubits] + ) + # Treat target_peak_width as the desired FWHM (matches the new fwhm semantics). + with np.errstate(divide="ignore", invalid="ignore"): + factor_cw = node.parameters.target_peak_width / fwhm_arr + sat_amp = factor_cw * used_amp / node.parameters.operation_amplitude_factor + x180_length = np.array([q.xy.operations["x180"].length * 1e-9 for q in qubits]) + factor_x180 = np.pi / (fwhm_arr * x180_length) + x180_amp = factor_x180 * used_amp + + # Tighten success further by the bound checks the legacy node used. + sat_amp_in_bounds = np.array( + [(not np.isnan(sat_amp[i])) and abs(sat_amp[i]) < limits[i].max_wf_amplitude for i in range(len(qubits))] + ) + success_arr = np.array(success_list) & sat_amp_in_bounds + + # --- Pack into ds_fit and FitParameters dict --- + ds_fit = ds_fit.assign( + { + "f0": xr.DataArray( + f0_arr, + coords={"qubit": qubit_names}, + dims="qubit", + attrs={"long_name": "qubit detuning at peak", "units": "Hz"}, + ), + "res_freq": xr.DataArray( + res_freq, + coords={"qubit": qubit_names}, + dims="qubit", + attrs={"long_name": "absolute qubit RF frequency", "units": "Hz"}, + ), + "fwhm": xr.DataArray( + fwhm_arr, + coords={"qubit": qubit_names}, + dims="qubit", + attrs={"long_name": "peak FWHM", "units": "Hz"}, + ), + "amplitude": xr.DataArray( + np.array(amp_list), + coords={"qubit": qubit_names}, + dims="qubit", + ), + "bg0": xr.DataArray( + np.array(bg0_list), + coords={"qubit": qubit_names}, + dims="qubit", + ), + "bg1": xr.DataArray( + np.array(bg1_list), + coords={"qubit": qubit_names}, + dims="qubit", + ), + "r2": xr.DataArray( + np.array(r2_list), + coords={"qubit": qubit_names}, + dims="qubit", + attrs={"long_name": "R²"}, + ), + "contrast": xr.DataArray( + np.array(contrast_list), + coords={"qubit": qubit_names}, + dims="qubit", + attrs={"long_name": "contrast = amp / |baseline|"}, + ), + "popt": xr.DataArray( + popt_arr, + coords={"qubit": qubit_names, "param": np.arange(5)}, + dims=["qubit", "param"], + attrs={"long_name": "fit parameters [f0, fwhm, amp, bg0, bg1]"}, + ), + "fit_window_half_hz": xr.DataArray( + np.array(fit_half_list), + coords={"qubit": qubit_names}, + dims="qubit", + attrs={"long_name": "half-width of the refit window", "units": "Hz"}, + ), + "iw_angle": xr.DataArray( + iw_angle_final, + coords={"qubit": qubit_names}, + dims="qubit", + attrs={"long_name": "integration-weight angle", "units": "rad"}, + ), + "saturation_amplitude": xr.DataArray( + sat_amp, + coords={"qubit": qubit_names}, + dims="qubit", + attrs={"units": "V"}, + ), + "x180_amplitude": xr.DataArray( + x180_amp, + coords={"qubit": qubit_names}, + dims="qubit", + attrs={"units": "V"}, + ), + "success": xr.DataArray( + success_arr, + coords={"qubit": qubit_names}, + dims="qubit", + ), + } ) - factor_cw = node.parameters.target_peak_width / fit.width - fit = fit.assign({"saturation_amplitude": factor_cw * used_amp / node.parameters.operation_amplitude_factor}) - # get expected x180 amplitude - factor_x180 = np.pi / (fit.width * x180_length) - fit = fit.assign({"x180_amplitude": factor_x180 * used_amp}) - - # Assess whether the fit was successful or not - freq_success = np.abs(res_freq) < node.parameters.frequency_span_in_mhz * 1e6 + full_freq - fwhm_success = np.abs(fwhm) < node.parameters.frequency_span_in_mhz * 1e6 + full_freq - saturation_amp_success = np.abs(fit.saturation_amplitude) < limits[0].max_wf_amplitude - # x180amp_success = np.abs(fit.x180_amplitude.data) < limits[0].max_x180_wf_amplitude - success_criteria = freq_success & fwhm_success & saturation_amp_success - fit = fit.assign({"success": success_criteria}) fit_results = { q: FitParameters( - frequency=fit.sel(qubit=q).res_freq.values.__float__(), - relative_freq=fit.sel(qubit=q).relative_freq.values.__float__(), - fwhm=fit.sel(qubit=q).fwhm.values.__float__(), - iw_angle=fit.sel(qubit=q).iw_angle.values.__float__(), - saturation_amp=fit.sel(qubit=q).saturation_amplitude.values.__float__(), - x180_amp=fit.sel(qubit=q).x180_amplitude.values.__float__(), - success=fit.sel(qubit=q).success.values.__bool__(), + frequency=float(ds_fit.sel(qubit=q).res_freq.values), + relative_freq=float(ds_fit.sel(qubit=q).f0.values), + fwhm=float(ds_fit.sel(qubit=q).fwhm.values), + r2=float(ds_fit.sel(qubit=q).r2.values), + contrast=float(ds_fit.sel(qubit=q).contrast.values), + iw_angle=float(ds_fit.sel(qubit=q).iw_angle.values), + saturation_amp=float(ds_fit.sel(qubit=q).saturation_amplitude.values), + x180_amp=float(ds_fit.sel(qubit=q).x180_amplitude.values), + success=bool(ds_fit.sel(qubit=q).success.values), ) - for q in fit.qubit.values + for q in qubit_names } - return fit, fit_results + return ds_fit, fit_results diff --git a/qualibration_graphs/superconducting/calibration_utils/qubit_spectroscopy/node.py b/qualibration_graphs/superconducting/calibration_utils/qubit_spectroscopy/node.py index 5fb88ee78e..aebc67a63d 100644 --- a/qualibration_graphs/superconducting/calibration_utils/qubit_spectroscopy/node.py +++ b/qualibration_graphs/superconducting/calibration_utils/qubit_spectroscopy/node.py @@ -1,3 +1,5 @@ +"""Node utilities for qubit spectroscopy calibration.""" + from typing import List from quam_builder.architecture.superconducting.qubit import AnyTransmon @@ -5,7 +7,7 @@ def get_optional_pulse_duration(qubits: List[AnyTransmon], node_parameters: Parameters): + """Get pulse duration from node parameters or qubit operations.""" if node_parameters.operation_len_in_ns is not None: return {q.name: node_parameters.operation_len_in_ns for q in qubits} - else: - return {q.name: q.xy.operations[node_parameters.operation].length for q in qubits} + return {q.name: q.xy.operations[node_parameters.operation].length for q in qubits} diff --git a/qualibration_graphs/superconducting/calibration_utils/qubit_spectroscopy/parameters.py b/qualibration_graphs/superconducting/calibration_utils/qubit_spectroscopy/parameters.py index 5e1e7d4fb4..1812e9a4ec 100644 --- a/qualibration_graphs/superconducting/calibration_utils/qubit_spectroscopy/parameters.py +++ b/qualibration_graphs/superconducting/calibration_utils/qubit_spectroscopy/parameters.py @@ -1,3 +1,5 @@ +"""Parameters for qubit spectroscopy calibration (v2).""" + from typing import Optional from qualibrate import NodeParameters from qualibrate.core.parameters import RunnableParameters @@ -5,6 +7,8 @@ class NodeSpecificParameters(RunnableParameters): + """Node-specific parameters for qubit spectroscopy.""" + num_shots: int = 100 """Number of averages to perform. Default is 100.""" frequency_span_in_mhz: float = 100 @@ -18,9 +22,20 @@ class NodeSpecificParameters(RunnableParameters): operation_len_in_ns: Optional[int] = None """Length of the operation in nanoseconds. Default is the predefined pulse length.""" target_peak_width: float = 3e6 - """Target peak width in Hz. Default is 3e6 Hz.""" + """Target peak FWHM in Hz used to rescale the saturation amplitude. Default 3 MHz.""" update_pulses_amplitude: bool = False - """Whether to update the saturation pulse and x180/x90 pulse amplitudes based on the peak width. Default is False""" + """Whether to update the saturation pulse and x180/x90 pulse amplitudes based on the peak width. Default False.""" + + # --- Fit quality gates (new in v2) --- + r2_threshold: float = 0.75 + """Minimum coefficient of determination (R²) inside the wider refit window for the + fit to count as successful. Default 0.75 — the refit covers ±max(4×FWHM, 10 MHz) + so a true qubit line plus several FWHM of flat baseline reliably scores > 0.85; + 0.75 gives a safety margin for noise-limited scans. Raise to e.g. 0.9 to be stricter.""" + max_fwhm_mhz: float = 30.0 + """Reject fits whose FWHM exceeds this many MHz (typical qubit lines are 0.1–5 MHz).""" + min_contrast: float = 0.05 + """Reject fits whose fitted peak / fitted baseline contrast is below this fraction.""" class Parameters( @@ -29,4 +44,4 @@ class Parameters( NodeSpecificParameters, QubitsExperimentNodeParameters, ): - pass + """Combined parameters for qubit spectroscopy calibration.""" diff --git a/qualibration_graphs/superconducting/calibration_utils/qubit_spectroscopy/plotting.py b/qualibration_graphs/superconducting/calibration_utils/qubit_spectroscopy/plotting.py index 5b5121106f..d44d937315 100644 --- a/qualibration_graphs/superconducting/calibration_utils/qubit_spectroscopy/plotting.py +++ b/qualibration_graphs/superconducting/calibration_utils/qubit_spectroscopy/plotting.py @@ -1,87 +1,217 @@ +"""Plotting utilities for qubit spectroscopy calibration (v2). + +Per-qubit panel keeps the dual-axis (absolute RF + detuning) layout from the +legacy plot but overlays the explicit Lorentzian+linear-bg fit curve (from +the stored `popt`) and an in-panel values box listing every fit-derived +value that `update_state` may write to QUAM: f₀, FWHM, R², contrast, +integration-weight angle, saturation amplitude, x180 amplitude. + +Failed-fit panels get a red box border and a `(FAILED)` suffix on the title +so bad panels are obvious at a glance. +""" + from typing import List + +import numpy as np import xarray as xr from matplotlib.axes import Axes from matplotlib.figure import Figure from qualang_tools.units import unit from qualibration_libs.plotting import QubitGrid, grid_iter -from qualibration_libs.analysis import lorentzian_peak from quam_builder.architecture.superconducting.qubit import AnyTransmon +from .analysis import lorentzian_peak_linbg + u = unit(coerce_to_integer=True) -def plot_raw_data_with_fit(ds: xr.Dataset, qubits: List[AnyTransmon], fits: xr.Dataset): +_FS_SUPTITLE = 22 +_FS_TITLE = 16 +_FS_LABEL = 14 +_FS_TICK = 12 +_FS_LEGEND = 11 +_FS_BOX = 10 + +_TITLE_PAD = 6 + +# Per-panel sizing. Width per column accommodates the dual axis labels + +# values box; height per row leaves space for title and ticks. +_PER_COL_INCH = 6.0 +_PER_ROW_INCH = 4.5 +_MIN_WIDTH_INCH = 8.0 +_MIN_HEIGHT_INCH = 6.0 + + +def _setup_grid_figure(grid) -> None: + """Resize the QubitGrid figure to its panel count and apply tight_layout. + + constrained_layout was tried but it crashes (ZeroDivisionError in + matplotlib's _layoutgrid) when combined with `ax.axis("off")` cells + that QubitGrid uses for empty grid slots and `bbox_inches="tight"` that + qualang_tools' data handler applies at save time. """ - Plots the resonator spectroscopy amplitude IQ_abs with fitted curves for the given qubits. - - Parameters - ---------- - ds : xr.Dataset - The dataset containing the quadrature data. - qubits : list of AnyTransmon - A list of qubits to plot. - fits : xr.Dataset - The dataset containing the fit parameters. - - Returns - ------- - Figure - The matplotlib figure object containing the plots. - - Notes - ----- - - The function creates a grid of subplots, one for each qubit. - - Each subplot contains the raw data and the fitted curve. + nrows, ncols = grid.all_axes.shape + width = max(_MIN_WIDTH_INCH, ncols * _PER_COL_INCH) + height = max(_MIN_HEIGHT_INCH, nrows * _PER_ROW_INCH) + grid.fig.set_size_inches(width, height) + grid.fig.tight_layout(pad=1.5, w_pad=2.0, h_pad=2.0) + + +def _apply_tick_fontsize(ax: Axes, size: int = _FS_TICK) -> None: + ax.tick_params(axis="both", labelsize=size) + + +def _format_values_box(fit: xr.Dataset, update_pulses: bool) -> str: + """Build the per-panel multi-line values text box. + + Shows only the values that get written back to QUAM (or the + user-facing measurements). R² and contrast were intentionally removed + per user feedback — they are still computed in `analysis.py` and stored + in `ds_fit` / `data.json`; the figure stays uncluttered while the + failure indicator (red border + (FAILED) title) communicates the gate + result without needing the underlying metric on the panel. """ + f0_hz = float(fit.res_freq.values) + fwhm_hz = float(fit.fwhm.values) + iw = float(fit.iw_angle.values) + sat_amp = float(fit.saturation_amplitude.values) + x180_amp = float(fit.x180_amplitude.values) + + lines = [ + f"f₀ = {f0_hz / u.GHz:.5f} GHz", + f"FWHM = {fwhm_hz / u.MHz:.3f} MHz", + f"iw angle = {iw:+.3f} rad", + ] + if update_pulses: + lines.append(f"sat amp = {sat_amp * 1e3:.2f} mV") + lines.append(f"x180 amp = {x180_amp * 1e3:.2f} mV") + return "\n".join(lines) + + +def plot_raw_data_with_fit(ds: xr.Dataset, qubits: List[AnyTransmon], fits: xr.Dataset) -> Figure: + """Per-qubit rotated-I vs detuning panel with Lorentzian fit overlay + values box.""" grid = QubitGrid(ds, [q.grid_location for q in qubits]) + update_pulses = _detect_update_pulses_flag(fits) for ax, qubit in grid_iter(grid): - plot_individual_data_with_fit(ax, ds, qubit, fits.sel(qubit=qubit["qubit"])) - - grid.fig.suptitle("Qubit spectroscopy (rotated 'I' quadrature + fit)") - grid.fig.set_size_inches(15, 9) - grid.fig.tight_layout() + plot_individual_data_with_fit( + ax, + ds, + qubit, + fits.sel(qubit=qubit["qubit"]), + update_pulses, + ) + grid.fig.suptitle( + "Qubit spectroscopy (rotated I + Lorentzian fit)", + fontsize=_FS_SUPTITLE, + ) + _setup_grid_figure(grid) return grid.fig -def plot_individual_data_with_fit(ax: Axes, ds: xr.Dataset, qubit: dict[str, str], fit: xr.Dataset = None): - """ - Plots individual qubit data on a given axis with optional fit. - - Parameters - ---------- - ax : matplotlib.axes.Axes - The axis on which to plot the data. - ds : xr.Dataset - The dataset containing the quadrature data. - qubit : dict[str, str] - mapping to the qubit to plot. - fit : xr.Dataset, optional - The dataset containing the fit parameters (default is None). - - Notes - ----- - - If the fit dataset is provided, the fitted curve is plotted along with the raw data. - """ - if fit: - fitted_data = lorentzian_peak( - ds.detuning, - float(fit.amplitude.values), - float(fit.position.values), - float(fit.width.values) / 2, - float(fit.base_line.mean().values), - ) - else: - fitted_data = None +def _detect_update_pulses_flag(fits: xr.Dataset) -> bool: + """True iff every qubit has a finite saturation_amplitude (proxy for the param).""" + if "saturation_amplitude" not in fits: + return False + vals = fits.saturation_amplitude.values + return bool(np.all(np.isfinite(vals))) + + +def plot_individual_data_with_fit( + ax: Axes, + ds: xr.Dataset, + qubit: dict[str, str], + fit: xr.Dataset, + update_pulses: bool, +) -> None: + """Single-qubit panel with fit overlay, dual axes, and a fit-values text box.""" + ds_q = ds.sel(qubit=qubit["qubit"]) if "qubit" in ds.dims else ds.loc[qubit] + fit_q = fit # already qubit-selected + qubit_name = qubit["qubit"] + + detuning_hz = ds_q.detuning.values + full_freq_hz = ds_q.full_freq.values + i_rot_mv = fit_q.I_rot.values / u.mV if "I_rot" in fit_q else ds_q.I_rot.values / u.mV - # Create a first x-axis for full_freq_GHz - (fit.assign_coords(full_freq_GHz=fit.full_freq / u.GHz).I_rot / u.mV).plot(ax=ax, x="full_freq_GHz") - ax.set_xlabel("RF frequency [GHz]") - ax.set_ylabel("Rotated I [mV]") - # Create a second x-axis for detuning_MHz + # Primary axis: absolute RF [GHz] + ax.plot(full_freq_hz / u.GHz, i_rot_mv, color="steelblue", linewidth=1.0) + ax.set_xlabel("RF frequency [GHz]", fontsize=_FS_LABEL) + ax.set_ylabel("Rotated I [mV]", fontsize=_FS_LABEL) + _apply_tick_fontsize(ax) + + # Secondary x-axis: detuning [MHz] ax2 = ax.twiny() - (fit.assign_coords(detuning_MHz=fit.detuning / u.MHz).I_rot / u.mV).plot(ax=ax2, x="detuning_MHz", label="") - ax2.set_xlabel("Detuning [MHz]") - # Plot the fitted data - if fitted_data is not None: - ax2.plot(fit.detuning / u.MHz, fitted_data / u.mV, "r--") + ax2.set_xlim(detuning_hz[0] / u.MHz, detuning_hz[-1] / u.MHz) + ax2.set_xlabel("Detuning [MHz]", fontsize=_FS_LABEL) + _apply_tick_fontsize(ax2) + + success = bool(fit_q.success.values) if "success" in fit_q else False + popt = fit_q.popt.values if "popt" in fit_q else None + + # Fit overlay — drawn only within the window the fit was actually done over, + # so the linear-bg slope isn't visually extrapolated across the full scan + # (the previous behaviour produced a "tilted tail" that drifted from the + # data's flat baseline even when the fit was correct). + if popt is not None and not np.any(np.isnan(popt)): + f0_rel = float(fit_q.f0.values) if "f0" in fit_q else float(popt[0]) + fwhm = float(fit_q.fwhm.values) if "fwhm" in fit_q else float(popt[1]) + if "fit_window_half_hz" in fit_q and np.isfinite(float(fit_q.fit_window_half_hz.values)): + window_half = float(fit_q.fit_window_half_hz.values) + else: + window_half = max(4.0 * fwhm, 10e6) + mask = np.abs(detuning_hz - f0_rel) <= window_half + if mask.sum() >= 4: + fit_curve_mv = lorentzian_peak_linbg(detuning_hz[mask], *popt) / u.mV + ax.plot( + full_freq_hz[mask] / u.GHz, + fit_curve_mv, + "r--", + linewidth=1.5, + label="Lorentzian fit", + ) + # f₀ marker + f0_hz = float(fit_q.res_freq.values) + ax.axvline( + f0_hz / u.GHz, + color="red", + linestyle=":", + linewidth=1.0, + alpha=0.7, + ) + ax.legend(fontsize=_FS_LEGEND, loc="upper left") + + # Title (with FAILED suffix if appropriate) + title = qubit_name + (" (FAILED)" if not success else "") + ax.set_title(title, loc="center", fontsize=_FS_TITLE, pad=_TITLE_PAD, color=("crimson" if not success else "black")) + + # Values box anchored upper-right; red border on failure. + if popt is not None and not np.any(np.isnan(popt)): + text = _format_values_box(fit_q, update_pulses) + edge = "crimson" if not success else "gray" + ax.text( + 0.98, + 0.98, + text, + transform=ax.transAxes, + ha="right", + va="top", + fontsize=_FS_BOX, + family="monospace", + bbox=dict( + facecolor="white", + alpha=0.88, + edgecolor=edge, + linewidth=1.5, + ), + ) + else: + ax.text( + 0.5, + 0.5, + "NO FIT", + transform=ax.transAxes, + ha="center", + va="center", + fontsize=_FS_LABEL, + color="gray", + ) diff --git a/qualibration_graphs/superconducting/calibrations/1Q_calibrations/03a_qubit_spectroscopy.py b/qualibration_graphs/superconducting/calibrations/1Q_calibrations/03a_qubit_spectroscopy.py index 4aef904176..8723227c3b 100644 --- a/qualibration_graphs/superconducting/calibrations/1Q_calibrations/03a_qubit_spectroscopy.py +++ b/qualibration_graphs/superconducting/calibrations/1Q_calibrations/03a_qubit_spectroscopy.py @@ -1,8 +1,20 @@ +"""Qubit spectroscopy calibration (v2: improved Lorentzian-fit analysis). + +Same QUA program as the previous 03_qubit_spectroscopy node; uses the +improved analysis package `calibration_utils.qubit_spectroscopy` which +replaces the single-pass `peaks_dips` path with an explicit Lorentzian-peak ++ linear-bg `curve_fit`, smart initial guess, narrow-window refit, and +R²/FWHM/contrast quality gates. Every fit-derived value is shown in the plot +panel; failed fits get a red "(FAILED)" badge so update_state does not +silently overwrite QUAM state. +""" + # %% {Imports} -import matplotlib.pyplot as plt +from dataclasses import asdict + import numpy as np import xarray as xr -from dataclasses import asdict +import matplotlib.pyplot as plt from qm.qua import * @@ -12,7 +24,10 @@ from qualang_tools.units import unit from qualibrate import QualibrationNode -from quam_config import Quam +from qualibration_libs.parameters import get_qubits +from qualibration_libs.runtime import simulate_and_plot +from qualibration_libs.data import XarrayDataFetcher + from calibration_utils.qubit_spectroscopy import ( Parameters, process_raw_dataset, @@ -20,89 +35,95 @@ log_fitted_results, plot_raw_data_with_fit, ) -from qualibration_libs.parameters import get_qubits -from qualibration_libs.runtime import simulate_and_plot -from qualibration_libs.data import XarrayDataFetcher +from quam_config import Quam -# %% {Node initialisation} -description = """ - QUBIT SPECTROSCOPY -This sequence involves sending a saturation pulse to the qubit, placing it in a mixed state, -and then measuring the state of the resonator across various qubit drive frequencies. -In order to facilitate the qubit search, the qubit pulse duration and amplitude can be changed manually -from the node parameters. -The data is post-processed to determine the qubit resonance frequency and the width of the peak. +import sys +import os -Note that it can happen that the qubit is excited by the image sideband or LO leakage instead of the desired sideband. -This is why calibrating the qubit mixer is highly recommended when using external mixers or the Octave. +sys.path.insert(0, os.path.join(os.path.dirname(__file__), "../../../..")) + + +# %% {Node initialisation} +description = """ + QUBIT SPECTROSCOPY (v2) +A saturation pulse drives the qubit across a detuning sweep +around its stored RF frequency; the resonator state is read +out and demodulated to I/Q for every drive frequency. + +The post-processing now uses an explicit Lorentzian-peak + linear-background +`scipy.optimize.curve_fit` instead of the legacy `peaks_dips` call: + - Savgol-smoothed `find_peaks` with a noise-floor prominence threshold for + the initial f0 guess. + - Local linear detrend + half-max crossings for a smart FWHM init. + - First fit on a ±4×FWHM window, then a narrow-window refit using the + first-pass values. + - Quality gates: R² (default >= 0.85), FWHM bound (default 30 MHz), + contrast = amp / |baseline| (default >= 0.05), plus the legacy saturation + amplitude bound check. + +Failed fits set node.outcomes[q.name] = "failed" so update_state skips them. Prerequisites: - - Having calibrated the mixer or the Octave (nodes 01a or 01b). - - Having calibrated the readout parameters (nodes 02a, 02b and/or 02c). - - Having specified the desired flux point if relevant (qubit.z.flux_point). - -State update: - - The qubit 0->1 frequency: qubit.f_01 & qubit.xy.RF_frequency - - The integration weight angle to get the state discrimination along the 'I' quadrature: qubit.resonator.operations["readout"].integration_weights_angle. - - (optional) The saturation pulse amplitude to get the targeted fwhm: qubit.xy.operations["saturation"].amplitude. - - (optional) The guessed x180/x90 pulse amplitude: qubit.xy.operations["x180"/"x90"].amplitude. + - Mixer / Octave calibrated (01a or 01b). + - Readout calibrated (02a/02b/02c). + - Flux point specified (qubit.z.flux_point) if relevant. + +State update (only when fit passes the new quality gates): + - qubit.f_01 & qubit.xy.RF_frequency + - qubit.resonator.operations["readout"].integration_weights_angle + - (optional) qubit.xy.operations["saturation"].amplitude + - (optional) qubit.xy.operations["x180"/"x90"].amplitude """ -# Be sure to include [Parameters, Quam] so the node has proper type hinting node = QualibrationNode[Parameters, Quam]( - name="03a_qubit_spectroscopy", # Name should be unique - description=description, # Describe what the node is doing, which is also reflected in the QUAlibrate GUI - parameters=Parameters(), # Node parameters defined under quam_experiment/experiments/node_name + name="03a_qubit_spectroscopy", + description=description, + parameters=Parameters(), machine=Quam.load(), ) -# Any parameters that should change for debugging purposes only should go in here -# These parameters are ignored when run through the GUI or as part of a graph @node.run_action(skip_if=node.modes.external) def custom_param(node: QualibrationNode[Parameters, Quam]): - """Allow the user to locally set the node parameters for debugging purposes, or execution in the Python IDE.""" - # You can get type hinting in your IDE by typing node.parameters. - # node.parameters.qubits = ["q1", "q2"] + """Allow the user to locally set the node parameters for debugging.""" pass +# Instantiate the QUAM class from the state file +# node.machine = Quam.load() + + # %% {Create_QUA_program} @node.run_action(skip_if=node.parameters.load_data_id is not None) def create_qua_program(node: QualibrationNode[Parameters, Quam]): - """Create the sweep axes and generate the QUA program from the pulse sequence and the node parameters.""" - # Class containing tools to help handle units and conversions. + """Create the sweep axes and generate the QUA program.""" u = unit(coerce_to_integer=True) - # Get the active qubits from the node and organize them by batches node.namespace["qubits"] = qubits = get_qubits(node) num_qubits = len(qubits) - operation = node.parameters.operation # The qubit operation to play - n_avg = node.parameters.num_shots # The number of averages - # Adjust the pulse duration and amplitude to drive the qubit into a mixed state - can be None + operation = node.parameters.operation + n_avg = node.parameters.num_shots operation_len = node.parameters.operation_len_in_ns - # pre-factor to the value defined in the config - restricted to [-2; 2) operation_amp = node.parameters.operation_amplitude_factor - # Qubit detuning sweep with respect to their resonance frequencies span = node.parameters.frequency_span_in_mhz * u.MHz step = node.parameters.frequency_step_in_mhz * u.MHz dfs = np.arange(-span // 2, +span // 2, step) - # Register the sweep axes to be added to the dataset when fetching data node.namespace["sweep_axes"] = { "qubit": xr.DataArray(qubits.get_names()), - "detuning": xr.DataArray(dfs, attrs={"long_name": "readout frequency", "units": "Hz"}), + "detuning": xr.DataArray( + dfs, + attrs={"long_name": "readout frequency", "units": "Hz"}, + ), } with program() as node.namespace["qua_program"]: - # Macro to declare I, Q, n and their respective streams for a given number of qubit I, I_st, Q, Q_st, n, n_st = node.machine.declare_qua_variables() - df = declare(int) # QUA variable for the qubit frequency + df = declare(int) for multiplexed_qubits in qubits.batch(): - # Initialize the QPU in terms of flux points (flux tunable transmons and/or tunable couplers) for qubit in multiplexed_qubits.values(): node.machine.initialize_qpu(target=qubit) align() @@ -111,11 +132,8 @@ def create_qua_program(node: QualibrationNode[Parameters, Quam]): save(n, n_st) with for_(*from_array(df, dfs)): for i, qubit in multiplexed_qubits.items(): - # Get the duration of the operation from the node parameters or the state duration = operation_len if operation_len is not None else qubit.xy.operations[operation].length - # Update the qubit frequency qubit.xy.update_frequency(df + qubit.xy.intermediate_frequency) - # Play the saturation pulse qubit.xy.play( operation, amplitude_scale=operation_amp, @@ -124,11 +142,8 @@ def create_qua_program(node: QualibrationNode[Parameters, Quam]): align() for i, qubit in multiplexed_qubits.items(): - # readout the resonator qubit.resonator.measure("readout", qua_vars=(I[i], Q[i])) - # wait for the resonator to deplete qubit.resonator.wait(node.machine.depletion_time * u.ns) - # save data save(I[i], I_st[i]) save(Q[i], Q_st[i]) align() @@ -143,30 +158,26 @@ def create_qua_program(node: QualibrationNode[Parameters, Quam]): # %% {Simulate} @node.run_action(skip_if=node.parameters.load_data_id is not None or not node.parameters.simulate) def simulate_qua_program(node: QualibrationNode[Parameters, Quam]): - """Connect to the QOP and simulate the QUA program""" - # Connect to the QOP + """Connect to the QOP and simulate the QUA program.""" qmm = node.machine.connect() - # Get the config from the machine config = node.machine.generate_config() - # Simulate the QUA program, generate the waveform report and plot the simulated samples - samples, fig, wf_report = simulate_and_plot(qmm, config, node.namespace["qua_program"], node.parameters) - # Store the figure, waveform report and simulated samples + samples, fig, wf_report = simulate_and_plot( + qmm, + config, + node.namespace["qua_program"], + node.parameters, + ) node.results["simulation"] = {"figure": fig, "wf_report": wf_report, "samples": samples} # %% {Execute} @node.run_action(skip_if=node.parameters.load_data_id is not None or node.parameters.simulate) def execute_qua_program(node: QualibrationNode[Parameters, Quam]): - """Connect to the QOP, execute the QUA program and fetch the raw data and store it in a xarray dataset called "ds_raw".""" - # Connect to the QOP + """Execute the QUA program and fetch the raw dataset.""" qmm = node.machine.connect() - # Get the config from the machine config = node.machine.generate_config() - # Execute the QUA program only if the quantum machine is available (this is to avoid interrupting running jobs). with qm_session(qmm, config, timeout=node.parameters.timeout) as qm: - # The job is stored in the node namespace to be reused in the fetching_data run_action node.namespace["job"] = job = qm.execute(node.namespace["qua_program"]) - # Display the progress bar data_fetcher = XarrayDataFetcher(job, node.namespace["sweep_axes"]) for dataset in data_fetcher: progress_counter( @@ -174,9 +185,7 @@ def execute_qua_program(node: QualibrationNode[Parameters, Quam]): node.parameters.num_shots, start_time=data_fetcher.t_start, ) - # Display the execution report to expose possible runtime errors node.log(job.execution_report()) - # Register the raw dataset node.results["ds_raw"] = dataset @@ -185,22 +194,19 @@ def execute_qua_program(node: QualibrationNode[Parameters, Quam]): def load_data(node: QualibrationNode[Parameters, Quam]): """Load a previously acquired dataset.""" load_data_id = node.parameters.load_data_id - # Load the specified dataset node.load_from_id(node.parameters.load_data_id) node.parameters.load_data_id = load_data_id - # Get the active qubits from the loaded node parameters node.namespace["qubits"] = get_qubits(node) # %% {Analyse_data} @node.run_action(skip_if=node.parameters.simulate) def analyse_data(node: QualibrationNode[Parameters, Quam]): - """Analyse the raw data and store the fitted data in another xarray dataset "ds_fit" and the fitted results in the "fit_results" dictionary.""" + """Smart-init + curve_fit + narrow refit + R² gates.""" node.results["ds_raw"] = process_raw_dataset(node.results["ds_raw"], node) node.results["ds_fit"], fit_results = fit_raw_data(node.results["ds_raw"], node) node.results["fit_results"] = {k: asdict(v) for k, v in fit_results.items()} - # Log the relevant information extracted from the data analysis log_fitted_results(node.results["fit_results"], log_callable=node.log) node.outcomes = { qubit_name: ("successful" if fit_result["success"] else "failed") @@ -211,35 +217,30 @@ def analyse_data(node: QualibrationNode[Parameters, Quam]): # %% {Plot_data} @node.run_action(skip_if=node.parameters.simulate) def plot_data(node: QualibrationNode[Parameters, Quam]): - """Plot the raw and fitted data in specific figures whose shape is given by qubit.grid_location.""" - fig_raw_fit = plot_raw_data_with_fit(node.results["ds_raw"], node.namespace["qubits"], node.results["ds_fit"]) + """Per-qubit panel with fit overlay + values box + FAILED indicator.""" + fig_raw_fit = plot_raw_data_with_fit( + node.results["ds_raw"], + node.namespace["qubits"], + node.results["ds_fit"], + ) plt.show() - # Store the generated figures - node.results["figures"] = { - "amplitude": fig_raw_fit, - } + node.results["figures"] = {"amplitude": fig_raw_fit} # %% {Update_state} @node.run_action(skip_if=node.parameters.simulate) def update_state(node: QualibrationNode[Parameters, Quam]): - """Update the relevant parameters if the qubit data analysis was successful.""" + """Write fitted values to QUAM only for qubits whose fit passed quality gates.""" with node.record_state_updates(): for q in node.namespace["qubits"]: if node.outcomes[q.name] == "failed": continue - - # Update the readout frequency for the given flux point - q.f_01 = node.results["fit_results"][q.name]["frequency"] - q.xy.RF_frequency = node.results["fit_results"][q.name]["frequency"] - fit_result = node.results["fit_results"][q.name] - # Update the integration weight angle + q.f_01 = fit_result["frequency"] + q.xy.RF_frequency = fit_result["frequency"] q.resonator.operations["readout"].integration_weights_angle = fit_result["iw_angle"] if node.parameters.update_pulses_amplitude: - # Update the saturation amplitude q.xy.operations["saturation"].amplitude = fit_result["saturation_amp"] - # Update the x180 and x90 amplitudes q.xy.operations["x180"].amplitude = fit_result["x180_amp"] q.xy.operations["x90"].amplitude = fit_result["x180_amp"] / 2 @@ -247,4 +248,5 @@ def update_state(node: QualibrationNode[Parameters, Quam]): # %% {Save_results} @node.run_action() def save_results(node: QualibrationNode[Parameters, Quam]): + """Save node results.""" node.save()