From b7e2ad11ac8a367c3c18b72226dcfecd23d924b3 Mon Sep 17 00:00:00 2001 From: paulQM Date: Thu, 18 Jun 2026 17:45:45 +0200 Subject: [PATCH 01/10] feat: port jazz n and jazz 2n nodes --- .../calibration_utils/cz_jazz2_n/__init__.py | 21 ++ .../calibration_utils/cz_jazz2_n/analysis.py | 273 ++++++++++++++++ .../cz_jazz2_n/parameters.py | 50 +++ .../calibration_utils/cz_jazz2_n/plotting.py | 86 +++++ .../calibration_utils/cz_jazz_n/__init__.py | 21 ++ .../calibration_utils/cz_jazz_n/analysis.py | 309 ++++++++++++++++++ .../calibration_utils/cz_jazz_n/parameters.py | 45 +++ .../calibration_utils/cz_jazz_n/plotting.py | 89 +++++ .../32a_cz_conditional_phase.py | 8 +- .../CZ_calibrations/33b_JAZZ-N.py | 256 +++++++++++++++ .../CZ_calibrations/33c_JAZZ2-N.py | 290 ++++++++++++++++ 11 files changed, 1446 insertions(+), 2 deletions(-) create mode 100644 qualibration_graphs/superconducting/calibration_utils/cz_jazz2_n/__init__.py create mode 100644 qualibration_graphs/superconducting/calibration_utils/cz_jazz2_n/analysis.py create mode 100644 qualibration_graphs/superconducting/calibration_utils/cz_jazz2_n/parameters.py create mode 100644 qualibration_graphs/superconducting/calibration_utils/cz_jazz2_n/plotting.py create mode 100644 qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/__init__.py create mode 100644 qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/analysis.py create mode 100644 qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/parameters.py create mode 100644 qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/plotting.py create mode 100644 qualibration_graphs/superconducting/calibrations/CZ_calibrations/33b_JAZZ-N.py create mode 100644 qualibration_graphs/superconducting/calibrations/CZ_calibrations/33c_JAZZ2-N.py diff --git a/qualibration_graphs/superconducting/calibration_utils/cz_jazz2_n/__init__.py b/qualibration_graphs/superconducting/calibration_utils/cz_jazz2_n/__init__.py new file mode 100644 index 000000000..2b557dd6f --- /dev/null +++ b/qualibration_graphs/superconducting/calibration_utils/cz_jazz2_n/__init__.py @@ -0,0 +1,21 @@ +"""JAZZ2-N CZ amplitude calibration utilities.""" + +from .analysis import ( + FitResults, + coerce_to_even, + fit_raw_data, + log_fitted_results, + process_raw_dataset, +) +from .parameters import Parameters +from .plotting import plot_raw_data_with_fit + +__all__ = [ + "FitResults", + "Parameters", + "coerce_to_even", + "fit_raw_data", + "log_fitted_results", + "plot_raw_data_with_fit", + "process_raw_dataset", +] diff --git a/qualibration_graphs/superconducting/calibration_utils/cz_jazz2_n/analysis.py b/qualibration_graphs/superconducting/calibration_utils/cz_jazz2_n/analysis.py new file mode 100644 index 000000000..3ef979d44 --- /dev/null +++ b/qualibration_graphs/superconducting/calibration_utils/cz_jazz2_n/analysis.py @@ -0,0 +1,273 @@ +"""Analysis module for the JAZZ2-N CZ amplitude calibration. + +The protocol measures the joint P_|00> of the qubit pair as a function of the +CZ-pulse amplitude scale, for several repetitions N = 2k. Ignoring +decoherence, equation (40) of arXiv:2402.18926v3 predicts: + + P_|00>(amp, N) = (1 - cos((N + 1) * theta_CZ(amp))) / 2, N = 2k. + +The multipliers m = N + 1 are odd integers 1, 3, 5, ..., identical to the +JAZZ-N case. Averaging over all probed N collapses the N-dependence into +the central lobe of a sinc-like pattern centred on the optimal amplitude +(see the JAZZ-N analysis docstring for the Dirichlet derivation), and we +fit the same sinc model + + f(amp) = B + A * sinc(w * (amp - amp_0) / pi) + +with sinc(y) = sin(pi y) / (pi y) and theoretical (A, B) = (1/2, 1/2). +""" + +import logging +from dataclasses import dataclass, field +from typing import Dict, Tuple + +import numpy as np +import xarray as xr +from qualibrate import QualibrationNode +from scipy.optimize import curve_fit + + +@dataclass +class FitResults: + """JAZZ2-N fit results for a single qubit pair.""" + + optimal_amplitude: float + """Absolute optimal CZ amplitude (Volts), i.e. amp_scale_optimal * stored amplitude.""" + optimal_amplitude_scale: float + """Optimal amplitude scale factor (dimensionless; multiplied with stored amplitude).""" + success: bool + """Whether the sinc fit succeeded (with the parabolic fallback path counted as failed).""" + fit_method: str = "sinc" + """Either 'sinc' (primary) or 'parabolic' (fallback).""" + sinc_params: Dict[str, float] = field(default_factory=dict) + """Fitted sinc parameters {'A','B','amp_0','w'} when ``fit_method == 'sinc'``.""" + + +def coerce_to_even(n: int) -> int: + """Coerce an integer to the nearest even integer >= 0.""" + if n < 0: + return 0 + return 2 * int(round(n / 2.0)) + + +def log_fitted_results(fit_results: Dict[str, FitResults], log_callable=None): + """Log the JAZZ2-N fit results per qubit pair.""" + if log_callable is None: + log_callable = logging.getLogger(__name__).info + + for qp_name, fit_result in fit_results.items(): + header = f"Results for qubit pair {qp_name}: " + ("SUCCESS!\n" if fit_result.success else "FAIL!\n") + body = ( + f"\tOptimal CZ amplitude: {fit_result.optimal_amplitude:.6f} V " + f"(scale {fit_result.optimal_amplitude_scale:.6f}; method={fit_result.fit_method})" + ) + log_callable(header + body) + + +def process_raw_dataset(ds: xr.Dataset, node: QualibrationNode) -> xr.Dataset: + """Augment the raw dataset with an absolute-amplitude coordinate per qubit pair.""" + qubit_pairs = node.namespace["qubit_pairs"] + operation = node.parameters.operation + + def abs_amp(qp, amp): + return amp * qp.macros[operation].flux_pulse_qubit.amplitude + + ds = ds.assign_coords( + {"amp_full": (["qubit_pair", "amp"], np.array([abs_amp(qp, ds.amp.values) for qp in qubit_pairs]))} + ) + return ds + + +def _sinc_model(amp: np.ndarray, A: float, B: float, amp_0: float, w: float) -> np.ndarray: + """Sinc model used to localise the central peak of

_N(amp).""" + return B + A * np.sinc(w * (amp - amp_0) / np.pi) + + +def _parabolic_refine(y: np.ndarray, i: int) -> float: + """Three-point parabolic refinement around an extremum index, returning fractional index.""" + if i <= 0 or i >= len(y) - 1: + return float(i) + y1, y2, y3 = y[i - 1], y[i], y[i + 1] + denom = y1 - 2.0 * y2 + y3 + if denom == 0: + return float(i) + delta = 0.5 * (y1 - y3) / denom + return float(i) + float(np.clip(delta, -0.5, 0.5)) + + +def _argmax_with_refine(amp_values: np.ndarray, y: np.ndarray) -> float: + """Argmax of ``y`` over ``amp_values`` with 3-point parabolic refinement.""" + i = int(np.argmax(y)) + i_star = _parabolic_refine(y, i) + return float(np.interp(i_star, np.arange(len(amp_values)), amp_values)) + + +def _estimate_fwhm_around(amp_values: np.ndarray, y: np.ndarray, i_max: int) -> float: + """Estimate the full-width at half-maximum of the central peak around ``i_max``.""" + finite = y[np.isfinite(y)] + if finite.size == 0: + return float("nan") + y_max = float(y[i_max]) + y_min = float(np.min(finite)) + half = y_max - (y_max - y_min) / 2.0 + left = i_max + while left > 0 and y[left] > half: + left -= 1 + right = i_max + while right < len(y) - 1 and y[right] > half: + right += 1 + fwhm = float(amp_values[right] - amp_values[left]) + if fwhm > 0: + return fwhm + return float(amp_values[-1] - amp_values[0]) / 4.0 + + +def _fit_one_pair( + amp_values: np.ndarray, n_values: np.ndarray, p_curve: np.ndarray +) -> Tuple[float, bool, str, Dict[str, float], np.ndarray, np.ndarray]: + """Average ``p_curve`` over N and fit a sinc model to the central peak. + + Parameters + ---------- + amp_values : (n_amp,) amplitude-scale values (centred at 1.0). + n_values : (n_N,) repetition counts N (sorted ascending; not used by the + fit but retained for consistency with the earlier signature). + p_curve : (n_N, n_amp) joint P_|00> values. + + Returns + ------- + amp_seed, success, method, params, p_avg, fit_curve : see JAZZ-N analysis. + """ + if p_curve.ndim != 2 or p_curve.shape != (len(n_values), len(amp_values)): + return float("nan"), False, "none", {}, np.full_like(amp_values, np.nan), np.full_like(amp_values, np.nan) + + p_avg = np.nanmean(p_curve, axis=0) + if not np.any(np.isfinite(p_avg)): + return float("nan"), False, "none", {}, p_avg, np.full_like(amp_values, np.nan) + + p_avg_finite = np.where(np.isfinite(p_avg), p_avg, -np.inf) + i_max = int(np.argmax(p_avg_finite)) + + A_init = float(np.nanmax(p_avg) - np.nanmedian(p_avg)) + if not np.isfinite(A_init) or A_init <= 0: + A_init = 0.5 + B_init = float(np.nanmedian(p_avg)) + if not np.isfinite(B_init): + B_init = 0.5 + amp_seed_init = float(amp_values[i_max]) + fwhm = _estimate_fwhm_around(amp_values, p_avg, i_max) + if np.isfinite(fwhm) and fwhm > 0: + w_init = 3.79 / fwhm + else: + w_init = 3.79 / max(float(amp_values[-1] - amp_values[0]) / 4.0, 1e-9) + + amp_min, amp_max = float(np.min(amp_values)), float(np.max(amp_values)) + + try: + popt, _ = curve_fit( + _sinc_model, + amp_values, + p_avg, + p0=[A_init, B_init, amp_seed_init, w_init], + bounds=( + [0.0, -0.5, amp_min, w_init / 50.0], + [2.0, 1.5, amp_max, w_init * 50.0], + ), + maxfev=10000, + ) + A_fit, B_fit, amp_0_fit, w_fit = map(float, popt) + if not (amp_min <= amp_0_fit <= amp_max): + raise RuntimeError(f"Fitted amp_0 = {amp_0_fit} outside swept window.") + fit_curve = _sinc_model(amp_values, A_fit, B_fit, amp_0_fit, w_fit) + return ( + amp_0_fit, + True, + "sinc", + {"A": A_fit, "B": B_fit, "amp_0": amp_0_fit, "w": w_fit}, + p_avg, + fit_curve, + ) + except Exception: # pylint: disable=broad-except + amp_seed_refined = _argmax_with_refine(amp_values, p_avg_finite) + return ( + amp_seed_refined, + bool(np.isfinite(amp_seed_refined)), + "parabolic", + {}, + p_avg, + np.full_like(amp_values, np.nan), + ) + + +def fit_raw_data(ds: xr.Dataset, node: QualibrationNode) -> Tuple[xr.Dataset, Dict[str, FitResults]]: + """Fit the JAZZ2-N data per qubit pair, augmenting ``ds`` with the sinc-fit outputs.""" + qubit_pairs = node.namespace["qubit_pairs"] + operation = node.parameters.operation + + if "p" not in ds: + raise RuntimeError("JAZZ2-N analysis requires 'p' in the dataset (state discrimination).") + + amp_values = ds.amp.values + n_values = ds.N.values + + opt_amps_abs = [] + opt_amps_scale = [] + successes = [] + methods = [] + p_avg_rows = [] + fit_rows = [] + qp_names = ds.qubit_pair.values + fit_results: Dict[str, FitResults] = {} + + for qp_name in qp_names: + qp = next(qp for qp in qubit_pairs if qp.name == qp_name) + p00 = ds["p"].sel(qubit_pair=qp_name).transpose("N", "amp").values + amp_scale, success, method, params, p_avg, fit_curve = _fit_one_pair( + amp_values, np.asarray(n_values), np.asarray(p00) + ) + stored_amp = qp.macros[operation].flux_pulse_qubit.amplitude + amp_abs = amp_scale * stored_amp if np.isfinite(amp_scale) else np.nan + + opt_amps_abs.append(amp_abs) + opt_amps_scale.append(amp_scale) + successes.append(success) + methods.append(method) + p_avg_rows.append(p_avg) + fit_rows.append(fit_curve) + fit_results[str(qp_name)] = FitResults( + optimal_amplitude=float(amp_abs), + optimal_amplitude_scale=float(amp_scale), + success=bool(success), + fit_method=str(method), + sinc_params=params, + ) + + ds_fit = ds.assign( + { + "p_avg": xr.DataArray( + np.asarray(p_avg_rows), + dims=("qubit_pair", "amp"), + coords={"qubit_pair": qp_names, "amp": amp_values}, + name="p_avg", + ), + "sinc_fit": xr.DataArray( + np.asarray(fit_rows), + dims=("qubit_pair", "amp"), + coords={"qubit_pair": qp_names, "amp": amp_values}, + name="sinc_fit", + ), + } + ) + ds_fit = ds_fit.assign_coords( + { + "optimal_amplitude": ("qubit_pair", np.array(opt_amps_abs, dtype=float)), + "optimal_amplitude_scale": ("qubit_pair", np.array(opt_amps_scale, dtype=float)), + "success": ("qubit_pair", np.array(successes, dtype=bool)), + "fit_method": ("qubit_pair", np.array(methods, dtype=object)), + } + ) + ds_fit.optimal_amplitude.attrs = {"long_name": "optimal CZ amplitude", "units": "V"} + ds_fit.optimal_amplitude_scale.attrs = {"long_name": "optimal CZ amplitude scale", "units": "a.u."} + ds_fit.p_avg.attrs = {"long_name": ">_N", "units": "a.u."} + ds_fit.sinc_fit.attrs = {"long_name": "sinc fit", "units": "a.u."} + return ds_fit, fit_results diff --git a/qualibration_graphs/superconducting/calibration_utils/cz_jazz2_n/parameters.py b/qualibration_graphs/superconducting/calibration_utils/cz_jazz2_n/parameters.py new file mode 100644 index 000000000..827161ba2 --- /dev/null +++ b/qualibration_graphs/superconducting/calibration_utils/cz_jazz2_n/parameters.py @@ -0,0 +1,50 @@ +"""Parameters module for the JAZZ2-N CZ amplitude calibration. + +The JAZZ2-N protocol (Appendix I.1, Fig. 13(b) of arXiv:2402.18926v3) is the +two-qubit-superposition variant of JAZZ-N: both qubits receive a boundary +X_{pi/2} pulse, both are read out, and the metric is P_|00>. The repetition +N must satisfy N = 2k (k = 0, 1, 2, ...). The same X_pi refocused CZ train +gives the (2k+1)*theta_CZ phase accumulation, which is maximised when +theta_CZ = pi. Compared to JAZZ-N, the principal-peak fringe in amplitude +is roughly twice as dense for a given total pulse count, so this protocol is +intended as a finer follow-up amplitude calibration (and serves as the reward +signal for downstream Z-pulse shape optimisation). +""" + +# pylint: disable=too-few-public-methods + +from typing import ClassVar, Literal + +from qualibrate import NodeParameters +from qualibrate.core.parameters import RunnableParameters +from qualibration_libs.parameters import CommonNodeParameters, QubitPairExperimentNodeParameters + + +class NodeSpecificParameters(RunnableParameters): + """Node-specific parameters for the JAZZ2-N CZ amplitude calibration.""" + + num_shots: int = 100 + """Number of shots to average over. Default is 100.""" + amp_range: float = 0.010 + """Half-width of the amplitude-scale sweep around the stored CZ amplitude (center = 1.0). Default is 0.010.""" + amp_step: float = 0.001 + """Step of the amplitude-scale sweep. Default is 0.001.""" + N_min: int = 0 + """Minimum repetition count. Required form: N = 2k (k = 0, 1, 2, ...); auto-coerced if not. Default is 0.""" + N_max: int = 50 + """Maximum repetition count. Required form: N = 2k (k = 0, 1, 2, ...); auto-coerced if not. Default is 50.""" + operation: Literal["cz_flattop", "cz_unipolar", "cz_bipolar", "cz_flattop_erf", "cz_SNZ"] = "cz_unipolar" + """Name of the CZGate macro to drive in place of the bare Z pulse. Default is 'cz_unipolar'.""" + use_state_discrimination: bool = True + """JAZZ2-N reads the joint P_|00> of both qubits, which requires state discrimination. Setting this to False raises.""" + + +class Parameters( + NodeParameters, + CommonNodeParameters, + NodeSpecificParameters, + QubitPairExperimentNodeParameters, +): + """Combined parameters for the JAZZ2-N CZ amplitude calibration node.""" + + targets_name: ClassVar[str] = "qubit_pairs" diff --git a/qualibration_graphs/superconducting/calibration_utils/cz_jazz2_n/plotting.py b/qualibration_graphs/superconducting/calibration_utils/cz_jazz2_n/plotting.py new file mode 100644 index 000000000..9e65a2cf8 --- /dev/null +++ b/qualibration_graphs/superconducting/calibration_utils/cz_jazz2_n/plotting.py @@ -0,0 +1,86 @@ +"""Plotting module for the JAZZ2-N CZ amplitude calibration.""" + +import matplotlib.pyplot as plt +import numpy as np +import xarray as xr +from qualibration_libs.core import BatchableList + + +def plot_raw_data_with_fit(ds_fit: xr.Dataset, qubit_pairs: BatchableList) -> plt.Figure: + """Plot the JAZZ2-N data and sinc fit per qubit pair. + + For each qubit pair we show two stacked panels: + + * Top: ``p00`` as a 2D map versus the repetition count N = 2k and the + amplitude scale, with the fitted optimal amplitude drawn as a vertical + line. + * Bottom: the averaged-over-N curve ``p_avg(amp)`` together with the + fitted sinc model and the optimum, so the fit can be eyeballed. + """ + n_pairs = len(qubit_pairs) + cols = min(4, n_pairs) + rows = (n_pairs + cols - 1) // cols + fig, axes = plt.subplots(2 * rows, cols, figsize=(4.5 * cols, 5.5 * rows), squeeze=False) + + for i, qp in enumerate(qubit_pairs): + row, col = divmod(i, cols) + ax_map = axes[2 * row, col] + ax_avg = axes[2 * row + 1, col] + qp_name = qp.name + fr = ds_fit.sel(qubit_pair=qp_name) + + amps_scale = fr.amp.values + amps_abs = fr["amp_full"].values if "amp_full" in fr.coords else amps_scale + n_values = fr.N.values + + # --- Top: 2D heatmap of P_|00> --- + p_map = fr["p"].transpose("N", "amp") + xg, yg = np.meshgrid(amps_scale, n_values) + pcm = ax_map.pcolormesh(xg, yg, p_map.values, cmap="magma", shading="auto") + + opt_scale = float(fr.optimal_amplitude_scale.values) + opt_method = str(fr.fit_method.values) + if np.isfinite(opt_scale): + ax_map.axvline(opt_scale, color="lime", lw=2, label=f"opt = {opt_scale:.4f} ({opt_method})") + + def amp_scale_to_abs(s, abs_values=amps_abs, scale_values=amps_scale): + return np.interp(s, scale_values, abs_values) + + def amp_abs_to_scale(a, abs_values=amps_abs, scale_values=amps_scale): + return np.interp(a, abs_values, scale_values) + + secax = ax_map.secondary_xaxis("top", functions=(amp_scale_to_abs, amp_abs_to_scale)) + secax.set_xlabel("Amplitude (V)") + ax_map.set_title(qp_name) + ax_map.set_xlabel("Amplitude scale (a.u.)") + ax_map.set_ylabel("Repetition N = 2k") + ax_map.legend(loc="upper right", fontsize=8) + cbar = fig.colorbar(pcm, ax=ax_map, shrink=0.85) + cbar.set_label("$P_{|00\\rangle}$") + + # --- Bottom: averaged P_|00> with sinc fit --- + if "p_avg" in fr.data_vars: + ax_avg.plot(amps_scale, fr["p_avg"].values, "o", ms=3, color="C0", label=r"$\langle P_{|00\rangle}\rangle_N$") + if "sinc_fit" in fr.data_vars: + fit_vals = fr["sinc_fit"].values + if np.any(np.isfinite(fit_vals)): + ax_avg.plot(amps_scale, fit_vals, "-", lw=1.5, color="C3", label="sinc fit") + if np.isfinite(opt_scale): + ax_avg.axvline(opt_scale, color="lime", lw=1.5, label=f"opt = {opt_scale:.4f}") + ax_avg.set_xlabel("Amplitude scale (a.u.)") + ax_avg.set_ylabel(r"$\langle P_{|00\rangle}\rangle_N$") + ax_avg.legend(loc="upper right", fontsize=8) + + used = set() + for i in range(n_pairs): + row, col = divmod(i, cols) + used.add((2 * row, col)) + used.add((2 * row + 1, col)) + for r in range(axes.shape[0]): + for c in range(axes.shape[1]): + if (r, c) not in used: + axes[r, c].axis("off") + + fig.suptitle("JAZZ2-N CZ amplitude calibration") + fig.tight_layout(rect=(0, 0, 1, 0.97)) + return fig diff --git a/qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/__init__.py b/qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/__init__.py new file mode 100644 index 000000000..3ff54ea58 --- /dev/null +++ b/qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/__init__.py @@ -0,0 +1,21 @@ +"""JAZZ-N CZ amplitude calibration utilities.""" + +from .analysis import ( + FitResults, + coerce_to_4k_plus_1, + fit_raw_data, + log_fitted_results, + process_raw_dataset, +) +from .parameters import Parameters +from .plotting import plot_raw_data_with_fit + +__all__ = [ + "FitResults", + "Parameters", + "coerce_to_4k_plus_1", + "fit_raw_data", + "log_fitted_results", + "plot_raw_data_with_fit", + "process_raw_dataset", +] diff --git a/qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/analysis.py b/qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/analysis.py new file mode 100644 index 000000000..5f9a3818e --- /dev/null +++ b/qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/analysis.py @@ -0,0 +1,309 @@ +"""Analysis module for the JAZZ-N CZ amplitude calibration. + +The protocol measures the target qubit's |1> population as a function of the +CZ-pulse amplitude scale, for several echo repetitions N = 4k + 1. Ignoring +decoherence, equation (36) of arXiv:2402.18926v3 predicts: + + P_|1>(amp, N) = (1 - cos((2k+1) * theta_CZ(amp))) / 2, N = 4k + 1. + +Averaging this over all probed odd multipliers m = 2k + 1 collapses the +N-dependence into the central lobe of a sinc-like pattern centred on the +optimal amplitude. Using the Dirichlet identity + + sum_{k=k_min..k_max} cos((2k+1) x) = sin(2 K x) / (2 sin x) + +with K the number of N values, and writing x = pi + delta around the +optimum (delta ~ alpha * (amp - amp*)), one finds + + >_N(amp) = 1/2 + (1 / (4 K)) * sin(2 K delta) / sin(delta) + ~ 1/2 + (1/2) * sinc(w * (amp - amp*) / pi) (small delta) + +with sinc(y) = sin(pi y) / (pi y). The fit model used here is therefore + + f(amp) = B + A * sinc(w * (amp - amp_0) / pi) + +with theoretical (A, B) = (1/2, 1/2). The peak position amp_0 is the +optimal CZ amplitude. +""" + +import logging +from dataclasses import dataclass, field +from typing import Dict, Tuple + +import numpy as np +import xarray as xr +from qualibrate import QualibrationNode +from scipy.optimize import curve_fit + + +@dataclass +class FitResults: + """JAZZ-N fit results for a single qubit pair.""" + + optimal_amplitude: float + """Absolute optimal CZ amplitude (Volts), i.e. amp_scale_optimal * stored amplitude.""" + optimal_amplitude_scale: float + """Optimal amplitude scale factor (dimensionless; multiplied with stored amplitude).""" + success: bool + """Whether the sinc fit succeeded (with the parabolic fallback path counted as failed).""" + fit_method: str = "sinc" + """Either 'sinc' (primary) or 'parabolic' (fallback).""" + sinc_params: Dict[str, float] = field(default_factory=dict) + """Fitted sinc parameters {'A','B','amp_0','w'} when ``fit_method == 'sinc'``.""" + + +def coerce_to_4k_plus_1(n: int) -> int: + """Coerce an integer to the nearest value of the form 4k + 1 with k >= 0.""" + if n < 1: + return 1 + k = round((n - 1) / 4) + return 4 * max(0, k) + 1 + + +def log_fitted_results(fit_results: Dict[str, FitResults], log_callable=None): + """Log the JAZZ-N fit results per qubit pair.""" + if log_callable is None: + log_callable = logging.getLogger(__name__).info + + for qp_name, fit_result in fit_results.items(): + header = f"Results for qubit pair {qp_name}: " + ("SUCCESS!\n" if fit_result.success else "FAIL!\n") + body = ( + f"\tOptimal CZ amplitude: {fit_result.optimal_amplitude:.6f} V " + f"(scale {fit_result.optimal_amplitude_scale:.6f}; method={fit_result.fit_method})" + ) + log_callable(header + body) + + +def process_raw_dataset(ds: xr.Dataset, node: QualibrationNode) -> xr.Dataset: + """Augment the raw dataset with an absolute-amplitude coordinate per qubit pair.""" + qubit_pairs = node.namespace["qubit_pairs"] + operation = node.parameters.operation + + def abs_amp(qp, amp): + return amp * qp.macros[operation].flux_pulse_qubit.amplitude + + ds = ds.assign_coords( + {"amp_full": (["qubit_pair", "amp"], np.array([abs_amp(qp, ds.amp.values) for qp in qubit_pairs]))} + ) + return ds + + +def _sinc_model(amp: np.ndarray, A: float, B: float, amp_0: float, w: float) -> np.ndarray: + """Sinc model used to localise the central peak of

_N(amp). + + f(amp) = B + A * sin(w * (amp - amp_0)) / (w * (amp - amp_0)) + = B + A * np.sinc(w * (amp - amp_0) / pi). + """ + return B + A * np.sinc(w * (amp - amp_0) / np.pi) + + +def _parabolic_refine(y: np.ndarray, i: int) -> float: + """Three-point parabolic refinement around an extremum index, returning fractional index.""" + if i <= 0 or i >= len(y) - 1: + return float(i) + y1, y2, y3 = y[i - 1], y[i], y[i + 1] + denom = y1 - 2.0 * y2 + y3 + if denom == 0: + return float(i) + delta = 0.5 * (y1 - y3) / denom + return float(i) + float(np.clip(delta, -0.5, 0.5)) + + +def _argmax_with_refine(amp_values: np.ndarray, y: np.ndarray) -> float: + """Argmax of ``y`` over ``amp_values`` with 3-point parabolic refinement.""" + i = int(np.argmax(y)) + i_star = _parabolic_refine(y, i) + return float(np.interp(i_star, np.arange(len(amp_values)), amp_values)) + + +def _estimate_fwhm_around(amp_values: np.ndarray, y: np.ndarray, i_max: int) -> float: + """Estimate the full-width at half-maximum of the central peak around ``i_max``. + + Uses ``half = y_max - (y_max - y_min) / 2`` as the half-maximum reference + (i.e. half-prominence above the global minimum of ``y``). + """ + finite = y[np.isfinite(y)] + if finite.size == 0: + return float("nan") + y_max = float(y[i_max]) + y_min = float(np.min(finite)) + half = y_max - (y_max - y_min) / 2.0 + left = i_max + while left > 0 and y[left] > half: + left -= 1 + right = i_max + while right < len(y) - 1 and y[right] > half: + right += 1 + fwhm = float(amp_values[right] - amp_values[left]) + if fwhm > 0: + return fwhm + return float(amp_values[-1] - amp_values[0]) / 4.0 + + +def _fit_one_pair( + amp_values: np.ndarray, n_values: np.ndarray, p_curve: np.ndarray +) -> Tuple[float, bool, str, Dict[str, float], np.ndarray, np.ndarray]: + """Average ``p_curve`` over N and fit a sinc model to the central peak. + + Parameters + ---------- + amp_values : (n_amp,) amplitude-scale values (centred at 1.0). + n_values : (n_N,) echo counts N (sorted ascending; not used by the fit but + retained to keep the call signature consistent with the previous version). + p_curve : (n_N, n_amp) target P_|1> values. + + Returns + ------- + amp_seed : optimal amplitude scale (dimensionless). + success : True if the sinc fit converged inside the swept window. + method : 'sinc' if the primary fit succeeded, else 'parabolic'. + params : fitted sinc parameters (empty dict if parabolic fallback was used). + p_avg : (n_amp,) the averaged-over-N curve. + fit_curve : (n_amp,) the fitted-sinc evaluated at ``amp_values``; NaN if + parabolic fallback was used. + """ + if p_curve.ndim != 2 or p_curve.shape != (len(n_values), len(amp_values)): + return float("nan"), False, "none", {}, np.full_like(amp_values, np.nan), np.full_like(amp_values, np.nan) + + p_avg = np.nanmean(p_curve, axis=0) + if not np.any(np.isfinite(p_avg)): + return float("nan"), False, "none", {}, p_avg, np.full_like(amp_values, np.nan) + + p_avg_finite = np.where(np.isfinite(p_avg), p_avg, -np.inf) + i_max = int(np.argmax(p_avg_finite)) + + # Initial guesses for the sinc fit. + A_init = float(np.nanmax(p_avg) - np.nanmedian(p_avg)) + if not np.isfinite(A_init) or A_init <= 0: + A_init = 0.5 + B_init = float(np.nanmedian(p_avg)) + if not np.isfinite(B_init): + B_init = 0.5 + amp_seed_init = float(amp_values[i_max]) + fwhm = _estimate_fwhm_around(amp_values, p_avg, i_max) + # For f(amp) = B + A * sinc(w (amp - amp_0) / pi), sin(x)/x = 1/2 at x ~ 1.8955, + # so FWHM ~ 2 * 1.8955 / w ~ 3.79 / w, hence w ~ 3.79 / FWHM. + if np.isfinite(fwhm) and fwhm > 0: + w_init = 3.79 / fwhm + else: + w_init = 3.79 / max(float(amp_values[-1] - amp_values[0]) / 4.0, 1e-9) + + amp_min, amp_max = float(np.min(amp_values)), float(np.max(amp_values)) + + try: + popt, _ = curve_fit( + _sinc_model, + amp_values, + p_avg, + p0=[A_init, B_init, amp_seed_init, w_init], + bounds=( + [0.0, -0.5, amp_min, w_init / 50.0], + [2.0, 1.5, amp_max, w_init * 50.0], + ), + maxfev=10000, + ) + A_fit, B_fit, amp_0_fit, w_fit = map(float, popt) + in_range = amp_min <= amp_0_fit <= amp_max + if not in_range: + raise RuntimeError(f"Fitted amp_0 = {amp_0_fit} outside swept window.") + fit_curve = _sinc_model(amp_values, A_fit, B_fit, amp_0_fit, w_fit) + return ( + amp_0_fit, + True, + "sinc", + {"A": A_fit, "B": B_fit, "amp_0": amp_0_fit, "w": w_fit}, + p_avg, + fit_curve, + ) + except Exception: # pylint: disable=broad-except + amp_seed_refined = _argmax_with_refine(amp_values, p_avg_finite) + return ( + amp_seed_refined, + bool(np.isfinite(amp_seed_refined)), + "parabolic", + {}, + p_avg, + np.full_like(amp_values, np.nan), + ) + + +def fit_raw_data(ds: xr.Dataset, node: QualibrationNode) -> Tuple[xr.Dataset, Dict[str, FitResults]]: + """Fit the JAZZ-N data per qubit pair. + + The dataset is augmented with three new data variables (per qubit pair, on + the amp axis): ``p_avg`` (mean over N of the target |1> population) and + ``sinc_fit`` (the fitted sinc model evaluated on the amp grid), plus + ``optimal_amplitude`` / ``optimal_amplitude_scale`` / ``success`` / + ``fit_method`` as coordinates. + """ + qubit_pairs = node.namespace["qubit_pairs"] + operation = node.parameters.operation + + data_var = "state_target" if "state_target" in ds else None + if data_var is None: + raise RuntimeError("JAZZ-N analysis requires 'state_target' in the dataset (state discrimination).") + + amp_values = ds.amp.values + n_values = ds.N.values + + opt_amps_abs = [] + opt_amps_scale = [] + successes = [] + methods = [] + p_avg_rows = [] + fit_rows = [] + qp_names = ds.qubit_pair.values + fit_results: Dict[str, FitResults] = {} + + for qp_name in qp_names: + qp = next(qp for qp in qubit_pairs if qp.name == qp_name) + p = ds[data_var].sel(qubit_pair=qp_name).transpose("N", "amp").values + amp_scale, success, method, params, p_avg, fit_curve = _fit_one_pair( + amp_values, np.asarray(n_values), np.asarray(p) + ) + stored_amp = qp.macros[operation].flux_pulse_qubit.amplitude + amp_abs = amp_scale * stored_amp if np.isfinite(amp_scale) else np.nan + + opt_amps_abs.append(amp_abs) + opt_amps_scale.append(amp_scale) + successes.append(success) + methods.append(method) + p_avg_rows.append(p_avg) + fit_rows.append(fit_curve) + fit_results[str(qp_name)] = FitResults( + optimal_amplitude=float(amp_abs), + optimal_amplitude_scale=float(amp_scale), + success=bool(success), + fit_method=str(method), + sinc_params=params, + ) + + ds_fit = ds.assign( + { + "p_avg": xr.DataArray( + np.asarray(p_avg_rows), + dims=("qubit_pair", "amp"), + coords={"qubit_pair": qp_names, "amp": amp_values}, + name="p_avg", + ), + "sinc_fit": xr.DataArray( + np.asarray(fit_rows), + dims=("qubit_pair", "amp"), + coords={"qubit_pair": qp_names, "amp": amp_values}, + name="sinc_fit", + ), + } + ) + ds_fit = ds_fit.assign_coords( + { + "optimal_amplitude": ("qubit_pair", np.array(opt_amps_abs, dtype=float)), + "optimal_amplitude_scale": ("qubit_pair", np.array(opt_amps_scale, dtype=float)), + "success": ("qubit_pair", np.array(successes, dtype=bool)), + "fit_method": ("qubit_pair", np.array(methods, dtype=object)), + } + ) + ds_fit.optimal_amplitude.attrs = {"long_name": "optimal CZ amplitude", "units": "V"} + ds_fit.optimal_amplitude_scale.attrs = {"long_name": "optimal CZ amplitude scale", "units": "a.u."} + ds_fit.p_avg.attrs = {"long_name": ">_N", "units": "a.u."} + ds_fit.sinc_fit.attrs = {"long_name": "sinc fit", "units": "a.u."} + return ds_fit, fit_results diff --git a/qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/parameters.py b/qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/parameters.py new file mode 100644 index 000000000..0f03e8be1 --- /dev/null +++ b/qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/parameters.py @@ -0,0 +1,45 @@ +"""Parameters module for the JAZZ-N CZ amplitude calibration. + +The JAZZ-N protocol (Appendix I.1, Fig. 13(a) of arXiv:2402.18926v3) measures +P_|1> of the target qubit after a refocused train of CZ gates. The number of +X_pi echo pulses N must satisfy N = 4k + 1 (k = 0, 1, 2, ...), which gives a +clean (2k+1)*theta_CZ phase accumulation that peaks at theta_CZ = pi. +""" + +# pylint: disable=too-few-public-methods + +from typing import ClassVar, Literal + +from qualibrate import NodeParameters +from qualibrate.core.parameters import RunnableParameters +from qualibration_libs.parameters import CommonNodeParameters, QubitPairExperimentNodeParameters + + +class NodeSpecificParameters(RunnableParameters): + """Node-specific parameters for the JAZZ-N CZ amplitude calibration.""" + + num_shots: int = 100 + """Number of shots to average over. Default is 100.""" + amp_range: float = 0.010 + """Half-width of the amplitude-scale sweep around the stored CZ amplitude (center = 1.0). Default is 0.010.""" + amp_step: float = 0.001 + """Step of the amplitude-scale sweep. Default is 0.001.""" + N_min: int = 1 + """Minimum number of X_pi echo pulses. Required form: N = 4k + 1; auto-coerced if not. Default is 1.""" + N_max: int = 101 + """Maximum number of X_pi echo pulses. Required form: N = 4k + 1; auto-coerced if not. Default is 101.""" + operation: Literal["cz_flattop", "cz_unipolar", "cz_bipolar", "cz_flattop_erf", "cz_SNZ"] = "cz_unipolar" + """Name of the CZGate macro to drive in place of the bare Z pulse. Default is 'cz_unipolar'.""" + use_state_discrimination: bool = True + """JAZZ-N reads P_|1> of the target qubit, which requires state discrimination. Setting this to False raises.""" + + +class Parameters( + NodeParameters, + CommonNodeParameters, + NodeSpecificParameters, + QubitPairExperimentNodeParameters, +): + """Combined parameters for the JAZZ-N CZ amplitude calibration node.""" + + targets_name: ClassVar[str] = "qubit_pairs" diff --git a/qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/plotting.py b/qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/plotting.py new file mode 100644 index 000000000..dea67bb8a --- /dev/null +++ b/qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/plotting.py @@ -0,0 +1,89 @@ +"""Plotting module for the JAZZ-N CZ amplitude calibration.""" + +import matplotlib.pyplot as plt +import numpy as np +import xarray as xr +from qualibration_libs.core import BatchableList + + +def plot_raw_data_with_fit(ds_fit: xr.Dataset, qubit_pairs: BatchableList) -> plt.Figure: + """Plot the JAZZ-N data and sinc fit per qubit pair. + + For each qubit pair we show two stacked panels: + + * Top: ``state_target`` (target |1> population) as a 2D map versus the + echo count N and the amplitude scale, with the fitted optimal amplitude + drawn as a vertical line. + * Bottom: the averaged-over-N curve ``p_avg(amp)`` together with the + fitted sinc model and the optimum, so the fit can be eyeballed. + """ + n_pairs = len(qubit_pairs) + cols = min(4, n_pairs) + rows = (n_pairs + cols - 1) // cols + fig, axes = plt.subplots(2 * rows, cols, figsize=(4.5 * cols, 5.5 * rows), squeeze=False) + + for i, qp in enumerate(qubit_pairs): + row, col = divmod(i, cols) + ax_map = axes[2 * row, col] + ax_avg = axes[2 * row + 1, col] + qp_name = qp.name + fr = ds_fit.sel(qubit_pair=qp_name) + + amps_scale = fr.amp.values + amps_abs = fr["amp_full"].values if "amp_full" in fr.coords else amps_scale + n_values = fr.N.values + + # --- Top: 2D heatmap of P_|1> --- + p_map = fr["state_target"].transpose("N", "amp") + xg, yg = np.meshgrid(amps_scale, n_values) + pcm = ax_map.pcolormesh(xg, yg, p_map.values, cmap="magma", shading="auto") + + opt_scale = float(fr.optimal_amplitude_scale.values) + opt_method = str(fr.fit_method.values) + if np.isfinite(opt_scale): + ax_map.axvline(opt_scale, color="lime", lw=2, label=f"opt = {opt_scale:.4f} ({opt_method})") + + def amp_scale_to_abs(s, abs_values=amps_abs, scale_values=amps_scale): + return np.interp(s, scale_values, abs_values) + + def amp_abs_to_scale(a, abs_values=amps_abs, scale_values=amps_scale): + return np.interp(a, abs_values, scale_values) + + secax = ax_map.secondary_xaxis("top", functions=(amp_scale_to_abs, amp_abs_to_scale)) + secax.set_xlabel("Amplitude (V)") + ax_map.set_title(qp_name) + ax_map.set_xlabel("Amplitude scale (a.u.)") + ax_map.set_ylabel("Echo count N = 4k + 1") + ax_map.legend(loc="upper right", fontsize=8) + cbar = fig.colorbar(pcm, ax=ax_map, shrink=0.85) + cbar.set_label("$P_{|1\\rangle}$ of target") + + # --- Bottom: averaged P_|1> with sinc fit --- + if "p_avg" in fr.data_vars: + ax_avg.plot(amps_scale, fr["p_avg"].values, "o", ms=3, color="C0", label=r"$\langle P_{|1\rangle}\rangle_N$") + if "sinc_fit" in fr.data_vars: + fit_vals = fr["sinc_fit"].values + if np.any(np.isfinite(fit_vals)): + ax_avg.plot(amps_scale, fit_vals, "-", lw=1.5, color="C3", label="sinc fit") + if np.isfinite(opt_scale): + ax_avg.axvline(opt_scale, color="lime", lw=1.5, label=f"opt = {opt_scale:.4f}") + ax_avg.set_xlabel("Amplitude scale (a.u.)") + ax_avg.set_ylabel(r"$\langle P_{|1\rangle}\rangle_N$") + ax_avg.legend(loc="upper right", fontsize=8) + + # Hide unused panels. + total_axes = axes.flatten() + used = set() + for i in range(n_pairs): + row, col = divmod(i, cols) + used.add((2 * row, col)) + used.add((2 * row + 1, col)) + for r in range(axes.shape[0]): + for c in range(axes.shape[1]): + if (r, c) not in used: + axes[r, c].axis("off") + del total_axes + + fig.suptitle("JAZZ-N CZ amplitude calibration") + fig.tight_layout(rect=(0, 0, 1, 0.97)) + return fig diff --git a/qualibration_graphs/superconducting/calibrations/CZ_calibrations/32a_cz_conditional_phase.py b/qualibration_graphs/superconducting/calibrations/CZ_calibrations/32a_cz_conditional_phase.py index e7916ee15..46e0b96da 100644 --- a/qualibration_graphs/superconducting/calibrations/CZ_calibrations/32a_cz_conditional_phase.py +++ b/qualibration_graphs/superconducting/calibrations/CZ_calibrations/32a_cz_conditional_phase.py @@ -73,6 +73,10 @@ @node.run_action(skip_if=node.modes.external) def custom_param(node: QualibrationNode[Parameters, Quam]): # You can get type hinting in your IDE by typing node.parameters. + node.parameters.amp_range = 0.06 + node.parameters.reset_type = "active" + node.parameters.operation = "cz_unipolar" + node.parameters.qubit_pairs = ["coupler_q4_q5"] pass @@ -175,8 +179,8 @@ def create_qua_program(node: QualibrationNode[Parameters, Quam]): if node.parameters.use_state_discrimination: # measure g/e/f populations for both qubits - mq.readout_state_gef(state_mq[ii]) - sq.readout_state_gef(state_sq[ii]) + mq.readout_state_gef(state_mq[ii], pulse_name="readout") + sq.readout_state_gef(state_sq[ii], pulse_name="readout") with switch_(state_mq[ii]): with case_(0): wait(4) diff --git a/qualibration_graphs/superconducting/calibrations/CZ_calibrations/33b_JAZZ-N.py b/qualibration_graphs/superconducting/calibrations/CZ_calibrations/33b_JAZZ-N.py new file mode 100644 index 000000000..0cd52355d --- /dev/null +++ b/qualibration_graphs/superconducting/calibrations/CZ_calibrations/33b_JAZZ-N.py @@ -0,0 +1,256 @@ +# pylint: disable=R0801 +# pylint: disable=duplicate-code + +# %% {Imports} +from dataclasses import asdict + +import matplotlib.pyplot as plt +import numpy as np +import xarray as xr +from calibration_utils.cz_jazz_n import ( + Parameters, + coerce_to_4k_plus_1, + fit_raw_data, + log_fitted_results, + plot_raw_data_with_fit, + process_raw_dataset, +) +from qm.qua import * +from qualang_tools.loops import from_array +from qualang_tools.multi_user import qm_session +from qualang_tools.results import progress_counter +from qualang_tools.units import unit +from qualibrate import QualibrationNode +from qualibration_libs.data import XarrayDataFetcher +from qualibration_libs.parameters import get_qubit_pairs +from qualibration_libs.runtime import simulate_and_plot +from quam_config import Quam + + +# %% {Initialisation} +description = """ + JAZZ-N CZ AMPLITUDE CALIBRATION +This node calibrates the CZ-pulse amplitude using the JAZZ-N protocol +(arXiv:2402.18926v3, Appendix I.1, Fig. 13(a)). The pulse sequence is + + x90(target) -- CZ -- [X_pi(control) & X_pi(target) -- CZ] x N -- x90(target) -- measure(target) + +where N = 4k + 1 (k = 0, 1, 2, ...). With the X_pi refocusing pulses on both +qubits, the target |1> population evolves as + + P_|1>(target) = (1 - cos((2k+1) * theta_CZ)) / 2, + +independently of any virtual-Z (single-qubit) phase shifts inside the CZ +macro. The optimal CZ amplitude is the value where theta_CZ = pi, i.e. where +P_|1> is maximal. As N grows the peak around theta_CZ = pi becomes sharper, +giving finer amplitude resolution. + +The Z-pulse is supplied by the full CZGate macro selected via the +``operation`` parameter, so any qubit/coupler flux pulse shape that the +macro provides (cz_unipolar, cz_flattop, cz_bipolar, cz_flattop_erf, cz_SNZ) +can be calibrated. + +Prerequisites: + - Calibrated single-qubit gates (x90, x180) for both qubits in the pair. + - Calibrated, state-discriminating readout for the target qubit. + - An initial estimate of the CZ amplitude (e.g. from 32a_cz_conditional_phase + or 32b_cz_conditional_phase_error_amp). + +State update: + - qubit_pair.macros[operation].flux_pulse_qubit.amplitude (fitted optimal CZ amplitude). +""" + +node = QualibrationNode[Parameters, Quam]( + name="33b_JAZZ_N", + 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 local debug parameter overrides when running directly from IDE.""" + # node.parameters.qubit_pairs = ["q1-q2"] + pass + + +# %% {Create_QUA_program} +@node.run_action(skip_if=node.parameters.load_data_id is not None) +def create_qua_program(node: QualibrationNode[Parameters, Quam]): # pylint: disable=too-many-statements + """Create the sweep axes and generate the QUA program for the JAZZ-N sequence.""" + unit(coerce_to_integer=True) + + if not node.parameters.use_state_discrimination: + raise RuntimeError( + "JAZZ-N reads the target qubit |1> population and therefore requires " + "use_state_discrimination = True." + ) + + node.namespace["qubit_pairs"] = qubit_pairs = get_qubit_pairs(node) + num_qubit_pairs = len(qubit_pairs) + + # Coerce N_min / N_max to the nearest 4k + 1 (>= 1) and warn if changed. + n_min_req = int(node.parameters.N_min) + n_max_req = int(node.parameters.N_max) + n_min = coerce_to_4k_plus_1(n_min_req) + n_max = coerce_to_4k_plus_1(n_max_req) + if n_min > n_max: + n_min, n_max = n_max, n_min + if n_min != n_min_req: + node.log(f"N_min {n_min_req} coerced to nearest 4k+1 value: {n_min}.") + if n_max != n_max_req: + node.log(f"N_max {n_max_req} coerced to nearest 4k+1 value: {n_max}.") + + n_avg = node.parameters.num_shots + amplitudes = np.arange(1 - node.parameters.amp_range, 1 + node.parameters.amp_range, node.parameters.amp_step) + n_values = np.arange(n_min, n_max + 1, 4, dtype=int) + operation = node.parameters.operation + + node.namespace["sweep_axes"] = { + "qubit_pair": xr.DataArray(qubit_pairs.get_names()), + "N": xr.DataArray(n_values, attrs={"long_name": "echo count N = 4k+1"}), + "amp": xr.DataArray(amplitudes, attrs={"long_name": "amplitude scale", "units": "a.u."}), + } + + with program() as node.namespace["qua_program"]: + amp = declare(fixed) + n = declare(int) + n_op = declare(int) + count = declare(int) + n_st = declare_stream() + state_t = [declare(int) for _ in range(num_qubit_pairs)] + state_t_st = [declare_stream() for _ in range(num_qubit_pairs)] + + for multiplexed_qubit_pairs in qubit_pairs.batch(): + for qp in multiplexed_qubit_pairs.values(): + node.machine.initialize_qpu(target=qp.qubit_control) + node.machine.initialize_qpu(target=qp.qubit_target) + + with for_(n, 0, n < n_avg, n + 1): + save(n, n_st) + with for_(n_op, int(n_min), n_op <= int(n_max), n_op + 4): + with for_(*from_array(amp, amplitudes)): + for ii, qp in multiplexed_qubit_pairs.items(): + qp.qubit_control.reset(node.parameters.reset_type, node.parameters.simulate) + qp.qubit_target.reset(node.parameters.reset_type, node.parameters.simulate) + qp.align() + reset_frame(qp.qubit_target.xy.name) + reset_frame(qp.qubit_control.xy.name) + + # Initial pi/2 on target (Q2 role); control (Q1 role) stays in |0>. + qp.qubit_target.xy.play("x90") + qp.align() + + # First CZ (the "Z" preceding the (pi-Z)^N pattern). + qp.macros[operation].apply(amplitude_scale_qubit=amp) + + # N echoes interleaved with N more CZs: [X_pi X_pi, CZ] x N. + with for_(count, 0, count < n_op, count + 1): + qp.align() + qp.qubit_control.xy.play("x180") + qp.qubit_target.xy.play("x180") + qp.align() + qp.macros[operation].apply(amplitude_scale_qubit=amp) + + qp.align() + qp.qubit_target.xy.play("x90") + qp.align() + + qp.qubit_target.readout_state(state_t[ii]) + save(state_t[ii], state_t_st[ii]) + + align() + + with stream_processing(): + n_st.save("n") + for ii in range(num_qubit_pairs): + state_t_st[ii].buffer(len(amplitudes)).buffer(len(n_values)).average().save(f"state_target{ii + 1}") + + +# %% {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.""" + qmm = node.machine.connect() + config = node.machine.generate_config() + samples, fig, wf_report = simulate_and_plot(qmm, config, node.namespace["qua_program"], node.parameters) + node.results["simulation"] = {"figure": fig, "wf_report": wf_report.to_dict(), "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.""" + qmm = node.machine.connect() + config = node.machine.generate_config() + with qm_session(qmm, config, timeout=node.parameters.timeout) as qm: + node.namespace["job"] = job = qm.execute(node.namespace["qua_program"]) + data_fetcher = XarrayDataFetcher(job, node.namespace["sweep_axes"]) + for dataset in data_fetcher: + progress_counter( + data_fetcher.get("n", 0), + node.parameters.num_shots, + start_time=data_fetcher.t_start, + ) + node.log(job.execution_report()) + node.results["ds_raw"] = dataset + + +# %% {Load_data} +@node.run_action(skip_if=node.parameters.load_data_id is None) +def load_data(node: QualibrationNode[Parameters, Quam]): + """Load a previously acquired dataset.""" + load_data_id = node.parameters.load_data_id + node.load_from_id(node.parameters.load_data_id) + node.parameters.load_data_id = load_data_id + node.namespace["qubit_pairs"] = get_qubit_pairs(node) + + +# %% {Analyse_data} +@node.run_action(skip_if=node.parameters.simulate) +def analyse_data(node: QualibrationNode[Parameters, Quam]): + """Process the raw data, run the coarse-to-fine fit and set node outcomes.""" + 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_fitted_results(fit_results, log_callable=node.log) + node.outcomes = { + qp_name: ("successful" if fit_result.success else "failed") for qp_name, fit_result in fit_results.items() + } + + +# %% {Plot_data} +@node.run_action(skip_if=node.parameters.simulate) +def plot_data(node: QualibrationNode[Parameters, Quam]): + """Plot the JAZZ-N P_|1> map and the fitted optimal amplitude for each qubit pair.""" + fig = plot_raw_data_with_fit(node.results["ds_fit"], node.namespace["qubit_pairs"]) + plt.show() + node.results["figures"] = {"jazz_n_amplitude": fig} + + +# %% {Update_state} +@node.run_action(skip_if=node.parameters.simulate) +def update_state(node: QualibrationNode[Parameters, Quam]): + """Update the CZ flux-pulse amplitude for every successfully fitted qubit pair.""" + operation = node.parameters.operation + fit_results = node.results["fit_results"] + with node.record_state_updates(): + for qp in node.namespace["qubit_pairs"]: + if node.outcomes[qp.name] == "failed": + node.log(f"Skipping state update for {qp.name}: fit failed.") + continue + qp.macros[operation].flux_pulse_qubit.amplitude = fit_results[qp.name]["optimal_amplitude"] + + +# %% {Save_results} +@node.run_action() +def save_results(node: QualibrationNode[Parameters, Quam]): + """Save the calibration results.""" + node.save() + + +# %% diff --git a/qualibration_graphs/superconducting/calibrations/CZ_calibrations/33c_JAZZ2-N.py b/qualibration_graphs/superconducting/calibrations/CZ_calibrations/33c_JAZZ2-N.py new file mode 100644 index 000000000..4b507a299 --- /dev/null +++ b/qualibration_graphs/superconducting/calibrations/CZ_calibrations/33c_JAZZ2-N.py @@ -0,0 +1,290 @@ +# pylint: disable=R0801 +# pylint: disable=duplicate-code + +# %% {Imports} +from dataclasses import asdict + +import matplotlib.pyplot as plt +import numpy as np +import xarray as xr +from calibration_utils.cz_jazz2_n import ( + Parameters, + coerce_to_even, + fit_raw_data, + log_fitted_results, + plot_raw_data_with_fit, + process_raw_dataset, +) +from qm.qua import * +from qualang_tools.loops import from_array +from qualang_tools.multi_user import qm_session +from qualang_tools.results import progress_counter +from qualang_tools.units import unit +from qualibrate import QualibrationNode +from qualibration_libs.data import XarrayDataFetcher +from qualibration_libs.parameters import get_qubit_pairs +from qualibration_libs.runtime import simulate_and_plot +from quam_config import Quam + + +# %% {Initialisation} +description = """ + JAZZ2-N CZ AMPLITUDE CALIBRATION +This node calibrates the CZ-pulse amplitude using the JAZZ2-N protocol +(arXiv:2402.18926v3, Appendix I.1, Fig. 13(b)). The pulse sequence is + + x90(control) & x90(target) (X_{pi/2} X_{pi/2}) + CZ (initial Z) + [X_pi(control) & X_pi(target) -- CZ] x (2N + 1) + x90(control) & x90(target) (X_{pi/2} X_{pi/2}) + measure(control), measure(target) -> p00 = (1 - state_c) * (1 - state_t) + +where N = 2k (k = 0, 1, 2, ...). With the X_pi refocusing pulses on both +qubits, the joint ground-state probability evolves as + + P_|00>(amp, N) = (1 - cos((N + 1) * theta_CZ(amp))) / 2, + +independently of any virtual-Z (single-qubit) phase shifts inside the CZ +macro. The optimal CZ amplitude is the value where theta_CZ = pi, i.e. where +P_|00> is maximal. Compared to JAZZ-N, the principal-peak fringe is denser in +amplitude for a given total pulse count, so this node is a sharper follow-up +amplitude calibration; the same measurement can also be used downstream as +the reward signal for Z-pulse shape optimisation. + +The Z-pulse is supplied by the full CZGate macro selected via the +``operation`` parameter, so any qubit/coupler flux pulse shape that the +macro provides (cz_unipolar, cz_flattop, cz_bipolar, cz_flattop_erf, cz_SNZ) +can be calibrated. + +Prerequisites: + - Calibrated single-qubit gates (x90, x180) for both qubits in the pair. + - Calibrated, state-discriminating readout for BOTH qubits. + - An initial estimate of the CZ amplitude (e.g. from 33b_JAZZ-N). + +State update: + - qubit_pair.macros[operation].flux_pulse_qubit.amplitude (fitted optimal CZ amplitude). +""" + +node = QualibrationNode[Parameters, Quam]( + name="33c_JAZZ2_N", + 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 local debug parameter overrides when running directly from IDE.""" + # node.parameters.qubit_pairs = ["q1-q2"] + pass + + +# %% {Create_QUA_program} +@node.run_action(skip_if=node.parameters.load_data_id is not None) +def create_qua_program(node: QualibrationNode[Parameters, Quam]): # pylint: disable=too-many-statements + """Create the sweep axes and generate the QUA program for the JAZZ2-N sequence.""" + unit(coerce_to_integer=True) + + if not node.parameters.use_state_discrimination: + raise RuntimeError( + "JAZZ2-N reads the joint P_|00> of the qubit pair and therefore requires " + "use_state_discrimination = True." + ) + + node.namespace["qubit_pairs"] = qubit_pairs = get_qubit_pairs(node) + num_qubit_pairs = len(qubit_pairs) + + # Coerce N_min / N_max to the nearest even integer (>= 0) and warn if changed. + n_min_req = int(node.parameters.N_min) + n_max_req = int(node.parameters.N_max) + n_min = coerce_to_even(n_min_req) + n_max = coerce_to_even(n_max_req) + if n_min > n_max: + n_min, n_max = n_max, n_min + if n_min != n_min_req: + node.log(f"N_min {n_min_req} coerced to nearest even value: {n_min}.") + if n_max != n_max_req: + node.log(f"N_max {n_max_req} coerced to nearest even value: {n_max}.") + + n_avg = node.parameters.num_shots + amplitudes = np.arange(1 - node.parameters.amp_range, 1 + node.parameters.amp_range, node.parameters.amp_step) + n_values = np.arange(n_min, n_max + 1, 2, dtype=int) + operation = node.parameters.operation + + node.namespace["sweep_axes"] = { + "qubit_pair": xr.DataArray(qubit_pairs.get_names()), + "N": xr.DataArray(n_values, attrs={"long_name": "repetition count N = 2k"}), + "amp": xr.DataArray(amplitudes, attrs={"long_name": "amplitude scale", "units": "a.u."}), + } + + with program() as node.namespace["qua_program"]: + amp = declare(fixed) + n = declare(int) + n_op = declare(int) + count = declare(int) + n_st = declare_stream() + state_c = [declare(int) for _ in range(num_qubit_pairs)] + state_t = [declare(int) for _ in range(num_qubit_pairs)] + p00 = [declare(int) for _ in range(num_qubit_pairs)] + p00_st = [declare_stream() for _ in range(num_qubit_pairs)] + state_c_st = [declare_stream() for _ in range(num_qubit_pairs)] + state_t_st = [declare_stream() for _ in range(num_qubit_pairs)] + + for multiplexed_qubit_pairs in qubit_pairs.batch(): + for qp in multiplexed_qubit_pairs.values(): + node.machine.initialize_qpu(target=qp.qubit_control) + node.machine.initialize_qpu(target=qp.qubit_target) + + with for_(n, 0, n < n_avg, n + 1): + save(n, n_st) + with for_(n_op, n_min, n_op <= n_max, n_op + 2): + with for_(*from_array(amp, amplitudes)): + for ii, qp in multiplexed_qubit_pairs.items(): + qp.qubit_control.reset(node.parameters.reset_type, node.parameters.simulate) + qp.qubit_target.reset(node.parameters.reset_type, node.parameters.simulate) + qp.align() + reset_frame(qp.qubit_target.xy.name) + reset_frame(qp.qubit_control.xy.name) + + # Boundary X_{pi/2} X_{pi/2} (both qubits). + qp.qubit_control.xy.play("x90") + qp.qubit_target.xy.play("x90") + qp.align() + + # First CZ (the "Z" preceding the (pi-Z)^(2N+1) pattern). + qp.macros[operation].apply(amplitude_scale_qubit=amp) + + qp.qubit_control.xy.play("x180") + qp.qubit_target.xy.play("x180") + qp.align() + + # First CZ (the "Z" preceding the (pi-Z)^(2N+1) pattern). + qp.macros[operation].apply(amplitude_scale_qubit=amp) + + # (X_pi X_pi, CZ) x (2N + 1). + with for_(count, 1, count <= n_op, count + 1): + qp.qubit_control.xy.play("x180") + qp.qubit_target.xy.play("x180") + qp.align() + qp.macros[operation].apply(amplitude_scale_qubit=amp) + qp.qubit_control.xy.frame_rotation_2pi(0.5) + qp.qubit_target.xy.frame_rotation_2pi(0.5) + qp.qubit_control.xy.play("x180") + qp.qubit_target.xy.play("x180") + qp.align() + qp.macros[operation].apply(amplitude_scale_qubit=amp) + qp.qubit_control.xy.frame_rotation_2pi(-0.5) + qp.qubit_target.xy.frame_rotation_2pi(-0.5) + + qp.align() + # Boundary X_{pi/2} X_{pi/2} (both qubits). + qp.qubit_control.xy.play("x90") + qp.qubit_control.xy.play("x180") + qp.qubit_target.xy.play("x90") + qp.qubit_target.xy.play("x180") + qp.align() + + qp.qubit_control.readout_state(state_c[ii]) + qp.qubit_target.readout_state(state_t[ii]) + assign(p00[ii], (1 - state_c[ii]) * (1 - state_t[ii])) + save(p00[ii], p00_st[ii]) + save(state_c[ii], state_c_st[ii]) + save(state_t[ii], state_t_st[ii]) + + align() + + with stream_processing(): + n_st.save("n") + for ii in range(num_qubit_pairs): + p00_st[ii].buffer(len(amplitudes)).buffer(len(n_values)).average().save(f"p{ii + 1}") + state_c_st[ii].buffer(len(amplitudes)).buffer(len(n_values)).average().save(f"state_c{ii + 1}") + state_t_st[ii].buffer(len(amplitudes)).buffer(len(n_values)).average().save(f"state_t{ii + 1}") + + +# %% {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.""" + qmm = node.machine.connect() + config = node.machine.generate_config() + samples, fig, wf_report = simulate_and_plot(qmm, config, node.namespace["qua_program"], node.parameters) + node.results["simulation"] = {"figure": fig, "wf_report": wf_report.to_dict(), "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.""" + qmm = node.machine.connect() + config = node.machine.generate_config() + with qm_session(qmm, config, timeout=node.parameters.timeout) as qm: + node.namespace["job"] = job = qm.execute(node.namespace["qua_program"]) + data_fetcher = XarrayDataFetcher(job, node.namespace["sweep_axes"]) + for dataset in data_fetcher: + progress_counter( + data_fetcher.get("n", 0), + node.parameters.num_shots, + start_time=data_fetcher.t_start, + ) + node.log(job.execution_report()) + node.results["ds_raw"] = dataset + + +# %% {Load_data} +@node.run_action(skip_if=node.parameters.load_data_id is None) +def load_data(node: QualibrationNode[Parameters, Quam]): + """Load a previously acquired dataset.""" + load_data_id = node.parameters.load_data_id + node.load_from_id(node.parameters.load_data_id) + node.parameters.load_data_id = load_data_id + node.namespace["qubit_pairs"] = get_qubit_pairs(node) + + +# %% {Analyse_data} +@node.run_action(skip_if=node.parameters.simulate) +def analyse_data(node: QualibrationNode[Parameters, Quam]): + """Process the raw data, run the coarse-to-fine fit and set node outcomes.""" + 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_fitted_results(fit_results, log_callable=node.log) + node.outcomes = { + qp_name: ("successful" if fit_result.success else "failed") for qp_name, fit_result in fit_results.items() + } + + +# %% {Plot_data} +@node.run_action(skip_if=node.parameters.simulate) +def plot_data(node: QualibrationNode[Parameters, Quam]): + """Plot the JAZZ2-N P_|00> map and the fitted optimal amplitude for each qubit pair.""" + fig = plot_raw_data_with_fit(node.results["ds_fit"], node.namespace["qubit_pairs"]) + plt.show() + node.results["figures"] = {"jazz2_n_amplitude": fig} + + +# %% {Update_state} +@node.run_action(skip_if=node.parameters.simulate) +def update_state(node: QualibrationNode[Parameters, Quam]): + """Update the CZ flux-pulse amplitude for every successfully fitted qubit pair.""" + operation = node.parameters.operation + fit_results = node.results["fit_results"] + with node.record_state_updates(): + for qp in node.namespace["qubit_pairs"]: + if node.outcomes[qp.name] == "failed": + node.log(f"Skipping state update for {qp.name}: fit failed.") + continue + qp.macros[operation].flux_pulse_qubit.amplitude = fit_results[qp.name]["optimal_amplitude"] + + +# %% {Save_results} +@node.run_action() +def save_results(node: QualibrationNode[Parameters, Quam]): + """Save the calibration results.""" + node.save() + + +# %% From 8926ebaf0b8c8d600f40e4cf8995ebf79b412a71 Mon Sep 17 00:00:00 2001 From: paulQM Date: Thu, 18 Jun 2026 17:57:15 +0200 Subject: [PATCH 02/10] tested and new stream declaration --- .../calibrations/CZ_calibrations/33b_JAZZ-N.py | 4 ++-- .../calibrations/CZ_calibrations/33c_JAZZ2-N.py | 8 ++++---- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/qualibration_graphs/superconducting/calibrations/CZ_calibrations/33b_JAZZ-N.py b/qualibration_graphs/superconducting/calibrations/CZ_calibrations/33b_JAZZ-N.py index 0cd52355d..686f7fad5 100644 --- a/qualibration_graphs/superconducting/calibrations/CZ_calibrations/33b_JAZZ-N.py +++ b/qualibration_graphs/superconducting/calibrations/CZ_calibrations/33b_JAZZ-N.py @@ -120,9 +120,9 @@ def create_qua_program(node: QualibrationNode[Parameters, Quam]): # pylint: dis n = declare(int) n_op = declare(int) count = declare(int) - n_st = declare_stream() + n_st = declare_output_stream() state_t = [declare(int) for _ in range(num_qubit_pairs)] - state_t_st = [declare_stream() for _ in range(num_qubit_pairs)] + state_t_st = [declare_output_stream() for _ in range(num_qubit_pairs)] for multiplexed_qubit_pairs in qubit_pairs.batch(): for qp in multiplexed_qubit_pairs.values(): diff --git a/qualibration_graphs/superconducting/calibrations/CZ_calibrations/33c_JAZZ2-N.py b/qualibration_graphs/superconducting/calibrations/CZ_calibrations/33c_JAZZ2-N.py index 4b507a299..ef7b056e6 100644 --- a/qualibration_graphs/superconducting/calibrations/CZ_calibrations/33c_JAZZ2-N.py +++ b/qualibration_graphs/superconducting/calibrations/CZ_calibrations/33c_JAZZ2-N.py @@ -125,13 +125,13 @@ def create_qua_program(node: QualibrationNode[Parameters, Quam]): # pylint: dis n = declare(int) n_op = declare(int) count = declare(int) - n_st = declare_stream() + n_st = declare_output_stream() state_c = [declare(int) for _ in range(num_qubit_pairs)] state_t = [declare(int) for _ in range(num_qubit_pairs)] p00 = [declare(int) for _ in range(num_qubit_pairs)] - p00_st = [declare_stream() for _ in range(num_qubit_pairs)] - state_c_st = [declare_stream() for _ in range(num_qubit_pairs)] - state_t_st = [declare_stream() for _ in range(num_qubit_pairs)] + p00_st = [declare_output_stream() for _ in range(num_qubit_pairs)] + state_c_st = [declare_output_stream() for _ in range(num_qubit_pairs)] + state_t_st = [declare_output_stream() for _ in range(num_qubit_pairs)] for multiplexed_qubit_pairs in qubit_pairs.batch(): for qp in multiplexed_qubit_pairs.values(): From f34a15f130ede20263979e99f1ddafaea0529bc9 Mon Sep 17 00:00:00 2001 From: paulQM Date: Thu, 18 Jun 2026 18:06:27 +0200 Subject: [PATCH 03/10] revert change --- .../calibrations/CZ_calibrations/32a_cz_conditional_phase.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/qualibration_graphs/superconducting/calibrations/CZ_calibrations/32a_cz_conditional_phase.py b/qualibration_graphs/superconducting/calibrations/CZ_calibrations/32a_cz_conditional_phase.py index 46e0b96da..f4f76a6bb 100644 --- a/qualibration_graphs/superconducting/calibrations/CZ_calibrations/32a_cz_conditional_phase.py +++ b/qualibration_graphs/superconducting/calibrations/CZ_calibrations/32a_cz_conditional_phase.py @@ -179,8 +179,8 @@ def create_qua_program(node: QualibrationNode[Parameters, Quam]): if node.parameters.use_state_discrimination: # measure g/e/f populations for both qubits - mq.readout_state_gef(state_mq[ii], pulse_name="readout") - sq.readout_state_gef(state_sq[ii], pulse_name="readout") + mq.readout_state_gef(state_mq[ii]) + sq.readout_state_gef(state_sq[ii]) with switch_(state_mq[ii]): with case_(0): wait(4) From d39a1e9d1c290fc7005b591ca99b236c5f8a5655 Mon Sep 17 00:00:00 2001 From: paulQM Date: Thu, 18 Jun 2026 18:07:51 +0200 Subject: [PATCH 04/10] revert --- .../calibrations/CZ_calibrations/32a_cz_conditional_phase.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/qualibration_graphs/superconducting/calibrations/CZ_calibrations/32a_cz_conditional_phase.py b/qualibration_graphs/superconducting/calibrations/CZ_calibrations/32a_cz_conditional_phase.py index f4f76a6bb..e7916ee15 100644 --- a/qualibration_graphs/superconducting/calibrations/CZ_calibrations/32a_cz_conditional_phase.py +++ b/qualibration_graphs/superconducting/calibrations/CZ_calibrations/32a_cz_conditional_phase.py @@ -73,10 +73,6 @@ @node.run_action(skip_if=node.modes.external) def custom_param(node: QualibrationNode[Parameters, Quam]): # You can get type hinting in your IDE by typing node.parameters. - node.parameters.amp_range = 0.06 - node.parameters.reset_type = "active" - node.parameters.operation = "cz_unipolar" - node.parameters.qubit_pairs = ["coupler_q4_q5"] pass From 9ceeb04546769f93d262e3b8418992841cb767fd Mon Sep 17 00:00:00 2001 From: paulQM Date: Tue, 30 Jun 2026 10:49:28 +0200 Subject: [PATCH 05/10] fix: black formatting --- .../superconducting/calibration_utils/cz_jazz2_n/plotting.py | 4 +++- .../superconducting/calibration_utils/cz_jazz_n/plotting.py | 4 +++- .../calibrations/CZ_calibrations/33b_JAZZ-N.py | 4 +--- 3 files changed, 7 insertions(+), 5 deletions(-) diff --git a/qualibration_graphs/superconducting/calibration_utils/cz_jazz2_n/plotting.py b/qualibration_graphs/superconducting/calibration_utils/cz_jazz2_n/plotting.py index 9e65a2cf8..1288829a3 100644 --- a/qualibration_graphs/superconducting/calibration_utils/cz_jazz2_n/plotting.py +++ b/qualibration_graphs/superconducting/calibration_utils/cz_jazz2_n/plotting.py @@ -60,7 +60,9 @@ def amp_abs_to_scale(a, abs_values=amps_abs, scale_values=amps_scale): # --- Bottom: averaged P_|00> with sinc fit --- if "p_avg" in fr.data_vars: - ax_avg.plot(amps_scale, fr["p_avg"].values, "o", ms=3, color="C0", label=r"$\langle P_{|00\rangle}\rangle_N$") + ax_avg.plot( + amps_scale, fr["p_avg"].values, "o", ms=3, color="C0", label=r"$\langle P_{|00\rangle}\rangle_N$" + ) if "sinc_fit" in fr.data_vars: fit_vals = fr["sinc_fit"].values if np.any(np.isfinite(fit_vals)): diff --git a/qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/plotting.py b/qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/plotting.py index dea67bb8a..be548011a 100644 --- a/qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/plotting.py +++ b/qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/plotting.py @@ -60,7 +60,9 @@ def amp_abs_to_scale(a, abs_values=amps_abs, scale_values=amps_scale): # --- Bottom: averaged P_|1> with sinc fit --- if "p_avg" in fr.data_vars: - ax_avg.plot(amps_scale, fr["p_avg"].values, "o", ms=3, color="C0", label=r"$\langle P_{|1\rangle}\rangle_N$") + ax_avg.plot( + amps_scale, fr["p_avg"].values, "o", ms=3, color="C0", label=r"$\langle P_{|1\rangle}\rangle_N$" + ) if "sinc_fit" in fr.data_vars: fit_vals = fr["sinc_fit"].values if np.any(np.isfinite(fit_vals)): diff --git a/qualibration_graphs/superconducting/calibrations/CZ_calibrations/33b_JAZZ-N.py b/qualibration_graphs/superconducting/calibrations/CZ_calibrations/33b_JAZZ-N.py index 686f7fad5..ce1e20f70 100644 --- a/qualibration_graphs/superconducting/calibrations/CZ_calibrations/33b_JAZZ-N.py +++ b/qualibration_graphs/superconducting/calibrations/CZ_calibrations/33b_JAZZ-N.py @@ -26,7 +26,6 @@ from qualibration_libs.runtime import simulate_and_plot from quam_config import Quam - # %% {Initialisation} description = """ JAZZ-N CZ AMPLITUDE CALIBRATION @@ -85,8 +84,7 @@ def create_qua_program(node: QualibrationNode[Parameters, Quam]): # pylint: dis if not node.parameters.use_state_discrimination: raise RuntimeError( - "JAZZ-N reads the target qubit |1> population and therefore requires " - "use_state_discrimination = True." + "JAZZ-N reads the target qubit |1> population and therefore requires " "use_state_discrimination = True." ) node.namespace["qubit_pairs"] = qubit_pairs = get_qubit_pairs(node) From 3ea6da3d7ff32c2893b59def792004f5e4785243 Mon Sep 17 00:00:00 2001 From: Deepak Khurana <119570568+Deepakkhurrana@users.noreply.github.com> Date: Wed, 1 Jul 2026 17:47:55 +0200 Subject: [PATCH 06/10] renaming and introducing moving and staionary qubits --- .../{33b_JAZZ-N.py => 32c_JAZZ-N.py} | 92 +++++++++---- .../{33c_JAZZ2-N.py => 32d_JAZZ2-N.py} | 124 ++++++++++++------ 2 files changed, 147 insertions(+), 69 deletions(-) rename qualibration_graphs/superconducting/calibrations/CZ_calibrations/{33b_JAZZ-N.py => 32c_JAZZ-N.py} (70%) rename qualibration_graphs/superconducting/calibrations/CZ_calibrations/{33c_JAZZ2-N.py => 32d_JAZZ2-N.py} (68%) diff --git a/qualibration_graphs/superconducting/calibrations/CZ_calibrations/33b_JAZZ-N.py b/qualibration_graphs/superconducting/calibrations/CZ_calibrations/32c_JAZZ-N.py similarity index 70% rename from qualibration_graphs/superconducting/calibrations/CZ_calibrations/33b_JAZZ-N.py rename to qualibration_graphs/superconducting/calibrations/CZ_calibrations/32c_JAZZ-N.py index ce1e20f70..6588aa29c 100644 --- a/qualibration_graphs/superconducting/calibrations/CZ_calibrations/33b_JAZZ-N.py +++ b/qualibration_graphs/superconducting/calibrations/CZ_calibrations/32c_JAZZ-N.py @@ -9,11 +9,13 @@ import xarray as xr from calibration_utils.cz_jazz_n import ( Parameters, + QubitRoles, coerce_to_4k_plus_1, fit_raw_data, log_fitted_results, plot_raw_data_with_fit, process_raw_dataset, + verify_moving_qubit, ) from qm.qua import * from qualang_tools.loops import from_array @@ -32,12 +34,12 @@ This node calibrates the CZ-pulse amplitude using the JAZZ-N protocol (arXiv:2402.18926v3, Appendix I.1, Fig. 13(a)). The pulse sequence is - x90(target) -- CZ -- [X_pi(control) & X_pi(target) -- CZ] x N -- x90(target) -- measure(target) + x90(stationary) -- CZ -- [X_pi(moving) & X_pi(stationary) -- CZ] x N -- x90(stationary) -- measure(stationary) where N = 4k + 1 (k = 0, 1, 2, ...). With the X_pi refocusing pulses on both -qubits, the target |1> population evolves as +qubits, the stationary |1> population evolves as - P_|1>(target) = (1 - cos((2k+1) * theta_CZ)) / 2, + P_|1>(stationary) = (1 - cos((2k+1) * theta_CZ)) / 2, independently of any virtual-Z (single-qubit) phase shifts inside the CZ macro. The optimal CZ amplitude is the value where theta_CZ = pi, i.e. where @@ -49,9 +51,19 @@ macro provides (cz_unipolar, cz_flattop, cz_bipolar, cz_flattop_erf, cz_SNZ) can be calibrated. +JAZZ-N is a precision fine-tuning upgrade over 32b. The X_pi refocusing +pulses echo out ordinary single-qubit phase (residual detuning, AC-Stark +shifts) accumulated over the sequence, so the extracted phase is purely +theta_CZ = theta_11 - theta_10 - theta_01 + theta_00, immune to control's +frequency calibration. This also means control's |0> and |1> conditions +don't need to be prepared and measured as two separate calibration runs -- +the echo sweeps control through both states within a single sequence and +cancels the state-independent part automatically, making the experiment +faster than the naive conditioned-tomography approach. + Prerequisites: - Calibrated single-qubit gates (x90, x180) for both qubits in the pair. - - Calibrated, state-discriminating readout for the target qubit. + - Calibrated, state-discriminating readout for the stationary qubit. - An initial estimate of the CZ amplitude (e.g. from 32a_cz_conditional_phase or 32b_cz_conditional_phase_error_amp). @@ -60,7 +72,7 @@ """ node = QualibrationNode[Parameters, Quam]( - name="33b_JAZZ_N", + name="32c_JAZZ_N", description=description, parameters=Parameters(), machine=Quam.load(), @@ -84,12 +96,19 @@ def create_qua_program(node: QualibrationNode[Parameters, Quam]): # pylint: dis if not node.parameters.use_state_discrimination: raise RuntimeError( - "JAZZ-N reads the target qubit |1> population and therefore requires " "use_state_discrimination = True." + "JAZZ-N reads the stationary qubit |1> population and therefore requires " + "use_state_discrimination = True." ) node.namespace["qubit_pairs"] = qubit_pairs = get_qubit_pairs(node) num_qubit_pairs = len(qubit_pairs) + qubit_roles_map = {} + for qp in qubit_pairs: + verify_moving_qubit(qp, operation=node.parameters.operation, log_callable=node.log) + qubit_roles_map[qp.name] = QubitRoles.resolve(qp) + node.namespace["qubit_roles_map"] = qubit_roles_map + # Coerce N_min / N_max to the nearest 4k + 1 (>= 1) and warn if changed. n_min_req = int(node.parameters.N_min) n_max_req = int(node.parameters.N_max) @@ -119,27 +138,31 @@ def create_qua_program(node: QualibrationNode[Parameters, Quam]): # pylint: dis n_op = declare(int) count = declare(int) n_st = declare_output_stream() - state_t = [declare(int) for _ in range(num_qubit_pairs)] - state_t_st = [declare_output_stream() for _ in range(num_qubit_pairs)] + state_sq = [declare(int) for _ in range(num_qubit_pairs)] + state_sq_st = [declare_output_stream() for _ in range(num_qubit_pairs)] for multiplexed_qubit_pairs in qubit_pairs.batch(): for qp in multiplexed_qubit_pairs.values(): - node.machine.initialize_qpu(target=qp.qubit_control) - node.machine.initialize_qpu(target=qp.qubit_target) + qubit_role = qubit_roles_map[qp.name] + mq, sq = qubit_role.moving, qubit_role.stationary + node.machine.initialize_qpu(target=mq) + node.machine.initialize_qpu(target=sq) with for_(n, 0, n < n_avg, n + 1): save(n, n_st) with for_(n_op, int(n_min), n_op <= int(n_max), n_op + 4): with for_(*from_array(amp, amplitudes)): for ii, qp in multiplexed_qubit_pairs.items(): - qp.qubit_control.reset(node.parameters.reset_type, node.parameters.simulate) - qp.qubit_target.reset(node.parameters.reset_type, node.parameters.simulate) + qubit_role = qubit_roles_map[qp.name] + mq, sq = qubit_role.moving, qubit_role.stationary + mq.reset(node.parameters.reset_type, node.parameters.simulate) + sq.reset(node.parameters.reset_type, node.parameters.simulate) qp.align() - reset_frame(qp.qubit_target.xy.name) - reset_frame(qp.qubit_control.xy.name) + reset_frame(sq.xy.name) + reset_frame(mq.xy.name) - # Initial pi/2 on target (Q2 role); control (Q1 role) stays in |0>. - qp.qubit_target.xy.play("x90") + # Initial pi/2 on stationary; moving stays in |0>. + sq.xy.play("x90") qp.align() # First CZ (the "Z" preceding the (pi-Z)^N pattern). @@ -148,24 +171,24 @@ def create_qua_program(node: QualibrationNode[Parameters, Quam]): # pylint: dis # N echoes interleaved with N more CZs: [X_pi X_pi, CZ] x N. with for_(count, 0, count < n_op, count + 1): qp.align() - qp.qubit_control.xy.play("x180") - qp.qubit_target.xy.play("x180") + mq.xy.play("x180") + sq.xy.play("x180") qp.align() qp.macros[operation].apply(amplitude_scale_qubit=amp) qp.align() - qp.qubit_target.xy.play("x90") + sq.xy.play("x90") qp.align() - qp.qubit_target.readout_state(state_t[ii]) - save(state_t[ii], state_t_st[ii]) + sq.readout_state(state_sq[ii]) + save(state_sq[ii], state_sq_st[ii]) align() with stream_processing(): n_st.save("n") for ii in range(num_qubit_pairs): - state_t_st[ii].buffer(len(amplitudes)).buffer(len(n_values)).average().save(f"state_target{ii + 1}") + state_sq_st[ii].buffer(len(amplitudes)).buffer(len(n_values)).average().save(f"state_stationary{ii + 1}") # %% {Simulate} @@ -195,6 +218,10 @@ def execute_qua_program(node: QualibrationNode[Parameters, Quam]): ) node.log(job.execution_report()) node.results["ds_raw"] = dataset + qubit_roles_map = node.namespace["qubit_roles_map"] + node.results["qubit_roles"] = { + name: {field: getattr(role, field).name for field in role._fields} for name, role in qubit_roles_map.items() + } # %% {Load_data} @@ -205,6 +232,15 @@ def load_data(node: QualibrationNode[Parameters, Quam]): node.load_from_id(node.parameters.load_data_id) node.parameters.load_data_id = load_data_id node.namespace["qubit_pairs"] = get_qubit_pairs(node) + if "qubit_roles" in node.results: + node.namespace["qubit_roles_map"] = { + name: QubitRoles(**{field: node.machine.qubits[qname] for field, qname in roles.items()}) + for name, roles in node.results["qubit_roles"].items() + } + else: + node.namespace["qubit_roles_map"] = { + qp.name: QubitRoles.resolve(qp) for qp in node.namespace["qubit_pairs"] + } # %% {Analyse_data} @@ -224,10 +260,14 @@ 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 JAZZ-N P_|1> map and the fitted optimal amplitude for each qubit pair.""" - fig = plot_raw_data_with_fit(node.results["ds_fit"], node.namespace["qubit_pairs"]) - plt.show() - node.results["figures"] = {"jazz_n_amplitude": fig} + """Plot the raw and fitted data in a specific figure whose shape is given by qubit pair grid locations.""" + figures = plot_raw_data_with_fit(node.results["ds_fit"], node.namespace["qubit_pairs"]) + for fig in figures.values(): + plt.show() + node.results["figures"] = { + "jazz_n_map": figures["map"], + "jazz_n_avg": figures["avg"], + } # %% {Update_state} diff --git a/qualibration_graphs/superconducting/calibrations/CZ_calibrations/33c_JAZZ2-N.py b/qualibration_graphs/superconducting/calibrations/CZ_calibrations/32d_JAZZ2-N.py similarity index 68% rename from qualibration_graphs/superconducting/calibrations/CZ_calibrations/33c_JAZZ2-N.py rename to qualibration_graphs/superconducting/calibrations/CZ_calibrations/32d_JAZZ2-N.py index ef7b056e6..c794768e5 100644 --- a/qualibration_graphs/superconducting/calibrations/CZ_calibrations/33c_JAZZ2-N.py +++ b/qualibration_graphs/superconducting/calibrations/CZ_calibrations/32d_JAZZ2-N.py @@ -9,11 +9,13 @@ import xarray as xr from calibration_utils.cz_jazz2_n import ( Parameters, + QubitRoles, coerce_to_even, fit_raw_data, log_fitted_results, plot_raw_data_with_fit, process_raw_dataset, + verify_moving_qubit, ) from qm.qua import * from qualang_tools.loops import from_array @@ -33,11 +35,11 @@ This node calibrates the CZ-pulse amplitude using the JAZZ2-N protocol (arXiv:2402.18926v3, Appendix I.1, Fig. 13(b)). The pulse sequence is - x90(control) & x90(target) (X_{pi/2} X_{pi/2}) + x90(moving) & x90(stationary) (X_{pi/2} X_{pi/2}) CZ (initial Z) - [X_pi(control) & X_pi(target) -- CZ] x (2N + 1) - x90(control) & x90(target) (X_{pi/2} X_{pi/2}) - measure(control), measure(target) -> p00 = (1 - state_c) * (1 - state_t) + [X_pi(moving) & X_pi(stationary) -- CZ] x (2N + 1) + x90(moving) & x90(stationary) (X_{pi/2} X_{pi/2}) + measure(moving), measure(stationary) -> p00 = (1 - state_moving) * (1 - state_stationary) where N = 2k (k = 0, 1, 2, ...). With the X_pi refocusing pulses on both qubits, the joint ground-state probability evolves as @@ -56,17 +58,24 @@ macro provides (cz_unipolar, cz_flattop, cz_bipolar, cz_flattop_erf, cz_SNZ) can be calibrated. +Compare to 32c, measuring both qubits in superposition together, rather than reading out +stationary qubit alone, makes the extracted phase more robust to single-qubit gate +errors: an imperfect x90/X_pi on either qubit is folded symmetrically into +the joint correlator instead of being dumped entirely onto one qubit's +readout, so small single-qubit miscalibration partially cancels rather +than biasing theta_CZ directly. + Prerequisites: - Calibrated single-qubit gates (x90, x180) for both qubits in the pair. - Calibrated, state-discriminating readout for BOTH qubits. - - An initial estimate of the CZ amplitude (e.g. from 33b_JAZZ-N). + - An initial estimate of the CZ amplitude (e.g. from 32a_cz_conditional_phase or 32b_cz_conditional_phase_error_amp). State update: - qubit_pair.macros[operation].flux_pulse_qubit.amplitude (fitted optimal CZ amplitude). """ node = QualibrationNode[Parameters, Quam]( - name="33c_JAZZ2_N", + name="32d_JAZZ2_N", description=description, parameters=Parameters(), machine=Quam.load(), @@ -97,6 +106,12 @@ def create_qua_program(node: QualibrationNode[Parameters, Quam]): # pylint: dis node.namespace["qubit_pairs"] = qubit_pairs = get_qubit_pairs(node) num_qubit_pairs = len(qubit_pairs) + qubit_roles_map = {} + for qp in qubit_pairs: + verify_moving_qubit(qp, operation=node.parameters.operation, log_callable=node.log) + qubit_roles_map[qp.name] = QubitRoles.resolve(qp) + node.namespace["qubit_roles_map"] = qubit_roles_map + # Coerce N_min / N_max to the nearest even integer (>= 0) and warn if changed. n_min_req = int(node.parameters.N_min) n_max_req = int(node.parameters.N_max) @@ -126,39 +141,43 @@ def create_qua_program(node: QualibrationNode[Parameters, Quam]): # pylint: dis n_op = declare(int) count = declare(int) n_st = declare_output_stream() - state_c = [declare(int) for _ in range(num_qubit_pairs)] - state_t = [declare(int) for _ in range(num_qubit_pairs)] + state_mq = [declare(int) for _ in range(num_qubit_pairs)] + state_sq = [declare(int) for _ in range(num_qubit_pairs)] p00 = [declare(int) for _ in range(num_qubit_pairs)] p00_st = [declare_output_stream() for _ in range(num_qubit_pairs)] - state_c_st = [declare_output_stream() for _ in range(num_qubit_pairs)] - state_t_st = [declare_output_stream() for _ in range(num_qubit_pairs)] + state_mq_st = [declare_output_stream() for _ in range(num_qubit_pairs)] + state_sq_st = [declare_output_stream() for _ in range(num_qubit_pairs)] for multiplexed_qubit_pairs in qubit_pairs.batch(): for qp in multiplexed_qubit_pairs.values(): - node.machine.initialize_qpu(target=qp.qubit_control) - node.machine.initialize_qpu(target=qp.qubit_target) + qubit_role = qubit_roles_map[qp.name] + mq, sq = qubit_role.moving, qubit_role.stationary + node.machine.initialize_qpu(target=mq) + node.machine.initialize_qpu(target=sq) with for_(n, 0, n < n_avg, n + 1): save(n, n_st) with for_(n_op, n_min, n_op <= n_max, n_op + 2): with for_(*from_array(amp, amplitudes)): for ii, qp in multiplexed_qubit_pairs.items(): - qp.qubit_control.reset(node.parameters.reset_type, node.parameters.simulate) - qp.qubit_target.reset(node.parameters.reset_type, node.parameters.simulate) + qubit_role = qubit_roles_map[qp.name] + mq, sq = qubit_role.moving, qubit_role.stationary + mq.reset(node.parameters.reset_type, node.parameters.simulate) + sq.reset(node.parameters.reset_type, node.parameters.simulate) qp.align() - reset_frame(qp.qubit_target.xy.name) - reset_frame(qp.qubit_control.xy.name) + reset_frame(sq.xy.name) + reset_frame(mq.xy.name) # Boundary X_{pi/2} X_{pi/2} (both qubits). - qp.qubit_control.xy.play("x90") - qp.qubit_target.xy.play("x90") + mq.xy.play("x90") + sq.xy.play("x90") qp.align() # First CZ (the "Z" preceding the (pi-Z)^(2N+1) pattern). qp.macros[operation].apply(amplitude_scale_qubit=amp) - qp.qubit_control.xy.play("x180") - qp.qubit_target.xy.play("x180") + mq.xy.play("x180") + sq.xy.play("x180") qp.align() # First CZ (the "Z" preceding the (pi-Z)^(2N+1) pattern). @@ -166,33 +185,33 @@ def create_qua_program(node: QualibrationNode[Parameters, Quam]): # pylint: dis # (X_pi X_pi, CZ) x (2N + 1). with for_(count, 1, count <= n_op, count + 1): - qp.qubit_control.xy.play("x180") - qp.qubit_target.xy.play("x180") + mq.xy.play("x180") + sq.xy.play("x180") qp.align() qp.macros[operation].apply(amplitude_scale_qubit=amp) - qp.qubit_control.xy.frame_rotation_2pi(0.5) - qp.qubit_target.xy.frame_rotation_2pi(0.5) - qp.qubit_control.xy.play("x180") - qp.qubit_target.xy.play("x180") + mq.xy.frame_rotation_2pi(0.5) + sq.xy.frame_rotation_2pi(0.5) + mq.xy.play("x180") + sq.xy.play("x180") qp.align() qp.macros[operation].apply(amplitude_scale_qubit=amp) - qp.qubit_control.xy.frame_rotation_2pi(-0.5) - qp.qubit_target.xy.frame_rotation_2pi(-0.5) + mq.xy.frame_rotation_2pi(-0.5) + sq.xy.frame_rotation_2pi(-0.5) qp.align() # Boundary X_{pi/2} X_{pi/2} (both qubits). - qp.qubit_control.xy.play("x90") - qp.qubit_control.xy.play("x180") - qp.qubit_target.xy.play("x90") - qp.qubit_target.xy.play("x180") + mq.xy.play("x90") + mq.xy.play("x180") + sq.xy.play("x90") + sq.xy.play("x180") qp.align() - qp.qubit_control.readout_state(state_c[ii]) - qp.qubit_target.readout_state(state_t[ii]) - assign(p00[ii], (1 - state_c[ii]) * (1 - state_t[ii])) + mq.readout_state(state_mq[ii]) + sq.readout_state(state_sq[ii]) + assign(p00[ii], (1 - state_mq[ii]) * (1 - state_sq[ii])) save(p00[ii], p00_st[ii]) - save(state_c[ii], state_c_st[ii]) - save(state_t[ii], state_t_st[ii]) + save(state_mq[ii], state_mq_st[ii]) + save(state_sq[ii], state_sq_st[ii]) align() @@ -200,8 +219,10 @@ def create_qua_program(node: QualibrationNode[Parameters, Quam]): # pylint: dis n_st.save("n") for ii in range(num_qubit_pairs): p00_st[ii].buffer(len(amplitudes)).buffer(len(n_values)).average().save(f"p{ii + 1}") - state_c_st[ii].buffer(len(amplitudes)).buffer(len(n_values)).average().save(f"state_c{ii + 1}") - state_t_st[ii].buffer(len(amplitudes)).buffer(len(n_values)).average().save(f"state_t{ii + 1}") + state_mq_st[ii].buffer(len(amplitudes)).buffer(len(n_values)).average().save(f"state_moving{ii + 1}") + state_sq_st[ii].buffer(len(amplitudes)).buffer(len(n_values)).average().save( + f"state_stationary{ii + 1}" + ) # %% {Simulate} @@ -231,6 +252,10 @@ def execute_qua_program(node: QualibrationNode[Parameters, Quam]): ) node.log(job.execution_report()) node.results["ds_raw"] = dataset + qubit_roles_map = node.namespace["qubit_roles_map"] + node.results["qubit_roles"] = { + name: {field: getattr(role, field).name for field in role._fields} for name, role in qubit_roles_map.items() + } # %% {Load_data} @@ -241,6 +266,15 @@ def load_data(node: QualibrationNode[Parameters, Quam]): node.load_from_id(node.parameters.load_data_id) node.parameters.load_data_id = load_data_id node.namespace["qubit_pairs"] = get_qubit_pairs(node) + if "qubit_roles" in node.results: + node.namespace["qubit_roles_map"] = { + name: QubitRoles(**{field: node.machine.qubits[qname] for field, qname in roles.items()}) + for name, roles in node.results["qubit_roles"].items() + } + else: + node.namespace["qubit_roles_map"] = { + qp.name: QubitRoles.resolve(qp) for qp in node.namespace["qubit_pairs"] + } # %% {Analyse_data} @@ -260,10 +294,14 @@ 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 JAZZ2-N P_|00> map and the fitted optimal amplitude for each qubit pair.""" - fig = plot_raw_data_with_fit(node.results["ds_fit"], node.namespace["qubit_pairs"]) - plt.show() - node.results["figures"] = {"jazz2_n_amplitude": fig} + """Plot the raw and fitted data in a specific figure whose shape is given by qubit pair grid locations.""" + figures = plot_raw_data_with_fit(node.results["ds_fit"], node.namespace["qubit_pairs"]) + for fig in figures.values(): + plt.show() + node.results["figures"] = { + "jazz2_n_map": figures["map"], + "jazz2_n_avg": figures["avg"], + } # %% {Update_state} From 7b16908d63590142e9eae314da8b3c7e959f2689 Mon Sep 17 00:00:00 2001 From: Deepak Khurana <119570568+Deepakkhurrana@users.noreply.github.com> Date: Wed, 1 Jul 2026 17:48:31 +0200 Subject: [PATCH 07/10] introducing moving and staionary qubits --- .../calibration_utils/cz_jazz2_n/__init__.py | 4 + .../calibration_utils/cz_jazz_n/__init__.py | 4 + .../calibration_utils/cz_jazz_n/analysis.py | 13 +- .../calibration_utils/cz_jazz_n/parameters.py | 4 +- .../two_qubit_interleaved_rb/analysis.py | 170 ++++++++++++++++++ 5 files changed, 186 insertions(+), 9 deletions(-) create mode 100644 qualibration_graphs/superconducting/calibration_utils/two_qubit_interleaved_rb/analysis.py diff --git a/qualibration_graphs/superconducting/calibration_utils/cz_jazz2_n/__init__.py b/qualibration_graphs/superconducting/calibration_utils/cz_jazz2_n/__init__.py index 2b557dd6f..8513a8138 100644 --- a/qualibration_graphs/superconducting/calibration_utils/cz_jazz2_n/__init__.py +++ b/qualibration_graphs/superconducting/calibration_utils/cz_jazz2_n/__init__.py @@ -1,5 +1,7 @@ """JAZZ2-N CZ amplitude calibration utilities.""" +from calibration_utils.cz_iswap_flux_bootstrap.parameters import QubitRoles, verify_moving_qubit # noqa: F401 + from .analysis import ( FitResults, coerce_to_even, @@ -13,9 +15,11 @@ __all__ = [ "FitResults", "Parameters", + "QubitRoles", "coerce_to_even", "fit_raw_data", "log_fitted_results", "plot_raw_data_with_fit", "process_raw_dataset", + "verify_moving_qubit", ] diff --git a/qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/__init__.py b/qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/__init__.py index 3ff54ea58..8593b5b70 100644 --- a/qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/__init__.py +++ b/qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/__init__.py @@ -1,5 +1,7 @@ """JAZZ-N CZ amplitude calibration utilities.""" +from calibration_utils.cz_iswap_flux_bootstrap.parameters import QubitRoles, verify_moving_qubit # noqa: F401 + from .analysis import ( FitResults, coerce_to_4k_plus_1, @@ -13,9 +15,11 @@ __all__ = [ "FitResults", "Parameters", + "QubitRoles", "coerce_to_4k_plus_1", "fit_raw_data", "log_fitted_results", "plot_raw_data_with_fit", "process_raw_dataset", + "verify_moving_qubit", ] diff --git a/qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/analysis.py b/qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/analysis.py index 5f9a3818e..25eb0d942 100644 --- a/qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/analysis.py +++ b/qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/analysis.py @@ -1,6 +1,6 @@ """Analysis module for the JAZZ-N CZ amplitude calibration. -The protocol measures the target qubit's |1> population as a function of the +The protocol measures the stationary qubit's |1> population as a function of the CZ-pulse amplitude scale, for several echo repetitions N = 4k + 1. Ignoring decoherence, equation (36) of arXiv:2402.18926v3 predicts: @@ -150,7 +150,7 @@ def _fit_one_pair( amp_values : (n_amp,) amplitude-scale values (centred at 1.0). n_values : (n_N,) echo counts N (sorted ascending; not used by the fit but retained to keep the call signature consistent with the previous version). - p_curve : (n_N, n_amp) target P_|1> values. + p_curve : (n_N, n_amp) stationary P_|1> values. Returns ------- @@ -231,7 +231,7 @@ def fit_raw_data(ds: xr.Dataset, node: QualibrationNode) -> Tuple[xr.Dataset, Di """Fit the JAZZ-N data per qubit pair. The dataset is augmented with three new data variables (per qubit pair, on - the amp axis): ``p_avg`` (mean over N of the target |1> population) and + the amp axis): ``p_avg`` (mean over N of the stationary |1> population) and ``sinc_fit`` (the fitted sinc model evaluated on the amp grid), plus ``optimal_amplitude`` / ``optimal_amplitude_scale`` / ``success`` / ``fit_method`` as coordinates. @@ -239,9 +239,8 @@ def fit_raw_data(ds: xr.Dataset, node: QualibrationNode) -> Tuple[xr.Dataset, Di qubit_pairs = node.namespace["qubit_pairs"] operation = node.parameters.operation - data_var = "state_target" if "state_target" in ds else None - if data_var is None: - raise RuntimeError("JAZZ-N analysis requires 'state_target' in the dataset (state discrimination).") + if "state_stationary" not in ds: + raise RuntimeError("JAZZ-N analysis requires 'state_stationary' in the dataset (state discrimination).") amp_values = ds.amp.values n_values = ds.N.values @@ -257,7 +256,7 @@ def fit_raw_data(ds: xr.Dataset, node: QualibrationNode) -> Tuple[xr.Dataset, Di for qp_name in qp_names: qp = next(qp for qp in qubit_pairs if qp.name == qp_name) - p = ds[data_var].sel(qubit_pair=qp_name).transpose("N", "amp").values + p = ds.state_stationary.sel(qubit_pair=qp_name).transpose("N", "amp").values amp_scale, success, method, params, p_avg, fit_curve = _fit_one_pair( amp_values, np.asarray(n_values), np.asarray(p) ) diff --git a/qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/parameters.py b/qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/parameters.py index 0f03e8be1..b833cb4ab 100644 --- a/qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/parameters.py +++ b/qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/parameters.py @@ -1,7 +1,7 @@ """Parameters module for the JAZZ-N CZ amplitude calibration. The JAZZ-N protocol (Appendix I.1, Fig. 13(a) of arXiv:2402.18926v3) measures -P_|1> of the target qubit after a refocused train of CZ gates. The number of +P_|1> of the stationary qubit after a refocused train of CZ gates. The number of X_pi echo pulses N must satisfy N = 4k + 1 (k = 0, 1, 2, ...), which gives a clean (2k+1)*theta_CZ phase accumulation that peaks at theta_CZ = pi. """ @@ -31,7 +31,7 @@ class NodeSpecificParameters(RunnableParameters): operation: Literal["cz_flattop", "cz_unipolar", "cz_bipolar", "cz_flattop_erf", "cz_SNZ"] = "cz_unipolar" """Name of the CZGate macro to drive in place of the bare Z pulse. Default is 'cz_unipolar'.""" use_state_discrimination: bool = True - """JAZZ-N reads P_|1> of the target qubit, which requires state discrimination. Setting this to False raises.""" + """JAZZ-N reads P_|1> of the stationary qubit, which requires state discrimination. Setting this to False raises.""" class Parameters( diff --git a/qualibration_graphs/superconducting/calibration_utils/two_qubit_interleaved_rb/analysis.py b/qualibration_graphs/superconducting/calibration_utils/two_qubit_interleaved_rb/analysis.py new file mode 100644 index 000000000..d81a2414d --- /dev/null +++ b/qualibration_graphs/superconducting/calibration_utils/two_qubit_interleaved_rb/analysis.py @@ -0,0 +1,170 @@ +"""Analysis utilities for two-qubit randomized benchmarking experiments. + +This module provides functions for processing and analyzing raw RB data, +including dataset processing and result logging. +""" + +import logging +from dataclasses import dataclass +from typing import Dict, Tuple + +import numpy as np +import xarray as xr +from qualibrate import QualibrationNode + +# @dataclass +# class FitResults: +# """Stores the relevant fit parameters for a single qubit pair in an RB experiment""" + +# optimal_amplitude: float +# success: bool + + +def log_fitted_results(fit_results: Dict[str, float], log_callable=None): + """ + Logs the node-specific fitted results for all qubit pairs. + + Parameters: + ----------- + fit_results : Dict[str, float] + Dictionary containing floats for each qubit pair. + log_callable : callable, optional + Logger for logging the fitted results. If None, a default logger is used. + """ + if log_callable is None: + log_callable = logging.getLogger(__name__).info + + for qp_name, fit_result in fit_results.items(): + s_qubit = f"Results for qubit pair {qp_name}: " + + s_alpha = f"\tFitted alpha: {fit_result['alpha']:.6f} a.u." + s_fidelity = f"\tFitted fidelity: {100*fit_result['fidelity']:.6f} %" + + if fit_result["success"]: + s_qubit += "SUCCESS!\n" + else: + s_qubit += "FAIL!\n" + + log_message = s_qubit + s_alpha + s_fidelity + + log_callable(log_message) + + +def process_raw_dataset(ds: xr.Dataset, node: QualibrationNode): + """ + Process the raw dataset by adding amplitude and detuning coordinates. + + Parameters: + ----------- + ds : xr.Dataset + Raw dataset from the experiment + node : QualibrationNode + The calibration node containing qubit pairs information + + Returns: + -------- + xr.Dataset + Processed dataset with additional coordinates + """ + ds = node.results["ds_raw"] + + rename_map = {"shots": "average", "sequence": "repeat", "depths": "circuit_depth"} + rename_map = {k: v for k, v in rename_map.items() if k in ds.dims} + + # Assume ds is your input dataset and ds['state'] is your DataArray + state = ds["state"] # shape: (qubit, shots, sequence, depths) + + # Outcome labels for 2-qubit states + labels = ["00", "01", "10", "11"] + + # Create a list of DataArrays: one for each outcome + probs = [state == i for i in range(4)] + + # Stack along a new outcome dimension + probs = xr.concat(probs, dim="outcome") + + # Assign outcome labels + probs = probs.assign_coords(outcome=("outcome", labels)) + + probs_00 = probs.sel(outcome="00") + if rename_map: + probs_00 = probs_00.rename(rename_map) + probs_00 = probs_00.transpose("qubit_pair", "repeat", "circuit_depth", "average") + + probs_00 = probs_00.astype(int) + + if rename_map: + ds_transposed = ds.rename(rename_map) + else: + ds_transposed = ds + ds_transposed = ds_transposed.transpose("qubit_pair", "repeat", "circuit_depth", "average") + + return ds_transposed + + +# def fit_raw_data(ds: xr.Dataset, node: QualibrationNode) -> Tuple[xr.Dataset, Dict[str, FitResults]]: +# """ +# Fit the CZ conditional phase data for each qubit pair. + +# Parameters: +# ----------- +# ds : xr.Dataset +# Dataset containing the processed data. +# node : QualibrationNode +# The calibration node containing parameters and qubit pairs. + +# Returns: +# -------- +# Tuple[xr.Dataset, Dict[str, FitResults]] +# Dataset with fit results and dictionary of fit results for each qubit pair. +# """ +# # For RB analysis, no fitting routine is currently implemented. +# # Pass the dataset through unchanged, or implement RB-specific fitting here if needed. +# ds_fit = ds + +# # Extract the relevant fitted parameters +# ds_fit, fit_results = _extract_relevant_parameters(ds_fit, node) + +# return ds_fit, fit_results + + +# def _extract_relevant_parameters( +# ds_fit: xr.Dataset, node: QualibrationNode +# ) -> Tuple[xr.Dataset, Dict[str, FitResults]]: +# """ +# Extract relevant fit parameters and create FitResults for each qubit pair. + +# Parameters: +# ----------- +# ds_fit : xr.Dataset +# Dataset containing the fit results from fit_routine. +# node : QualibrationNode +# The calibration node containing parameters and qubit pairs. + +# Returns: +# -------- +# Tuple[xr.Dataset, Dict[str, FitResults]] +# Dataset with additional metadata and dictionary of FitResults for each qubit pair. +# """ +# qubit_pairs = node.namespace["qubit_pairs"] + +# # Add metadata attributes to the dataset +# if "optimal_amplitude" in ds_fit.data_vars: +# ds_fit.optimal_amplitude.attrs = {"long_name": "optimal CZ amplitude", "units": "a.u."} +# if "phase_diff" in ds_fit.data_vars: +# ds_fit.phase_diff.attrs = {"long_name": "phase difference", "units": "2π"} +# if "fitted_curve" in ds_fit.data_vars: +# ds_fit.fitted_curve.attrs = {"long_name": "fitted tanh curve", "units": "2π"} + +# # Create FitResults for each qubit pair +# fit_results = {} +# for qp in qubit_pairs: +# qp_name = qp.name +# qp_data = ds_fit.sel(qubit_pair=qp_name) + +# fit_results[qp_name] = FitResults( +# optimal_amplitude=float(qp_data.optimal_amplitude.values), +# success=bool(qp_data.success.values), +# ) + +# return ds_fit, fit_results From c80cc9f933a7d4db737be932879517ede93f659a Mon Sep 17 00:00:00 2001 From: Deepak Khurana <119570568+Deepakkhurrana@users.noreply.github.com> Date: Wed, 1 Jul 2026 17:49:02 +0200 Subject: [PATCH 08/10] Making plotting consustent with pair grid --- .../calibration_utils/cz_jazz2_n/plotting.py | 199 +++++++++++------- .../calibration_utils/cz_jazz_n/plotting.py | 198 ++++++++++------- 2 files changed, 238 insertions(+), 159 deletions(-) diff --git a/qualibration_graphs/superconducting/calibration_utils/cz_jazz2_n/plotting.py b/qualibration_graphs/superconducting/calibration_utils/cz_jazz2_n/plotting.py index 1288829a3..cfb63302f 100644 --- a/qualibration_graphs/superconducting/calibration_utils/cz_jazz2_n/plotting.py +++ b/qualibration_graphs/superconducting/calibration_utils/cz_jazz2_n/plotting.py @@ -1,88 +1,129 @@ """Plotting module for the JAZZ2-N CZ amplitude calibration.""" -import matplotlib.pyplot as plt +from typing import Dict + import numpy as np import xarray as xr -from qualibration_libs.core import BatchableList +from matplotlib.axes import Axes +from matplotlib.figure import Figure +from qualibration_libs.plotting import grid_iter + +from calibration_utils.pair_grid import QubitPairGrid, grid_pair_names + +def plot_raw_data_with_fit( + ds_fit: xr.Dataset, + qubit_pairs: list, + title_prefix: str = "JAZZ2-N CZ amplitude calibration", +) -> Dict[str, Figure]: + """Plot the JAZZ2-N data and sinc fit per qubit pair on chip-topology grids. -def plot_raw_data_with_fit(ds_fit: xr.Dataset, qubit_pairs: BatchableList) -> plt.Figure: - """Plot the JAZZ2-N data and sinc fit per qubit pair. + For each qubit pair we show two separate figures laid out by chip topology: - For each qubit pair we show two stacked panels: + * ``"map"``: joint ``p`` (P_|00>) as a 2D map versus repetition count N = 2k + and amplitude scale, with the fitted optimal amplitude drawn as a vertical line. + * ``"avg"``: the averaged-over-N curve ``p_avg(amp)`` together with the + fitted sinc model and the optimum. - * Top: ``p00`` as a 2D map versus the repetition count N = 2k and the - amplitude scale, with the fitted optimal amplitude drawn as a vertical - line. - * Bottom: the averaged-over-N curve ``p_avg(amp)`` together with the - fitted sinc model and the optimum, so the fit can be eyeballed. + Parameters + ---------- + ds_fit : xr.Dataset + Fit dataset containing ``p``, ``p_avg``, ``sinc_fit``, + ``optimal_amplitude_scale``, and ``fit_method``. + qubit_pairs : list + Qubit pair objects used for grid placement. + title_prefix : str + Prefix for figure suptitles. + + Returns + ------- + dict[str, Figure] + ``"map"`` contains the 2D heatmaps and ``"avg"`` contains the averaged + curves with sinc fits; both include optimal-amplitude markers. """ - n_pairs = len(qubit_pairs) - cols = min(4, n_pairs) - rows = (n_pairs + cols - 1) // cols - fig, axes = plt.subplots(2 * rows, cols, figsize=(4.5 * cols, 5.5 * rows), squeeze=False) - - for i, qp in enumerate(qubit_pairs): - row, col = divmod(i, cols) - ax_map = axes[2 * row, col] - ax_avg = axes[2 * row + 1, col] - qp_name = qp.name - fr = ds_fit.sel(qubit_pair=qp_name) - - amps_scale = fr.amp.values - amps_abs = fr["amp_full"].values if "amp_full" in fr.coords else amps_scale - n_values = fr.N.values - - # --- Top: 2D heatmap of P_|00> --- - p_map = fr["p"].transpose("N", "amp") - xg, yg = np.meshgrid(amps_scale, n_values) - pcm = ax_map.pcolormesh(xg, yg, p_map.values, cmap="magma", shading="auto") - - opt_scale = float(fr.optimal_amplitude_scale.values) - opt_method = str(fr.fit_method.values) - if np.isfinite(opt_scale): - ax_map.axvline(opt_scale, color="lime", lw=2, label=f"opt = {opt_scale:.4f} ({opt_method})") - - def amp_scale_to_abs(s, abs_values=amps_abs, scale_values=amps_scale): - return np.interp(s, scale_values, abs_values) - - def amp_abs_to_scale(a, abs_values=amps_abs, scale_values=amps_scale): - return np.interp(a, abs_values, scale_values) - - secax = ax_map.secondary_xaxis("top", functions=(amp_scale_to_abs, amp_abs_to_scale)) - secax.set_xlabel("Amplitude (V)") - ax_map.set_title(qp_name) - ax_map.set_xlabel("Amplitude scale (a.u.)") - ax_map.set_ylabel("Repetition N = 2k") - ax_map.legend(loc="upper right", fontsize=8) - cbar = fig.colorbar(pcm, ax=ax_map, shrink=0.85) - cbar.set_label("$P_{|00\\rangle}$") - - # --- Bottom: averaged P_|00> with sinc fit --- - if "p_avg" in fr.data_vars: - ax_avg.plot( - amps_scale, fr["p_avg"].values, "o", ms=3, color="C0", label=r"$\langle P_{|00\rangle}\rangle_N$" - ) - if "sinc_fit" in fr.data_vars: - fit_vals = fr["sinc_fit"].values - if np.any(np.isfinite(fit_vals)): - ax_avg.plot(amps_scale, fit_vals, "-", lw=1.5, color="C3", label="sinc fit") - if np.isfinite(opt_scale): - ax_avg.axvline(opt_scale, color="lime", lw=1.5, label=f"opt = {opt_scale:.4f}") - ax_avg.set_xlabel("Amplitude scale (a.u.)") - ax_avg.set_ylabel(r"$\langle P_{|00\rangle}\rangle_N$") - ax_avg.legend(loc="upper right", fontsize=8) - - used = set() - for i in range(n_pairs): - row, col = divmod(i, cols) - used.add((2 * row, col)) - used.add((2 * row + 1, col)) - for r in range(axes.shape[0]): - for c in range(axes.shape[1]): - if (r, c) not in used: - axes[r, c].axis("off") - - fig.suptitle("JAZZ2-N CZ amplitude calibration") - fig.tight_layout(rect=(0, 0, 1, 0.97)) - return fig + grid_names, pair_names = grid_pair_names(qubit_pairs) + figures = {} + + map_grid = QubitPairGrid(grid_names, pair_names) + for ax, qubit in grid_iter(map_grid): + qp_name = qubit["qubit"] + plot_individual_map_with_fit(ax, ds_fit, qp_name) + map_grid.fig.suptitle(fr"{title_prefix} — $P_{{|00\rangle}}$ vs $N$ and amplitude") + map_grid.fig.tight_layout() + figures["map"] = map_grid.fig + + avg_grid = QubitPairGrid(grid_names, pair_names) + for ax, qubit in grid_iter(avg_grid): + qp_name = qubit["qubit"] + plot_individual_avg_with_fit(ax, ds_fit, qp_name) + avg_grid.fig.suptitle(fr"{title_prefix} — averaged $P_{{|00\rangle}}$ and sinc fit") + avg_grid.fig.tight_layout() + figures["avg"] = avg_grid.fig + + return figures + + +def plot_individual_map_with_fit(ax: Axes, ds_fit: xr.Dataset, qp_name: str) -> None: + """Plot one qubit-pair JAZZ2-N heatmap of P_|00>.""" + fr = ds_fit.sel(qubit_pair=qp_name) + + if "p" not in fr: + ax.text(0.5, 0.5, "No P(00) data", ha="center", va="center", transform=ax.transAxes) + ax.set_title(qp_name) + return + + amps_scale = fr.amp.values + amps_abs = fr["amp_full"].values if "amp_full" in fr.coords else amps_scale + n_values = fr.N.values + + p_map = fr["p"].transpose("N", "amp") + xg, yg = np.meshgrid(amps_scale, n_values) + pcm = ax.pcolormesh(xg, yg, p_map.values, cmap="magma", shading="auto") + + opt_scale = float(fr.optimal_amplitude_scale.values) + opt_method = str(fr.fit_method.values) + success = "success" in fr.coords and bool(fr.success) and np.isfinite(opt_scale) + if success: + ax.axvline(opt_scale, color="lime", lw=2, label=f"opt = {opt_scale:.4f} ({opt_method})") + + def amp_scale_to_abs(s, abs_values=amps_abs, scale_values=amps_scale): + return np.interp(s, scale_values, abs_values) + + def amp_abs_to_scale(a, abs_values=amps_abs, scale_values=amps_scale): + return np.interp(a, abs_values, scale_values) + + secax = ax.secondary_xaxis("top", functions=(amp_scale_to_abs, amp_abs_to_scale)) + secax.set_xlabel("Amplitude (V)") + ax.set_title(qp_name if success else f"{qp_name} — fit failed") + ax.set_xlabel("Amplitude scale (a.u.)") + ax.set_ylabel("Repetition N = 2k") + if success: + ax.legend(loc="upper right", fontsize=8) + ax.figure.colorbar(pcm, ax=ax, shrink=0.85).set_label("$P_{|00\\rangle}$") + + +def plot_individual_avg_with_fit(ax: Axes, ds_fit: xr.Dataset, qp_name: str) -> None: + """Plot one qubit-pair averaged P_|00> curve with sinc fit.""" + fr = ds_fit.sel(qubit_pair=qp_name) + + if "p_avg" not in fr: + ax.text(0.5, 0.5, "No averaged data", ha="center", va="center", transform=ax.transAxes) + ax.set_title(qp_name) + return + + amps_scale = fr.amp.values + opt_scale = float(fr.optimal_amplitude_scale.values) + success = "success" in fr.coords and bool(fr.success) and np.isfinite(opt_scale) + + ax.plot(amps_scale, fr["p_avg"].values, "o", ms=3, color="C0", label=r"$\langle P_{|00\rangle}\rangle_N$") + if "sinc_fit" in fr: + fit_vals = fr["sinc_fit"].values + if np.any(np.isfinite(fit_vals)): + ax.plot(amps_scale, fit_vals, "-", lw=1.5, color="C3", label="sinc fit") + if success: + ax.axvline(opt_scale, color="lime", lw=1.5, label=f"opt = {opt_scale:.4f}") + + ax.set_title(qp_name if success else f"{qp_name} — fit failed") + ax.set_xlabel("Amplitude scale (a.u.)") + ax.set_ylabel(r"$\langle P_{|00\rangle}\rangle_N$") + ax.legend(loc="upper right", fontsize=8) diff --git a/qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/plotting.py b/qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/plotting.py index be548011a..0d26bd7aa 100644 --- a/qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/plotting.py +++ b/qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/plotting.py @@ -1,91 +1,129 @@ """Plotting module for the JAZZ-N CZ amplitude calibration.""" -import matplotlib.pyplot as plt +from typing import Dict + import numpy as np import xarray as xr -from qualibration_libs.core import BatchableList +from matplotlib.axes import Axes +from matplotlib.figure import Figure +from qualibration_libs.plotting import grid_iter + +from calibration_utils.pair_grid import QubitPairGrid, grid_pair_names -def plot_raw_data_with_fit(ds_fit: xr.Dataset, qubit_pairs: BatchableList) -> plt.Figure: - """Plot the JAZZ-N data and sinc fit per qubit pair. +def plot_raw_data_with_fit( + ds_fit: xr.Dataset, + qubit_pairs: list, + title_prefix: str = "JAZZ-N CZ amplitude calibration", +) -> Dict[str, Figure]: + """Plot the JAZZ-N data and sinc fit per qubit pair on chip-topology grids. - For each qubit pair we show two stacked panels: + For each qubit pair we show two separate figures laid out by chip topology: - * Top: ``state_target`` (target |1> population) as a 2D map versus the + * ``"map"``: ``state_stationary`` (stationary |1> population) as a 2D map versus the echo count N and the amplitude scale, with the fitted optimal amplitude drawn as a vertical line. - * Bottom: the averaged-over-N curve ``p_avg(amp)`` together with the - fitted sinc model and the optimum, so the fit can be eyeballed. + * ``"avg"``: the averaged-over-N curve ``p_avg(amp)`` together with the + fitted sinc model and the optimum. + + Parameters + ---------- + ds_fit : xr.Dataset + Fit dataset containing ``state_stationary``, ``p_avg``, ``sinc_fit``, + ``optimal_amplitude_scale``, and ``fit_method``. + qubit_pairs : list + Qubit pair objects used for grid placement. + title_prefix : str + Prefix for figure suptitles. + + Returns + ------- + dict[str, Figure] + ``"map"`` contains the 2D heatmaps and ``"avg"`` contains the averaged + curves with sinc fits; both include optimal-amplitude markers. """ - n_pairs = len(qubit_pairs) - cols = min(4, n_pairs) - rows = (n_pairs + cols - 1) // cols - fig, axes = plt.subplots(2 * rows, cols, figsize=(4.5 * cols, 5.5 * rows), squeeze=False) - - for i, qp in enumerate(qubit_pairs): - row, col = divmod(i, cols) - ax_map = axes[2 * row, col] - ax_avg = axes[2 * row + 1, col] - qp_name = qp.name - fr = ds_fit.sel(qubit_pair=qp_name) - - amps_scale = fr.amp.values - amps_abs = fr["amp_full"].values if "amp_full" in fr.coords else amps_scale - n_values = fr.N.values - - # --- Top: 2D heatmap of P_|1> --- - p_map = fr["state_target"].transpose("N", "amp") - xg, yg = np.meshgrid(amps_scale, n_values) - pcm = ax_map.pcolormesh(xg, yg, p_map.values, cmap="magma", shading="auto") - - opt_scale = float(fr.optimal_amplitude_scale.values) - opt_method = str(fr.fit_method.values) - if np.isfinite(opt_scale): - ax_map.axvline(opt_scale, color="lime", lw=2, label=f"opt = {opt_scale:.4f} ({opt_method})") - - def amp_scale_to_abs(s, abs_values=amps_abs, scale_values=amps_scale): - return np.interp(s, scale_values, abs_values) - - def amp_abs_to_scale(a, abs_values=amps_abs, scale_values=amps_scale): - return np.interp(a, abs_values, scale_values) - - secax = ax_map.secondary_xaxis("top", functions=(amp_scale_to_abs, amp_abs_to_scale)) - secax.set_xlabel("Amplitude (V)") - ax_map.set_title(qp_name) - ax_map.set_xlabel("Amplitude scale (a.u.)") - ax_map.set_ylabel("Echo count N = 4k + 1") - ax_map.legend(loc="upper right", fontsize=8) - cbar = fig.colorbar(pcm, ax=ax_map, shrink=0.85) - cbar.set_label("$P_{|1\\rangle}$ of target") - - # --- Bottom: averaged P_|1> with sinc fit --- - if "p_avg" in fr.data_vars: - ax_avg.plot( - amps_scale, fr["p_avg"].values, "o", ms=3, color="C0", label=r"$\langle P_{|1\rangle}\rangle_N$" - ) - if "sinc_fit" in fr.data_vars: - fit_vals = fr["sinc_fit"].values - if np.any(np.isfinite(fit_vals)): - ax_avg.plot(amps_scale, fit_vals, "-", lw=1.5, color="C3", label="sinc fit") - if np.isfinite(opt_scale): - ax_avg.axvline(opt_scale, color="lime", lw=1.5, label=f"opt = {opt_scale:.4f}") - ax_avg.set_xlabel("Amplitude scale (a.u.)") - ax_avg.set_ylabel(r"$\langle P_{|1\rangle}\rangle_N$") - ax_avg.legend(loc="upper right", fontsize=8) - - # Hide unused panels. - total_axes = axes.flatten() - used = set() - for i in range(n_pairs): - row, col = divmod(i, cols) - used.add((2 * row, col)) - used.add((2 * row + 1, col)) - for r in range(axes.shape[0]): - for c in range(axes.shape[1]): - if (r, c) not in used: - axes[r, c].axis("off") - del total_axes - - fig.suptitle("JAZZ-N CZ amplitude calibration") - fig.tight_layout(rect=(0, 0, 1, 0.97)) - return fig + grid_names, pair_names = grid_pair_names(qubit_pairs) + figures = {} + + map_grid = QubitPairGrid(grid_names, pair_names) + for ax, qubit in grid_iter(map_grid): + qp_name = qubit["qubit"] + plot_individual_map_with_fit(ax, ds_fit, qp_name) + map_grid.fig.suptitle(fr"{title_prefix} — $P_{{|1\rangle}}$ vs $N$ and amplitude") + map_grid.fig.tight_layout() + figures["map"] = map_grid.fig + + avg_grid = QubitPairGrid(grid_names, pair_names) + for ax, qubit in grid_iter(avg_grid): + qp_name = qubit["qubit"] + plot_individual_avg_with_fit(ax, ds_fit, qp_name) + avg_grid.fig.suptitle(fr"{title_prefix} — averaged $\langle P_{{|1\rangle}}\rangle_N$ and sinc fit") + avg_grid.fig.tight_layout() + figures["avg"] = avg_grid.fig + + return figures + + +def plot_individual_map_with_fit(ax: Axes, ds_fit: xr.Dataset, qp_name: str) -> None: + """Plot one qubit-pair JAZZ-N heatmap of stationary |1> population.""" + fr = ds_fit.sel(qubit_pair=qp_name) + if "state_stationary" not in fr: + ax.text(0.5, 0.5, "No state data", ha="center", va="center", transform=ax.transAxes) + ax.set_title(qp_name) + return + + amps_scale = fr.amp.values + amps_abs = fr["amp_full"].values if "amp_full" in fr.coords else amps_scale + n_values = fr.N.values + + p_map = fr["state_stationary"].transpose("N", "amp") + xg, yg = np.meshgrid(amps_scale, n_values) + pcm = ax.pcolormesh(xg, yg, p_map.values, cmap="magma", shading="auto") + + opt_scale = float(fr.optimal_amplitude_scale.values) + opt_method = str(fr.fit_method.values) + success = "success" in fr.coords and bool(fr.success) and np.isfinite(opt_scale) + if success: + ax.axvline(opt_scale, color="lime", lw=2, label=f"opt = {opt_scale:.4f} ({opt_method})") + + def amp_scale_to_abs(s, abs_values=amps_abs, scale_values=amps_scale): + return np.interp(s, scale_values, abs_values) + + def amp_abs_to_scale(a, abs_values=amps_abs, scale_values=amps_scale): + return np.interp(a, abs_values, scale_values) + + secax = ax.secondary_xaxis("top", functions=(amp_scale_to_abs, amp_abs_to_scale)) + secax.set_xlabel("Amplitude (V)") + ax.set_title(qp_name if success else f"{qp_name} — fit failed") + ax.set_xlabel("Amplitude scale (a.u.)") + ax.set_ylabel("Echo count N = 4k + 1") + if success: + ax.legend(loc="upper right", fontsize=8) + ax.figure.colorbar(pcm, ax=ax, shrink=0.85).set_label("$P_{|1\\rangle}$ of stationary") + + +def plot_individual_avg_with_fit(ax: Axes, ds_fit: xr.Dataset, qp_name: str) -> None: + """Plot one qubit-pair averaged P_|1> curve with sinc fit.""" + fr = ds_fit.sel(qubit_pair=qp_name) + + if "p_avg" not in fr: + ax.text(0.5, 0.5, "No averaged data", ha="center", va="center", transform=ax.transAxes) + ax.set_title(qp_name) + return + + amps_scale = fr.amp.values + opt_scale = float(fr.optimal_amplitude_scale.values) + success = "success" in fr.coords and bool(fr.success) and np.isfinite(opt_scale) + + ax.plot(amps_scale, fr["p_avg"].values, "o", ms=3, color="C0", label=r"$\langle P_{|1\rangle}\rangle_N$") + if "sinc_fit" in fr: + fit_vals = fr["sinc_fit"].values + if np.any(np.isfinite(fit_vals)): + ax.plot(amps_scale, fit_vals, "-", lw=1.5, color="C3", label="sinc fit") + if success: + ax.axvline(opt_scale, color="lime", lw=1.5, label=f"opt = {opt_scale:.4f}") + + ax.set_title(qp_name if success else f"{qp_name} — fit failed") + ax.set_xlabel("Amplitude scale (a.u.)") + ax.set_ylabel(r"$\langle P_{|1\rangle}\rangle_N$") + ax.legend(loc="upper right", fontsize=8) From 3505d6acefe21a00604b8bfcbecb8da8ede11b71 Mon Sep 17 00:00:00 2001 From: Deepak Khurana <119570568+Deepakkhurrana@users.noreply.github.com> Date: Wed, 1 Jul 2026 17:51:09 +0200 Subject: [PATCH 09/10] Black formatting --- .../calibration_utils/cz_jazz2_n/plotting.py | 4 ++-- .../calibration_utils/cz_jazz_n/plotting.py | 4 ++-- .../calibrations/CZ_calibrations/32c_JAZZ-N.py | 8 ++++---- .../calibrations/CZ_calibrations/32d_JAZZ2-N.py | 5 +---- 4 files changed, 9 insertions(+), 12 deletions(-) diff --git a/qualibration_graphs/superconducting/calibration_utils/cz_jazz2_n/plotting.py b/qualibration_graphs/superconducting/calibration_utils/cz_jazz2_n/plotting.py index cfb63302f..6d11789ca 100644 --- a/qualibration_graphs/superconducting/calibration_utils/cz_jazz2_n/plotting.py +++ b/qualibration_graphs/superconducting/calibration_utils/cz_jazz2_n/plotting.py @@ -48,7 +48,7 @@ def plot_raw_data_with_fit( for ax, qubit in grid_iter(map_grid): qp_name = qubit["qubit"] plot_individual_map_with_fit(ax, ds_fit, qp_name) - map_grid.fig.suptitle(fr"{title_prefix} — $P_{{|00\rangle}}$ vs $N$ and amplitude") + map_grid.fig.suptitle(rf"{title_prefix} — $P_{{|00\rangle}}$ vs $N$ and amplitude") map_grid.fig.tight_layout() figures["map"] = map_grid.fig @@ -56,7 +56,7 @@ def plot_raw_data_with_fit( for ax, qubit in grid_iter(avg_grid): qp_name = qubit["qubit"] plot_individual_avg_with_fit(ax, ds_fit, qp_name) - avg_grid.fig.suptitle(fr"{title_prefix} — averaged $P_{{|00\rangle}}$ and sinc fit") + avg_grid.fig.suptitle(rf"{title_prefix} — averaged $P_{{|00\rangle}}$ and sinc fit") avg_grid.fig.tight_layout() figures["avg"] = avg_grid.fig diff --git a/qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/plotting.py b/qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/plotting.py index 0d26bd7aa..3fa4857cf 100644 --- a/qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/plotting.py +++ b/qualibration_graphs/superconducting/calibration_utils/cz_jazz_n/plotting.py @@ -49,7 +49,7 @@ def plot_raw_data_with_fit( for ax, qubit in grid_iter(map_grid): qp_name = qubit["qubit"] plot_individual_map_with_fit(ax, ds_fit, qp_name) - map_grid.fig.suptitle(fr"{title_prefix} — $P_{{|1\rangle}}$ vs $N$ and amplitude") + map_grid.fig.suptitle(rf"{title_prefix} — $P_{{|1\rangle}}$ vs $N$ and amplitude") map_grid.fig.tight_layout() figures["map"] = map_grid.fig @@ -57,7 +57,7 @@ def plot_raw_data_with_fit( for ax, qubit in grid_iter(avg_grid): qp_name = qubit["qubit"] plot_individual_avg_with_fit(ax, ds_fit, qp_name) - avg_grid.fig.suptitle(fr"{title_prefix} — averaged $\langle P_{{|1\rangle}}\rangle_N$ and sinc fit") + avg_grid.fig.suptitle(rf"{title_prefix} — averaged $\langle P_{{|1\rangle}}\rangle_N$ and sinc fit") avg_grid.fig.tight_layout() figures["avg"] = avg_grid.fig diff --git a/qualibration_graphs/superconducting/calibrations/CZ_calibrations/32c_JAZZ-N.py b/qualibration_graphs/superconducting/calibrations/CZ_calibrations/32c_JAZZ-N.py index 6588aa29c..c1d7ac051 100644 --- a/qualibration_graphs/superconducting/calibrations/CZ_calibrations/32c_JAZZ-N.py +++ b/qualibration_graphs/superconducting/calibrations/CZ_calibrations/32c_JAZZ-N.py @@ -188,7 +188,9 @@ def create_qua_program(node: QualibrationNode[Parameters, Quam]): # pylint: dis with stream_processing(): n_st.save("n") for ii in range(num_qubit_pairs): - state_sq_st[ii].buffer(len(amplitudes)).buffer(len(n_values)).average().save(f"state_stationary{ii + 1}") + state_sq_st[ii].buffer(len(amplitudes)).buffer(len(n_values)).average().save( + f"state_stationary{ii + 1}" + ) # %% {Simulate} @@ -238,9 +240,7 @@ def load_data(node: QualibrationNode[Parameters, Quam]): for name, roles in node.results["qubit_roles"].items() } else: - node.namespace["qubit_roles_map"] = { - qp.name: QubitRoles.resolve(qp) for qp in node.namespace["qubit_pairs"] - } + node.namespace["qubit_roles_map"] = {qp.name: QubitRoles.resolve(qp) for qp in node.namespace["qubit_pairs"]} # %% {Analyse_data} diff --git a/qualibration_graphs/superconducting/calibrations/CZ_calibrations/32d_JAZZ2-N.py b/qualibration_graphs/superconducting/calibrations/CZ_calibrations/32d_JAZZ2-N.py index c794768e5..99f1eba64 100644 --- a/qualibration_graphs/superconducting/calibrations/CZ_calibrations/32d_JAZZ2-N.py +++ b/qualibration_graphs/superconducting/calibrations/CZ_calibrations/32d_JAZZ2-N.py @@ -28,7 +28,6 @@ from qualibration_libs.runtime import simulate_and_plot from quam_config import Quam - # %% {Initialisation} description = """ JAZZ2-N CZ AMPLITUDE CALIBRATION @@ -272,9 +271,7 @@ def load_data(node: QualibrationNode[Parameters, Quam]): for name, roles in node.results["qubit_roles"].items() } else: - node.namespace["qubit_roles_map"] = { - qp.name: QubitRoles.resolve(qp) for qp in node.namespace["qubit_pairs"] - } + node.namespace["qubit_roles_map"] = {qp.name: QubitRoles.resolve(qp) for qp in node.namespace["qubit_pairs"]} # %% {Analyse_data} From 53f4d0545855f45bed48ab7e71d17b3dca72bf13 Mon Sep 17 00:00:00 2001 From: Deepak Khurana <119570568+Deepakkhurrana@users.noreply.github.com> Date: Thu, 2 Jul 2026 09:42:03 +0200 Subject: [PATCH 10/10] reverting unwanted changes --- .../two_qubit_interleaved_rb/analysis.py | 170 ------------------ 1 file changed, 170 deletions(-) delete mode 100644 qualibration_graphs/superconducting/calibration_utils/two_qubit_interleaved_rb/analysis.py diff --git a/qualibration_graphs/superconducting/calibration_utils/two_qubit_interleaved_rb/analysis.py b/qualibration_graphs/superconducting/calibration_utils/two_qubit_interleaved_rb/analysis.py deleted file mode 100644 index d81a2414d..000000000 --- a/qualibration_graphs/superconducting/calibration_utils/two_qubit_interleaved_rb/analysis.py +++ /dev/null @@ -1,170 +0,0 @@ -"""Analysis utilities for two-qubit randomized benchmarking experiments. - -This module provides functions for processing and analyzing raw RB data, -including dataset processing and result logging. -""" - -import logging -from dataclasses import dataclass -from typing import Dict, Tuple - -import numpy as np -import xarray as xr -from qualibrate import QualibrationNode - -# @dataclass -# class FitResults: -# """Stores the relevant fit parameters for a single qubit pair in an RB experiment""" - -# optimal_amplitude: float -# success: bool - - -def log_fitted_results(fit_results: Dict[str, float], log_callable=None): - """ - Logs the node-specific fitted results for all qubit pairs. - - Parameters: - ----------- - fit_results : Dict[str, float] - Dictionary containing floats for each qubit pair. - log_callable : callable, optional - Logger for logging the fitted results. If None, a default logger is used. - """ - if log_callable is None: - log_callable = logging.getLogger(__name__).info - - for qp_name, fit_result in fit_results.items(): - s_qubit = f"Results for qubit pair {qp_name}: " - - s_alpha = f"\tFitted alpha: {fit_result['alpha']:.6f} a.u." - s_fidelity = f"\tFitted fidelity: {100*fit_result['fidelity']:.6f} %" - - if fit_result["success"]: - s_qubit += "SUCCESS!\n" - else: - s_qubit += "FAIL!\n" - - log_message = s_qubit + s_alpha + s_fidelity - - log_callable(log_message) - - -def process_raw_dataset(ds: xr.Dataset, node: QualibrationNode): - """ - Process the raw dataset by adding amplitude and detuning coordinates. - - Parameters: - ----------- - ds : xr.Dataset - Raw dataset from the experiment - node : QualibrationNode - The calibration node containing qubit pairs information - - Returns: - -------- - xr.Dataset - Processed dataset with additional coordinates - """ - ds = node.results["ds_raw"] - - rename_map = {"shots": "average", "sequence": "repeat", "depths": "circuit_depth"} - rename_map = {k: v for k, v in rename_map.items() if k in ds.dims} - - # Assume ds is your input dataset and ds['state'] is your DataArray - state = ds["state"] # shape: (qubit, shots, sequence, depths) - - # Outcome labels for 2-qubit states - labels = ["00", "01", "10", "11"] - - # Create a list of DataArrays: one for each outcome - probs = [state == i for i in range(4)] - - # Stack along a new outcome dimension - probs = xr.concat(probs, dim="outcome") - - # Assign outcome labels - probs = probs.assign_coords(outcome=("outcome", labels)) - - probs_00 = probs.sel(outcome="00") - if rename_map: - probs_00 = probs_00.rename(rename_map) - probs_00 = probs_00.transpose("qubit_pair", "repeat", "circuit_depth", "average") - - probs_00 = probs_00.astype(int) - - if rename_map: - ds_transposed = ds.rename(rename_map) - else: - ds_transposed = ds - ds_transposed = ds_transposed.transpose("qubit_pair", "repeat", "circuit_depth", "average") - - return ds_transposed - - -# def fit_raw_data(ds: xr.Dataset, node: QualibrationNode) -> Tuple[xr.Dataset, Dict[str, FitResults]]: -# """ -# Fit the CZ conditional phase data for each qubit pair. - -# Parameters: -# ----------- -# ds : xr.Dataset -# Dataset containing the processed data. -# node : QualibrationNode -# The calibration node containing parameters and qubit pairs. - -# Returns: -# -------- -# Tuple[xr.Dataset, Dict[str, FitResults]] -# Dataset with fit results and dictionary of fit results for each qubit pair. -# """ -# # For RB analysis, no fitting routine is currently implemented. -# # Pass the dataset through unchanged, or implement RB-specific fitting here if needed. -# ds_fit = ds - -# # Extract the relevant fitted parameters -# ds_fit, fit_results = _extract_relevant_parameters(ds_fit, node) - -# return ds_fit, fit_results - - -# def _extract_relevant_parameters( -# ds_fit: xr.Dataset, node: QualibrationNode -# ) -> Tuple[xr.Dataset, Dict[str, FitResults]]: -# """ -# Extract relevant fit parameters and create FitResults for each qubit pair. - -# Parameters: -# ----------- -# ds_fit : xr.Dataset -# Dataset containing the fit results from fit_routine. -# node : QualibrationNode -# The calibration node containing parameters and qubit pairs. - -# Returns: -# -------- -# Tuple[xr.Dataset, Dict[str, FitResults]] -# Dataset with additional metadata and dictionary of FitResults for each qubit pair. -# """ -# qubit_pairs = node.namespace["qubit_pairs"] - -# # Add metadata attributes to the dataset -# if "optimal_amplitude" in ds_fit.data_vars: -# ds_fit.optimal_amplitude.attrs = {"long_name": "optimal CZ amplitude", "units": "a.u."} -# if "phase_diff" in ds_fit.data_vars: -# ds_fit.phase_diff.attrs = {"long_name": "phase difference", "units": "2π"} -# if "fitted_curve" in ds_fit.data_vars: -# ds_fit.fitted_curve.attrs = {"long_name": "fitted tanh curve", "units": "2π"} - -# # Create FitResults for each qubit pair -# fit_results = {} -# for qp in qubit_pairs: -# qp_name = qp.name -# qp_data = ds_fit.sel(qubit_pair=qp_name) - -# fit_results[qp_name] = FitResults( -# optimal_amplitude=float(qp_data.optimal_amplitude.values), -# success=bool(qp_data.success.values), -# ) - -# return ds_fit, fit_results