diff --git a/src/cqedtoolbox/protocols/configs/parameter-manager-opx_config.json b/src/cqedtoolbox/protocols/configs/parameter-manager-opx_config.json index 14d06a1..c735229 100644 --- a/src/cqedtoolbox/protocols/configs/parameter-manager-opx_config.json +++ b/src/cqedtoolbox/protocols/configs/parameter-manager-opx_config.json @@ -35,6 +35,10 @@ "unit": "", "value": null }, + "parameter_manager.q01.chi": { + "unit": "Hz", + "value": 1e6 + }, "parameter_manager.q01.IF": { "unit": "Hz", "value": 124085680.0 @@ -431,10 +435,6 @@ "unit": "", "value": 0.75 }, - "parameter_manager.corrections.t2e.max_echos": { - "unit": "", - "value": 3 - }, "parameter_manager.corrections.t2e.averaging_factor": { "unit": "", "value": 2.0 diff --git a/src/cqedtoolbox/protocols/operations/__init__.py b/src/cqedtoolbox/protocols/operations/__init__.py index 00f9d28..ccf132f 100644 --- a/src/cqedtoolbox/protocols/operations/__init__.py +++ b/src/cqedtoolbox/protocols/operations/__init__.py @@ -1,3 +1,26 @@ +from enum import Enum + +from cqedtoolbox.fitfuncs.resonators import ( + HangerResponseBruno, + ReflectionResponse, + TransmissionResponse, +) + + +class ResonatorGeometry(Enum): + HANGER = "hanger" + REFLECTION = "reflection" + TRANSMISSION = "transmission" + + @property + def fit_cls(self): + if self is ResonatorGeometry.HANGER: + return HangerResponseBruno + if self is ResonatorGeometry.REFLECTION: + return ReflectionResponse + return TransmissionResponse + + from cqedtoolbox.protocols.operations.single_qubit.res_spec import ResonatorSpectroscopy from cqedtoolbox.protocols.operations.single_qubit.res_spec_vs_gain import ResonatorSpectroscopyVsGain from cqedtoolbox.protocols.operations.fluxonium.res_spec_vs_flux import ResonatorSpectroscopyVsFlux @@ -8,4 +31,4 @@ from cqedtoolbox.protocols.operations.single_qubit.sat_spec import SaturationSpectroscopy from cqedtoolbox.protocols.operations.single_qubit.t1 import T1Operation from cqedtoolbox.protocols.operations.single_qubit.t2e import T2EOperation -from cqedtoolbox.protocols.operations.single_qubit.t2r import T2ROperation \ No newline at end of file +from cqedtoolbox.protocols.operations.single_qubit.t2r import T2ROperation diff --git a/src/cqedtoolbox/protocols/operations/single_qubit/pi_spec.py b/src/cqedtoolbox/protocols/operations/single_qubit/pi_spec.py index a597770..d6c3f04 100644 --- a/src/cqedtoolbox/protocols/operations/single_qubit/pi_spec.py +++ b/src/cqedtoolbox/protocols/operations/single_qubit/pi_spec.py @@ -28,6 +28,7 @@ ) from cqedtoolbox.measurement_lib.opx.advanced.qubit_tuneup import measure_pi_spec from cqedtoolbox.measurement_lib.qick.single_transmon_v2 import PiSpecProgram +from cqedtoolbox.readout.qubit_readout import rotate_complex_qubit_data logger = logging.getLogger(__name__) @@ -166,22 +167,14 @@ def __init__(self, params): self.repetitions, self.averaging_factor, self.max_averaging_increases ) self._register_check("quality_check", self._check_quality, self._increase_averaging) - self._register_success_update(self.qubit_freq, lambda: self.winner_fit.params["x0"].value) + self._register_success_update(self.qubit_freq, lambda: self.fit_result.params["x0"].value) self.independents = {"frequencies": []} self.dependents = {"signal": []} - self.fit_result_re = None - self.fit_result_imag = None - self.fit_result_mag = None - self.snr_re = None - self.snr_imag = None - self.snr_mag = None - self.winner_name = None - self.winner_snr = None - self.winner_fit = None - self.winner_key = None - self.sorted_components = None + self.fit_result = None + self.residuals = None + self.snr = None def _measure_dummy(self) -> Path: logger.info("Starting dummy pi spectroscopy measurement") @@ -198,12 +191,10 @@ def generate(frequencies): return loc def _load_data_dummy(self): - path = self.data_loc / "data.ddh5" - if not path.exists(): - raise FileNotFoundError(f"File {path} does not exist") - data = datadict_from_hdf5(path) - self.independents["frequencies"] = data["frequencies"]["values"] - self.dependents["signal"] = data["signal"]["values"] + data = load_as_xr(self.data_loc) + rotated = rotate_complex_qubit_data(data)[0] + self.independents["frequencies"] = rotated["frequencies"].values + self.dependents["signal"] = rotated["signal"].values def _measure_qick(self) -> Path: logger.info("Starting qick pi spectroscopy measurement") @@ -222,145 +213,65 @@ def _measure_opx(self) -> Path: return loc def _load_data_qick(self): - path = self.data_loc / "data.ddh5" - if not path.exists(): - raise FileNotFoundError(f"File {path} does not exist") - data = datadict_from_hdf5(path) - - self.independents["frequencies"] = data["freq"]["values"] - self.dependents["signal"] = data["signal"]["values"] + data = load_as_xr(self.data_loc) + rotated = rotate_complex_qubit_data(data)[0] + self.independents["frequencies"] = rotated["freq"].values + self.dependents["signal"] = rotated["signal"].values def _load_data_opx(self): - data = load_as_xr(self.data_loc).mean("repetition") + data = load_as_xr(self.data_loc) + if "repetition" in data.dims: + data = data.mean("repetition") + data, _ = rotate_complex_qubit_data(data) self.independents["frequencies"] = data["ssb_frequency"].values - self.dependents["signal"] = data["signal_Re"].values + 1j * data["signal_Im"].values - - def _fit_gaussian_components(self, frequencies, signal, fig_title="") -> tuple: - """ - Fit real, imaginary, and magnitude components with Gaussian fits. - Returns (fit_result_re, fit_result_imag, fit_result_mag, fig_re, fig_imag, fig_mag) - """ - signal_re = signal.real - signal_imag = signal.imag - signal_mag = np.abs(signal) - - # Fit real part - fit_re = Gaussian(frequencies, signal_re) - fit_result_re = fit_re.run(fit_re) - fit_curve_re = fit_result_re.eval() - residuals_re = signal_re - fit_curve_re - amp_re = fit_result_re.params["A"].value - noise_re = np.std(residuals_re) - snr_re = np.abs(amp_re / (4 * noise_re)) - - # Fit imaginary part - fit_imag = Gaussian(frequencies, signal_imag) - fit_result_imag = fit_imag.run(fit_imag) - fit_curve_imag = fit_result_imag.eval() - residuals_imag = signal_imag - fit_curve_imag - amp_imag = fit_result_imag.params["A"].value - noise_imag = np.std(residuals_imag) - snr_imag = np.abs(amp_imag / (4 * noise_imag)) - - # Fit magnitude - fit_mag = Gaussian(frequencies, signal_mag) - fit_result_mag = fit_mag.run(fit_mag) - fit_curve_mag = fit_result_mag.eval() - residuals_mag = signal_mag - fit_curve_mag - amp_mag = fit_result_mag.params["A"].value - noise_mag = np.std(residuals_mag) - snr_mag = np.abs(amp_mag / (4 * noise_mag)) - - # Create three separate figures - # Real plot - fig_re, ax_re = plt.subplots() - ax_re.set_title(f"{fig_title} - Real") - ax_re.set_xlabel("Frequency (MHz)") - ax_re.set_ylabel("Signal Real (A.U)") - ax_re.plot(frequencies, signal_re, label="Data") - ax_re.plot(frequencies, fit_curve_re, label="Fit") - ax_re.legend() - - # Imaginary plot - fig_imag, ax_imag = plt.subplots() - ax_imag.set_title(f"{fig_title} - Imaginary") - ax_imag.set_xlabel("Frequency (MHz)") - ax_imag.set_ylabel("Signal Imaginary (A.U)") - ax_imag.plot(frequencies, signal_imag, label="Data") - ax_imag.plot(frequencies, fit_curve_imag, label="Fit") - ax_imag.legend() - - # Magnitude plot - fig_mag, ax_mag = plt.subplots() - ax_mag.set_title(f"{fig_title} - Magnitude") - ax_mag.set_xlabel("Frequency (MHz)") - ax_mag.set_ylabel("Signal Magnitude (A.U)") - ax_mag.plot(frequencies, signal_mag, label="Data") - ax_mag.plot(frequencies, fit_curve_mag, label="Fit") - ax_mag.legend() - - return ( - (fit_result_re, residuals_re, snr_re), - (fit_result_imag, residuals_imag, snr_imag), - (fit_result_mag, residuals_mag, snr_mag), - fig_re, - fig_imag, - fig_mag - ) + self.dependents["signal"] = data["signal"].values + + def _fit_gaussian(self, frequencies, signal, fig_title="") -> tuple: + fit = Gaussian(frequencies, signal) + fit_result = fit.run(fit) + fit_curve = fit_result.eval() + residuals = signal - fit_curve + amp = fit_result.params["A"].value + noise = np.std(residuals) + snr = np.abs(amp / (4 * noise)) + + fig, ax = plt.subplots() + ax.set_title(fig_title) + ax.set_xlabel("Frequency (MHz)") + ax.set_ylabel("Rotated Signal (A.U)") + ax.plot(frequencies, signal, label="Data") + ax.plot(frequencies, fit_curve, label="Fit") + ax.legend() + + return fit_result, residuals, snr, fig def analyze(self): with DatasetAnalysis(self.data_loc, self.name) as ds: - result_re, result_imag, result_mag, fig_re, fig_imag, fig_mag = self._fit_gaussian_components( + self.fit_result, self.residuals, self.snr, fig = self._fit_gaussian( self.independents["frequencies"], self.dependents["signal"], "Pi Spectroscopy" ) - self.fit_result_re, residuals_re, self.snr_re = result_re - self.fit_result_imag, residuals_imag, self.snr_imag = result_imag - self.fit_result_mag, residuals_mag, self.snr_mag = result_mag - # Save all fit results ds.add( - fit_result_re=self.fit_result_re, - params_re=serialize_fit_params(self.fit_result_re.params), - snr_re=float(self.snr_re), - fit_result_imag=self.fit_result_imag, - params_imag=serialize_fit_params(self.fit_result_imag.params), - snr_imag=float(self.snr_imag), - fit_result_mag=self.fit_result_mag, - params_mag=serialize_fit_params(self.fit_result_mag.params), - snr_mag=float(self.snr_mag) + fit_result=self.fit_result, + params=serialize_fit_params(self.fit_result.params), + snr=float(self.snr) ) - # Save all three figures separately - ds.add_figure(f"{self.name}_real", fig=fig_re) - image_path_re = ds._new_file_path(ds.savefolders[1], f"{self.name}_real", suffix="png") - self.figure_paths.append(image_path_re) - - ds.add_figure(f"{self.name}_imag", fig=fig_imag) - image_path_imag = ds._new_file_path(ds.savefolders[1], f"{self.name}_imag", suffix="png") - self.figure_paths.append(image_path_imag) - - ds.add_figure(f"{self.name}_mag", fig=fig_mag) - image_path_mag = ds._new_file_path(ds.savefolders[1], f"{self.name}_mag", suffix="png") - self.figure_paths.append(image_path_mag) + ds.add_figure(self.name, fig=fig) + image_path = ds._new_file_path(ds.savefolders[1], self.name, suffix="png") + self.figure_paths.append(image_path) def _check_quality(self) -> CheckResult: - snr_dict = { - "Real": (self.snr_re, self.fit_result_re, "re"), - "Imaginary": (self.snr_imag, self.fit_result_imag, "imag"), - "Magnitude": (self.snr_mag, self.fit_result_mag, "mag"), - } - self.sorted_components = sorted(snr_dict.items(), key=lambda x: x[1][0], reverse=True) - self.winner_name, (self.winner_snr, self.winner_fit, self.winner_key) = self.sorted_components[0] - + # TODO: make sure that the fit frequency is inside the swept range threshold = self.snr_threshold() - snr_passed = self.winner_snr >= threshold + snr_passed = self.snr >= threshold max_error = self.max_fit_param_error() bad_params = [] - for pname, param in self.winner_fit.params.items(): + for pname, param in self.fit_result.params.items(): if param.stderr is None: bad_params.append(f"{pname}(no stderr)") elif param.value != 0 and abs(param.stderr / param.value) > max_error: @@ -368,20 +279,14 @@ def _check_quality(self) -> CheckResult: bad_params.append(f"{pname}({pct:.0f}%)") passed = snr_passed and len(bad_params) == 0 - parts = [f"SNR={self.winner_snr:.3f} (threshold={threshold:.3f}, component={self.winner_name})"] + parts = [f"SNR={self.snr:.3f} (threshold={threshold:.3f})"] if bad_params: parts.append(f"high-error params: {', '.join(bad_params)}") return CheckResult("quality_check", passed, "; ".join(parts)) def correct(self, result: EvaluateResult) -> EvaluateResult: - # Pull all three figures before super() can auto-append the last one. - # figure_paths order after analyze(): [0]=real, [1]=imag, [2]=mag - plot_map = {} - if len(self.figure_paths) >= 3: - plot_map["re"] = self.figure_paths[0].resolve() - plot_map["imag"] = self.figure_paths[1].resolve() - plot_map["mag"] = self.figure_paths[2].resolve() - self.figure_paths.clear() # prevent auto-append + figure = self.figure_paths[0].resolve() if self.figure_paths else None + self.figure_paths.clear() header = (f"## Pi Spectroscopy\n" f"Frequencies: {self.start_freq():.3f}–{self.end_freq():.3f} MHz\n" @@ -390,42 +295,12 @@ def correct(self, result: EvaluateResult) -> EvaluateResult: result = super().correct(result) # adds check table; no auto-figure since list is empty - if self.sorted_components: - winner_name, (winner_snr, winner_fit, winner_key) = self.sorted_components[0] - winner_report = f"**Fit Report:**\n```\n{str(winner_fit.lmfit_result.fit_report())}\n```\n\n" - - if result.status == OperationStatus.SUCCESS: - self.report_output.extend([ - f"### **{winner_name} Component (SELECTED)**\n" - f"Fit was **SUCCESSFUL** with {winner_name} SNR of {winner_snr:.3f}\n" - f"This component was selected because it has the highest SNR.\n\n", - plot_map.get(winner_key, ""), - winner_report, - ]) - for comp_name, (comp_snr, comp_fit, comp_key) in self.sorted_components[1:]: - threshold = self.snr_threshold() - if comp_snr >= threshold: - reason = f"SNR of {comp_snr:.3f} is above threshold but lower than {winner_name} (SNR={winner_snr:.3f})" - else: - reason = f"SNR of {comp_snr:.3f} is below threshold of {threshold:.3f}" - self.report_output.extend([ - f"### **{comp_name} Component (NOT SELECTED)**\n", - plot_map.get(comp_key, ""), - f"Not used: {reason}\n\n**Fit Report:**\n```\n{str(comp_fit.lmfit_result.fit_report())}\n```\n\n", - ]) - else: - self.report_output.extend([ - f"### **{winner_name} Component (Highest SNR)**\n" - f"Fit was **UNSUCCESSFUL** with {winner_name} SNR of {winner_snr:.3f}\n" - f"NO value has been changed.\n\n", - plot_map.get(winner_key, ""), - winner_report, - ]) - for comp_name, (comp_snr, comp_fit, comp_key) in self.sorted_components[1:]: - self.report_output.extend([ - f"### **{comp_name} Component**\n", - plot_map.get(comp_key, ""), - f"SNR of {comp_snr:.3f}\n\n**Fit Report:**\n```\n{str(comp_fit.lmfit_result.fit_report())}\n```\n\n", - ]) + status_line = "SUCCESSFUL" if result.status == OperationStatus.SUCCESS else "UNSUCCESSFUL" + self.report_output.extend([ + f"### Rotated Signal Fit\n" + f"Fit was **{status_line}** with SNR of {self.snr:.3f}\n\n", + figure or "", + f"**Fit Report:**\n```\n{str(self.fit_result.lmfit_result.fit_report())}\n```\n\n", + ]) return result diff --git a/src/cqedtoolbox/protocols/operations/single_qubit/power_rabi.py b/src/cqedtoolbox/protocols/operations/single_qubit/power_rabi.py index 23be70b..3615ebc 100644 --- a/src/cqedtoolbox/protocols/operations/single_qubit/power_rabi.py +++ b/src/cqedtoolbox/protocols/operations/single_qubit/power_rabi.py @@ -27,6 +27,7 @@ ) from cqedtoolbox.measurement_lib.opx.advanced.qubit_tuneup import measure_power_rabi from cqedtoolbox.measurement_lib.qick.single_transmon_v2 import AmplitudeRabiProgram +from cqedtoolbox.readout.qubit_readout import rotate_complex_qubit_data logger = logging.getLogger(__name__) @@ -362,23 +363,15 @@ def __init__(self, params): self._register_success_update( self.qubit_gain, - lambda: 1 / (2 * self._winner_fit.params["f"].value), + lambda: 1 / (2 * self.fit_result.params["f"].value), ) self.independents = {"gains": []} self.dependents = {"signal": []} - self.fit_result_re = None - self.fit_result_imag = None - self.fit_result_mag = None - self.snr_re = None - self.snr_imag = None - self.snr_mag = None - self._winner_fit = None - self._winner_snr = None - self._winner_key = None - self._winner_name = None - self._sorted_components = None + self.fit_result = None + self.residuals = None + self.snr = None def _measure_qick(self) -> Path: logger.info("Starting qick power rabi measurement") @@ -411,151 +404,70 @@ def _measure_dummy(self): return loc def _load_data_qick(self): - path = self.data_loc / "data.ddh5" - if not path.exists(): - raise FileNotFoundError(f"File {path} does not exist") - data = datadict_from_hdf5(path) - - self.independents["gains"] = data["gain"]["values"] - self.dependents["signal"] = data["signal"]["values"] + data = load_as_xr(self.data_loc) + rotated = rotate_complex_qubit_data(data)[0] + self.independents["gains"] = rotated["gain"].values + self.dependents["signal"] = rotated["signal"].values def _load_data_opx(self): - data = load_as_xr(self.data_loc).mean("repetition") + data = load_as_xr(self.data_loc) + if "repetition" in data.dims: + data = data.mean("repetition") + data, _ = rotate_complex_qubit_data(data) self.independents["gains"] = data["amplitude"].values - self.dependents["signal"] = data["signal_Re"].values + 1j * data["signal_Im"].values + self.dependents["signal"] = data["signal"].values def _load_data_dummy(self): - path = self.data_loc / "data.ddh5" - if not path.exists(): - raise FileNotFoundError(f"File {path} does not exist") - data = datadict_from_hdf5(path) - - self.independents["gains"] = data["gains"]["values"] - self.dependents["signal"] = data["signal"]["values"] - - def _fit_cosine_components(self, gains, signal, fig_title="") -> tuple: - """ - Fit real, imaginary, and magnitude components with Cosine fits. - Returns (fit_result_re, fit_result_imag, fit_result_mag, fig_re, fig_imag, fig_mag) - """ - signal_re = signal.real - signal_imag = signal.imag - signal_mag = np.abs(signal) - - # Fit real part - fit_re = Cosine(gains, signal_re) - fit_result_re = fit_re.run(fit_re) - fit_curve_re = fit_result_re.eval() - residuals_re = signal_re - fit_curve_re - amp_re = fit_result_re.params["A"].value - noise_re = np.std(residuals_re) - snr_re = np.abs(amp_re / (4 * noise_re)) - - # Fit imaginary part - fit_imag = Cosine(gains, signal_imag) - fit_result_imag = fit_imag.run(fit_imag) - fit_curve_imag = fit_result_imag.eval() - residuals_imag = signal_imag - fit_curve_imag - amp_imag = fit_result_imag.params["A"].value - noise_imag = np.std(residuals_imag) - snr_imag = np.abs(amp_imag / (4 * noise_imag)) - - # Fit magnitude - fit_mag = Cosine(gains, signal_mag) - fit_result_mag = fit_mag.run(fit_mag) - fit_curve_mag = fit_result_mag.eval() - residuals_mag = signal_mag - fit_curve_mag - amp_mag = fit_result_mag.params["A"].value - noise_mag = np.std(residuals_mag) - snr_mag = np.abs(amp_mag / (4 * noise_mag)) - - # Create three separate figures - fig_re, ax_re = plt.subplots() - ax_re.set_title(f"{fig_title} - Real") - ax_re.set_xlabel("Gain (A.U)") - ax_re.set_ylabel("Signal Real (A.U)") - ax_re.plot(gains, signal_re, label="Data") - ax_re.plot(gains, fit_curve_re, label="Fit") - ax_re.legend() - - fig_imag, ax_imag = plt.subplots() - ax_imag.set_title(f"{fig_title} - Imaginary") - ax_imag.set_xlabel("Gain (A.U)") - ax_imag.set_ylabel("Signal Imaginary (A.U)") - ax_imag.plot(gains, signal_imag, label="Data") - ax_imag.plot(gains, fit_curve_imag, label="Fit") - ax_imag.legend() - - fig_mag, ax_mag = plt.subplots() - ax_mag.set_title(f"{fig_title} - Magnitude") - ax_mag.set_xlabel("Gain (A.U)") - ax_mag.set_ylabel("Signal Magnitude (A.U)") - ax_mag.plot(gains, signal_mag, label="Data") - ax_mag.plot(gains, fit_curve_mag, label="Fit") - ax_mag.legend() - - return ( - (fit_result_re, residuals_re, snr_re), - (fit_result_imag, residuals_imag, snr_imag), - (fit_result_mag, residuals_mag, snr_mag), - fig_re, - fig_imag, - fig_mag - ) + data = load_as_xr(self.data_loc) + rotated = rotate_complex_qubit_data(data)[0] + self.independents["gains"] = rotated["gains"].values + self.dependents["signal"] = rotated["signal"].values + + def _fit_cosine(self, gains, signal, fig_title="") -> tuple: + fit = Cosine(gains, signal) + fit_result = fit.run(fit) + fit_curve = fit_result.eval() + residuals = signal - fit_curve + amp = fit_result.params["A"].value + noise = np.std(residuals) + snr = np.abs(amp / (4 * noise)) + + fig, ax = plt.subplots() + ax.set_title(fig_title) + ax.set_xlabel("Gain (A.U)") + ax.set_ylabel("Rotated Signal (A.U)") + ax.plot(gains, signal, label="Data") + ax.plot(gains, fit_curve, label="Fit") + ax.legend() + + return fit_result, residuals, snr, fig def analyze(self): with DatasetAnalysis(self.data_loc, self.name) as ds: - result_re, result_imag, result_mag, fig_re, fig_imag, fig_mag = self._fit_cosine_components( + self.fit_result, self.residuals, self.snr, fig = self._fit_cosine( self.independents["gains"], self.dependents["signal"], "Power Rabi" ) - self.fit_result_re, residuals_re, self.snr_re = result_re - self.fit_result_imag, residuals_imag, self.snr_imag = result_imag - self.fit_result_mag, residuals_mag, self.snr_mag = result_mag - - # Determine winner (highest SNR) for check and success update - snr_dict = { - "Real": (self.snr_re, self.fit_result_re, "re"), - "Imaginary": (self.snr_imag, self.fit_result_imag, "imag"), - "Magnitude": (self.snr_mag, self.fit_result_mag, "mag"), - } - self._sorted_components = sorted(snr_dict.items(), key=lambda x: x[1][0], reverse=True) - self._winner_name, (self._winner_snr, self._winner_fit, self._winner_key) = self._sorted_components[0] - # Save all fit results ds.add( - fit_result_re=self.fit_result_re, - params_re=serialize_fit_params(self.fit_result_re.params), - snr_re=float(self.snr_re), - fit_result_imag=self.fit_result_imag, - params_imag=serialize_fit_params(self.fit_result_imag.params), - snr_imag=float(self.snr_imag), - fit_result_mag=self.fit_result_mag, - params_mag=serialize_fit_params(self.fit_result_mag.params), - snr_mag=float(self.snr_mag) + fit_result=self.fit_result, + params=serialize_fit_params(self.fit_result.params), + snr=float(self.snr) ) - ds.add_figure(f"{self.name}_real", fig=fig_re) - image_path_re = ds._new_file_path(ds.savefolders[1], f"{self.name}_real", suffix="png") - self.figure_paths.append(image_path_re) - - ds.add_figure(f"{self.name}_imag", fig=fig_imag) - image_path_imag = ds._new_file_path(ds.savefolders[1], f"{self.name}_imag", suffix="png") - self.figure_paths.append(image_path_imag) - - ds.add_figure(f"{self.name}_mag", fig=fig_mag) - image_path_mag = ds._new_file_path(ds.savefolders[1], f"{self.name}_mag", suffix="png") - self.figure_paths.append(image_path_mag) + ds.add_figure(self.name, fig=fig) + image_path = ds._new_file_path(ds.savefolders[1], self.name, suffix="png") + self.figure_paths.append(image_path) def _check_quality(self) -> CheckResult: threshold = self.snr_threshold() - snr_passed = self._winner_snr >= threshold + snr_passed = self.snr >= threshold max_error = self.max_fit_param_error() bad_params = [] - for pname, param in self._winner_fit.params.items(): + for pname, param in self.fit_result.params.items(): if param.stderr is None: bad_params.append(f"{pname}(no stderr)") elif param.value == 0 or abs(param.stderr / param.value) > max_error: @@ -563,19 +475,14 @@ def _check_quality(self) -> CheckResult: bad_params.append(f"{pname}({pct:.0f}%)") passed = snr_passed and len(bad_params) == 0 - parts = [f"winner={self._winner_name}, SNR={self._winner_snr:.3f} (threshold={threshold:.3f})"] + parts = [f"SNR={self.snr:.3f} (threshold={threshold:.3f})"] if bad_params: parts.append(f"high-error params: {', '.join(bad_params)}") return CheckResult("quality_check", passed, "; ".join(parts)) def correct(self, result: EvaluateResult) -> EvaluateResult: - # Pop all figures before super() auto-appends the last one - fig_re = self.figure_paths.pop(0) if len(self.figure_paths) >= 3 else None - fig_imag = self.figure_paths.pop(0) if self.figure_paths else None - fig_mag = self.figure_paths.pop(0) if self.figure_paths else None - self.figure_paths.clear() # prevent auto-append - - plot_map = {"re": fig_re, "imag": fig_imag, "mag": fig_mag} + figure = self.figure_paths[0] if self.figure_paths else None + self.figure_paths.clear() self.report_output.append( f"## Power Rabi\n" @@ -584,15 +491,13 @@ def correct(self, result: EvaluateResult) -> EvaluateResult: f"Data Path: `{self.data_loc}`\n\n" ) - for i, (comp_name, (comp_snr, comp_fit, comp_key)) in enumerate(self._sorted_components): - tag = "(SELECTED)" if i == 0 else "(NOT SELECTED)" - self.report_output.append(f"### **{comp_name} Component {tag}**\n") - if plot_map[comp_key]: - self.report_output.append(plot_map[comp_key]) - self.report_output.append( - f"SNR={comp_snr:.3f}\n\n" - f"**Fit Report:**\n```\n{str(comp_fit.lmfit_result.fit_report())}\n```\n\n" - ) + self.report_output.append("### Rotated Signal Fit\n") + if figure: + self.report_output.append(figure) + self.report_output.append( + f"SNR={self.snr:.3f}\n\n" + f"**Fit Report:**\n```\n{str(self.fit_result.lmfit_result.fit_report())}\n```\n\n" + ) result = super().correct(result) # adds check table + success update line return result diff --git a/src/cqedtoolbox/protocols/operations/single_qubit/res_spec.py b/src/cqedtoolbox/protocols/operations/single_qubit/res_spec.py index 5b5b56b..e3b6ffe 100644 --- a/src/cqedtoolbox/protocols/operations/single_qubit/res_spec.py +++ b/src/cqedtoolbox/protocols/operations/single_qubit/res_spec.py @@ -6,20 +6,25 @@ from numpy.typing import ArrayLike import matplotlib.pyplot as plt +from scipy.signal import savgol_filter +from scipy.interpolate import CubicSpline +from scipy.optimize import curve_fit + from labcore.analysis import DatasetAnalysis, FitResult from labcore.measurement.storage import run_and_save_sweep from labcore.data.datadict_storage import datadict_from_hdf5, load_as_xr from labcore.measurement import sweep_parameter, record_as from labcore.protocols.base import (ProtocolOperation, OperationStatus, serialize_fit_params, - ParamImprovement, CorrectionParameter, CheckResult, Correction) + ParamImprovement, CorrectionParameter, CheckResult, Correction, PlatformTypes) +from cqedtoolbox.protocols.operations import ResonatorGeometry from cqedtoolbox.protocols.parameters import (Repetition, ResonatorSpecSteps, ReadoutGain, ReadoutLength, StartReadoutFrequency, EndReadoutFrequency, ReadoutFrequency, nestedAttributeFromString) from cqedtoolbox.measurement_lib.opx.advanced.qubit_tuneup import measure_pulse_resonator_spec from cqedtoolbox.measurement_lib.qick.single_transmon_v2 import FreqSweepProgram -from cqedtoolbox.fitfuncs.resonators import HangerResponseBruno +from cqedtoolbox.fitfuncs.resonators import HangerResponseBruno, ReflectionResponse, TransmissionResponse logger = logging.getLogger(__name__) @@ -33,7 +38,9 @@ class UnwindAndFitRet: fit_curve: ArrayLike fit_result: FitResult residuals: ArrayLike + residual_std: float snr: float + unwind_sign: int fig: plt.Figure ax: plt.Axes @@ -286,22 +293,94 @@ def apply(self) -> None: def report_output(self) -> str: return self._last_change + +def fit_sine(t, data): + """Fit a sine wave to the given data.""" + t = np.asarray(t, dtype=float) + data = np.asarray(data, dtype=float) -class ResonatorSpectroscopy(ProtocolOperation): + def sine_function(x, amplitude, frequency, phase, offset): + return amplitude * np.sin(2 * np.pi * frequency * x + phase) + offset + + y_offset = np.mean(data) + amplitude = (np.max(data) - np.min(data)) / 2 + + fft = np.fft.fft(data - y_offset) + freqs = np.fft.fftfreq(len(data), d=t[1] - t[0]) + positive = np.argmax(np.abs(fft[1:len(fft) // 2 + 1])) + 1 if len(fft) > 2 else 0 + frequency = np.abs(freqs[positive]) + phase = np.angle(fft[positive]) if positive else 0.0 + + initial_guess = (amplitude, frequency, phase, y_offset) + params, covariance = curve_fit(sine_function, t, data, p0=initial_guess, maxfev=5000) + return params, covariance + +def background_filter(x, y): + window_size = 8 + + def moving_window_variance(data): + half_window = window_size // 2 + padded_data = np.pad(data, (half_window, half_window), mode="reflect") + return np.array([np.var(padded_data[i:i + window_size]) for i in range(len(data))]) + + def filtered_variance(x, y): + smooth = savgol_filter(y, window_size, 1) + spline_derivative = CubicSpline(x, smooth)(y, 1) + return savgol_filter(moving_window_variance(spline_derivative), window_size, 1) + + s = filtered_variance(x, y.real) + 1j * filtered_variance(x, y.imag) + sabs = np.abs(s) + s0 = np.median(sabs) + sm, sp = s0 - 0.5 * sabs.std(), s0 + 0.5 * sabs.std() + return (sabs > sm) & (sabs < sp) +def unwind_signal(x, y, f=None, sign=1): + """Fit and remove the dominant sinusoidal loading from the complex signal.""" + x = np.asarray(x, dtype=float) + y = np.asarray(y) + + if f is None: + bg = background_filter(x, y) + if np.count_nonzero(bg) < 5: + bg = np.ones(len(x), dtype=bool) + + xbg, ybg = x[bg], y[bg] + ixbg = np.linspace(x[0], x[-1], x.size) + iybg = np.interp(ixbg, xbg, ybg) + + pr, _ = fit_sine(ixbg, iybg.real) + pi, _ = fit_sine(ixbg, iybg.imag) + f = np.mean([pr[1], pi[1]]) + + unwound = y * np.exp(1j * sign * 2 * np.pi * f * x) + return unwound.real, unwound.imag, f + + +class ResonatorSpectroscopy(ProtocolOperation): _SIM_F0 = 7e9 _SIM_QI = 20e3 _SIM_QC = 20e3 _SIM_A = 4.0 _SIM_PHI = 0.0 _SIM_NOISE_AMP = 0.05 - - def __init__(self, params): + def __init__(self, params, geometry: ResonatorGeometry | str): super().__init__() self.params = params + if isinstance(geometry, str): + try: + geometry = ResonatorGeometry(geometry.lower()) + except ValueError as err: + valid = ", ".join(g.value for g in ResonatorGeometry) + raise ValueError( + f"Unsupported resonator geometry '{geometry}'. Expected one of: {valid}" + ) from err + + self.geometry = geometry + self._fit_cls = geometry.fit_cls + self._register_inputs( repetitions=Repetition(params), steps=ResonatorSpecSteps(params), @@ -349,6 +428,10 @@ def __init__(self, params): self._check_quality, [self._window_shift, self._increase_sampling, self._increase_averaging], ) + self._register_check( + "fit_in_range", + self._check_fit_in_range, + ) self._register_success_update( self.readout_freq, @@ -357,12 +440,12 @@ def __init__(self, params): self._register_success_update( self.start_frequency, - lambda: self.fit_result.params["f_0"].value - 5, + lambda: self.fit_result.params["f_0"].value - 5 if self.platform_type == PlatformTypes.QICK else self.fit_result.params["f_0"].value - 5e6, ) self._register_success_update( self.end_frequency, - lambda: self.fit_result.params["f_0"].value + 5, + lambda: self.fit_result.params["f_0"].value + 5 if self.platform_type == PlatformTypes.QICK else self.fit_result.params["f_0"].value + 5e6, ) self.condition = f"Success if the SNR of the measurement is bigger than the current threshold of " # {self.SNR_THRESHOLD}" @@ -410,23 +493,58 @@ def _measure_dummy(self): logger.info("Dummy measurement complete") return loc - @staticmethod - def add_mag_and_unwind_and_fit(frequencies, signal_raw, fig_title="") -> UnwindAndFitRet: - phase_unwrap = np.unwrap(np.angle(signal_raw)) - phase_slope = np.polyfit(frequencies, phase_unwrap, 1)[0] + def _load_data_qick(self): + path = self.data_loc/"data.ddh5" + if not path.exists(): + raise FileNotFoundError(f"File {path} does not exist") + data = datadict_from_hdf5(path) - signal_unwind = signal_raw * np.exp(-1j * frequencies * phase_slope) - magnitude = np.abs(signal_raw) - phase = np.arctan2(signal_unwind.imag, signal_unwind.real) + self.independents["frequencies"] = data["freq"]["values"] + self.dependents["signal"] = data["signal"]["values"] + + def _load_data_opx(self): + data = load_as_xr(self.data_loc).mean("repetition") + q = nestedAttributeFromString(self.params, "active.qubit")() + lo = nestedAttributeFromString(self.params, f"{q}.readout.LO")() + self.independents["frequencies"] = data["ssb_frequency"].values + lo + self.dependents["signal"] = data["signal_Re"].values + 1j * data["signal_Im"].values + + def _load_data_dummy(self): + path = self.data_loc/"data.ddh5" + if not path.exists(): + raise FileNotFoundError(f"File {path} does not exist") + data = datadict_from_hdf5(path) - fit = HangerResponseBruno(frequencies, signal_unwind) - fit_result = fit.run(fit) - fit_curve = fit_result.eval() - residuals = signal_unwind - fit_curve + self.independents["frequencies"] = data["frequencies"]["values"] + self.dependents["signal"] = data["signal"]["values"] - amp = fit_result.params["A"].value - noise = np.std(residuals) - snr = np.abs(amp / (4 * noise)) + @staticmethod + def add_mag_and_unwind_and_fit(frequencies, signal_raw, platform_type, fit_cls, fig_title="") -> UnwindAndFitRet: + frequencies = np.asarray(frequencies, dtype=float) + signal_raw = np.asarray(signal_raw) + + magnitude = np.abs(signal_raw) + del platform_type + + def fit_candidate(sign: int): + unwound_real, unwound_imag, _ = unwind_signal( + frequencies, signal_raw, sign=sign + ) + signal_unwind = unwound_real + 1j * unwound_imag + fit = fit_cls(frequencies, signal_unwind) + fit_result = fit.run(fit) + fit_curve = fit_result.eval() + residuals = signal_unwind - fit_curve + residual_std = float(np.std(residuals)) + amp = fit_result.params["A"].value + snr = np.abs(amp / (4 * residual_std)) if residual_std > 0 else float("inf") + return signal_unwind, fit_result, fit_curve, residuals, residual_std, snr + + candidates = [fit_candidate(sign) for sign in (1, -1)] + best_idx = min(range(len(candidates)), key=lambda i: candidates[i][4]) + signal_unwind, fit_result, fit_curve, residuals, residual_std, snr = candidates[best_idx] + unwind_sign = 1 if best_idx == 0 else -1 + phase = np.angle(signal_unwind) fig, ax = plt.subplots() ax.set_title(fig_title) @@ -443,44 +561,24 @@ def add_mag_and_unwind_and_fit(frequencies, signal_raw, fig_title="") -> UnwindA fit_curve=fit_curve, fit_result=fit_result, residuals=residuals, + residual_std=residual_std, snr=snr, + unwind_sign=unwind_sign, fig=fig, ax=ax, ) return ret - def _load_data_qick(self): - path = self.data_loc/"data.ddh5" - if not path.exists(): - raise FileNotFoundError(f"File {path} does not exist") - data = datadict_from_hdf5(path) - - self.independents["frequencies"] = data["freq"]["values"] - self.dependents["signal"] = data["signal"]["values"] - - def _load_data_opx(self): - data = load_as_xr(self.data_loc).mean("repetition") - q = nestedAttributeFromString(self.params, "active.qubit")() - lo = nestedAttributeFromString(self.params, f"{q}.readout.LO")() - self.independents["frequencies"] = data["ssb_frequency"].values + lo - self.dependents["signal"] = data["signal_Re"].values + 1j * data["signal_Im"].values - - def _load_data_dummy(self): - path = self.data_loc/"data.ddh5" - if not path.exists(): - raise FileNotFoundError(f"File {path} does not exist") - data = datadict_from_hdf5(path) - - self.independents["frequencies"] = data["frequencies"]["values"] - self.dependents["signal"] = data["signal"]["values"] - def analyze(self): with DatasetAnalysis(self.data_loc, self.name) as ds: ret = self.add_mag_and_unwind_and_fit(self.independents["frequencies"], self.dependents["signal"], + self.platform_type, + self._fit_cls, "Resonator Spectroscopy") + self.unwind_signal = ret.signal_unwind self.magnitude = ret.magnitude self.phase = ret.phase self.snr = ret.snr @@ -521,3 +619,19 @@ def _check_quality(self) -> CheckResult: parts.append(f"high-error params: {', '.join(bad_params)}") return CheckResult("quality_check", passed, "; ".join(parts)) + + def _check_fit_in_range(self) -> CheckResult: + if self.fit_result is None: + raise RuntimeError("Fit result must be set before checking fit range") + + freqs = np.asarray(self.independents["frequencies"], dtype=float) + f0 = self.fit_result.params["f_0"].value + fmin = float(np.min(freqs)) + fmax = float(np.max(freqs)) + passed = fmin <= f0 <= fmax + + return CheckResult( + "fit_in_range", + passed, + f"f_0={f0:.6g}, sweep=[{fmin:.6g}, {fmax:.6g}]", + ) diff --git a/src/cqedtoolbox/protocols/operations/single_qubit/res_spec_after_pi.py b/src/cqedtoolbox/protocols/operations/single_qubit/res_spec_after_pi.py index 89f9b1c..334fd43 100644 --- a/src/cqedtoolbox/protocols/operations/single_qubit/res_spec_after_pi.py +++ b/src/cqedtoolbox/protocols/operations/single_qubit/res_spec_after_pi.py @@ -15,6 +15,7 @@ from labcore.protocols.base import (ProtocolOperation, OperationStatus, serialize_fit_params, CorrectionParameter, CheckResult, EvaluateResult) +from cqedtoolbox.protocols.operations import ResonatorGeometry from cqedtoolbox.protocols.parameters import ( Repetition, ResonatorSpecSteps, @@ -28,8 +29,10 @@ from cqedtoolbox.measurement_lib.opx.advanced.qubit_tuneup import measure_pulse_resonator_spec_after_pi_pulse from cqedtoolbox.measurement_lib.qick.single_transmon_v2 import FreqSweepProgram, ResProbeProgram -from cqedtoolbox.fitfuncs.resonators import HangerResponseBruno -from cqedtoolbox.protocols.operations.single_qubit.res_spec import SyntheticHangerResonatorData +from cqedtoolbox.protocols.operations.single_qubit.res_spec import ( + ResonatorSpectroscopy, + SyntheticHangerResonatorData, +) logger = logging.getLogger(__name__) @@ -92,11 +95,22 @@ def _opx_setter(self, value): class ResonatorSpectroscopyAfterPi(ProtocolOperation): - - def __init__(self, params): + def __init__(self, params, geometry: ResonatorGeometry | str): super().__init__() self.params = params + if isinstance(geometry, str): + try: + geometry = ResonatorGeometry(geometry.lower()) + except ValueError as err: + valid = ", ".join(g.value for g in ResonatorGeometry) + raise ValueError( + f"Unsupported resonator geometry '{geometry}'. Expected one of: {valid}" + ) from err + + self.geometry = geometry + self._fit_cls = geometry.fit_cls + self._register_inputs( repetitions=Repetition(params), steps=ResonatorSpecSteps(params), @@ -117,6 +131,8 @@ def __init__(self, params): self._register_check("quality_check_before", self._check_quality_before, None) self._register_check("quality_check_after", self._check_quality_after, None) + self._register_check("fit_in_range_before", self._check_fit_in_range_before, None) + self._register_check("fit_in_range_after", self._check_fit_in_range_after, None) self._register_check("detuning_check", self._check_detuning, None) self._register_success_update(self.detuning, lambda: self.chi) @@ -215,34 +231,6 @@ def _measure_opx(self) -> Path: logger.info("Measurement complete") return loc - def _add_mag_and_unwind_and_fit(self, frequencies, signal_raw, fig_title="") -> tuple: - """Unwind phase, calculate magnitude, and fit with hanger response""" - phase_unwrap = np.unwrap(np.angle(signal_raw)) - phase_slope = np.polyfit(frequencies, phase_unwrap, 1)[0] - - signal_unwind = signal_raw * np.exp(-1j * frequencies * phase_slope) - magnitude = np.abs(signal_raw) - phase = np.arctan2(signal_unwind.imag, signal_unwind.real) - - fit = HangerResponseBruno(frequencies, signal_unwind) - fit_result = fit.run(fit) - fit_curve = fit_result.eval() - residuals = signal_unwind - fit_curve - - amp = fit_result.params["A"].value - noise = np.std(residuals) - snr = np.abs(amp / (4 * noise)) - - fig, ax = plt.subplots() - ax.set_title(fig_title) - ax.set_xlabel("Frequency (MHz)") - ax.set_ylabel("Magnitude Signal (A.U)") - ax.plot(frequencies, magnitude, label="Data") - ax.plot(frequencies, np.abs(fit_curve), label="Fit") - ax.legend() - - return signal_unwind, magnitude, phase, fit_result, residuals, snr, fig - def _load_data_qick(self): # Load before data path_before = self.data_loc_before / "data.ddh5" @@ -278,48 +266,56 @@ def _load_data_opx(self): def analyze(self): # Analyze before measurement with DatasetAnalysis(self.data_loc_before, f"{self.name}_before") as ds: - ret_before = self._add_mag_and_unwind_and_fit( + ret_before = ResonatorSpectroscopy.add_mag_and_unwind_and_fit( self.independents_before["frequencies"], self.dependents_before["signal"], + self.platform_type, + self._fit_cls, "Resonator Spectroscopy Before Pi" ) - (self.unwind_signal_before, self.magnitude_before, phase_before, - self.fit_result_before, residuals_before, self.snr_before, fig_before) = ret_before + self.unwind_signal_before = ret_before.signal_unwind + self.magnitude_before = ret_before.magnitude + self.fit_result_before = ret_before.fit_result + self.snr_before = ret_before.snr self.f0_before = self.fit_result_before.params["f_0"].value ds.add( - fit_curve=self.fit_result_before.eval(), + fit_curve=ret_before.fit_curve, fit_result=self.fit_result_before, params=serialize_fit_params(self.fit_result_before.params), snr=float(self.snr_before) ) - ds.add_figure(f"{self.name}_before", fig=fig_before) + ds.add_figure(f"{self.name}_before", fig=ret_before.fig) image_path_before = ds._new_file_path(ds.savefolders[1], f"{self.name}_before", suffix="png") self.figure_paths.append(image_path_before) # Analyze after measurement with DatasetAnalysis(self.data_loc_after, f"{self.name}_after") as ds: - ret_after = self._add_mag_and_unwind_and_fit( + ret_after = ResonatorSpectroscopy.add_mag_and_unwind_and_fit( self.independents_after["frequencies"], self.dependents_after["signal"], + self.platform_type, + self._fit_cls, "Resonator Spectroscopy After Pi" ) - (self.unwind_signal_after, self.magnitude_after, phase_after, - self.fit_result_after, residuals_after, self.snr_after, fig_after) = ret_after + self.unwind_signal_after = ret_after.signal_unwind + self.magnitude_after = ret_after.magnitude + self.fit_result_after = ret_after.fit_result + self.snr_after = ret_after.snr self.f0_after = self.fit_result_after.params["f_0"].value ds.add( - fit_curve=self.fit_result_after.eval(), + fit_curve=ret_after.fit_curve, fit_result=self.fit_result_after, params=serialize_fit_params(self.fit_result_after.params), snr=float(self.snr_after) ) - ds.add_figure(f"{self.name}_after", fig=fig_after) + ds.add_figure(f"{self.name}_after", fig=ret_after.fig) image_path_after = ds._new_file_path(ds.savefolders[1], f"{self.name}_after", suffix="png") self.figure_paths.append(image_path_after) @@ -392,6 +388,28 @@ def _check_quality_before(self) -> CheckResult: def _check_quality_after(self) -> CheckResult: return self._check_fit_quality(self.snr_after, self.fit_result_after, "quality_check_after") + def _check_fit_in_range(self, freqs, fit_result, check_name) -> CheckResult: + freqs = np.asarray(freqs, dtype=float) + f0 = fit_result.params["f_0"].value + fmin = float(np.min(freqs)) + fmax = float(np.max(freqs)) + passed = fmin <= f0 <= fmax + return CheckResult(check_name, passed, f"f_0={f0:.6g}, sweep=[{fmin:.6g}, {fmax:.6g}]") + + def _check_fit_in_range_before(self) -> CheckResult: + return self._check_fit_in_range( + self.independents_before["frequencies"], + self.fit_result_before, + "fit_in_range_before", + ) + + def _check_fit_in_range_after(self) -> CheckResult: + return self._check_fit_in_range( + self.independents_after["frequencies"], + self.fit_result_after, + "fit_in_range_after", + ) + def _check_detuning(self) -> CheckResult: threshold = self.detuning_threshold() passed = abs(self.chi) >= threshold @@ -416,7 +434,8 @@ def correct(self, result: EvaluateResult) -> EvaluateResult: self.report_output.extend([header, plot_combined]) result = super().correct(result) # adds check table; no auto-figure since list is empty - + + # FIXME: units for qick and OPX need to be reconciled here self.report_output.append( f"**Detuning (χ): {self.chi:.3f} MHz** " f"(f_0 before: {self.f0_before:.3f} MHz, f_0 after: {self.f0_after:.3f} MHz)\n\n" diff --git a/src/cqedtoolbox/protocols/operations/single_qubit/res_spec_vs_gain.py b/src/cqedtoolbox/protocols/operations/single_qubit/res_spec_vs_gain.py index 5dcf8e2..97f412b 100644 --- a/src/cqedtoolbox/protocols/operations/single_qubit/res_spec_vs_gain.py +++ b/src/cqedtoolbox/protocols/operations/single_qubit/res_spec_vs_gain.py @@ -7,7 +7,6 @@ plt.switch_backend("agg") - from labcore.analysis import DatasetAnalysis from labcore.measurement.storage import run_and_save_sweep from labcore.data.datadict_storage import datadict_from_hdf5, load_as_xr @@ -15,7 +14,8 @@ from labcore.measurement.record import recording, dep, indep from labcore.protocols.base import (ProtocolOperation, OperationStatus, serialize_fit_params, - CorrectionParameter, CheckResult, Correction, EvaluateResult) + CorrectionParameter, CheckResult, Correction, EvaluateResult, PlatformTypes) +from cqedtoolbox.protocols.operations import ResonatorGeometry from cqedtoolbox.protocols.parameters import ( Repetition, StartReadoutFrequency, @@ -24,7 +24,10 @@ ReadoutLength, StartReadoutGain, EndReadoutGain, ResonatorSpecSteps, ResonatorSpecVsGainSteps, nestedAttributeFromString, ) -from cqedtoolbox.protocols.operations.single_qubit.res_spec import ResonatorSpectroscopy, SyntheticHangerResonatorData +from cqedtoolbox.protocols.operations.single_qubit.res_spec import ( + ResonatorSpectroscopy, + SyntheticHangerResonatorData, +) from cqedtoolbox.measurement_lib.opx.advanced.qubit_tuneup import measure_pulse_resonator_spec_vs_readout_amp from cqedtoolbox.measurement_lib.qick.single_transmon_v2 import FreqGainSweepProgram @@ -156,10 +159,22 @@ class ResonatorSpectroscopyVsGain(ProtocolOperation): _SIM_N_GAIN_STEPS = 11 - def __init__(self, params): + def __init__(self, params, geometry: ResonatorGeometry | str): super().__init__() self.params = params + if isinstance(geometry, str): + try: + geometry = ResonatorGeometry(geometry.lower()) + except ValueError as err: + valid = ", ".join(g.value for g in ResonatorGeometry) + raise ValueError( + f"Unsupported resonator geometry '{geometry}'. Expected one of: {valid}" + ) from err + + self.geometry = geometry + self._fit_cls = geometry.fit_cls + self._register_inputs( repetitions=Repetition(params), readout_steps=ResonatorSpecSteps(params), @@ -203,6 +218,35 @@ def __init__(self, params): self.deviations = [] self.fit_results = [] self.snr_values = [] + self.trace_valid = [] + self.trace_invalid_reasons = [] + + def _trace_in_range(self, freqs, fit_result) -> bool: + f0 = fit_result.params["f_0"].value + freqs = np.asarray(freqs, dtype=float) + return float(np.min(freqs)) <= f0 <= float(np.max(freqs)) + + def _trace_invalid_reason(self, freqs, fit_result, snr) -> str | None: + if not self._trace_in_range(freqs, fit_result): + f0 = fit_result.params["f_0"].value + fmin = float(np.min(freqs)) + fmax = float(np.max(freqs)) + return f"f_0={f0:.6g} outside sweep=[{fmin:.6g}, {fmax:.6g}]" + + if snr < self.snr_threshold(): + return f"SNR={snr:.3f} < {self.snr_threshold():.3f}" + + max_error = self.max_fit_param_error() + for pname, param in fit_result.params.items(): + if pname in ["transmission_slope", "phase_slope", "phase_offset"]: + continue + if param.stderr is None: + return f"{pname}: no stderr" + if param.value == 0 or abs(param.stderr / param.value) > max_error: + pct = abs(param.stderr / param.value) * 100 if param.value != 0 else float("inf") + return f"{pname}: {pct:.0f}% error" + + return None def _measure_qick(self) -> Path: logger.info("Starting qick resonator spectroscopy vs gain measurement") @@ -341,18 +385,27 @@ def analyze(self): self.figure_paths.append(image_path) # Analyze each gain trace individually - gains = self.independents["gains"][0] + gains = self.independents["gains"] res_f_arr = [] + self.fit_results = [] + self.snr_values = [] + self.trace_valid = [] + self.trace_invalid_reasons = [] for i, g in enumerate(gains): - trace_signal = self.dependents["signal"].T[i] # Transpose to achieve gain as axis 0 - freqs = self.independents["frequencies"].T[i] - folder_name = f"resonator_spec_vs_gain_i={i}_g={g}" + # FIXME: is this the correct pattern for qick? + if self.platform_type == PlatformTypes.OPX: + trace_signal = self.dependents["signal"][i] + freqs = self.independents["frequencies"] + else: + trace_signal = self.dependents["signal"].T[i] # Transpose to achieve gain as axis 0 + freqs = self.independents["frequencies"].T[i] + # Use the static method from ResonatorSpectroscopy ret = ResonatorSpectroscopy.add_mag_and_unwind_and_fit( - freqs, trace_signal, f"Gain = {g}" + freqs, trace_signal, self.platform_type, self._fit_cls, f"Gain = {g}" ) _excluded = {"transmission_slope", "phase_slope", "phase_offset"} @@ -366,12 +419,15 @@ def analyze(self): f"{_null_stderr_params} — re-fitting" ) ret = ResonatorSpectroscopy.add_mag_and_unwind_and_fit( - freqs, trace_signal, f"Gain = {g}" + freqs, trace_signal, self.platform_type, self._fit_cls, f"Gain = {g}" ) self.fit_results.append(ret.fit_result) self.snr_values.append(ret.snr) res_f_arr.append(ret.fit_result.params["f_0"].value) + invalid_reason = self._trace_invalid_reason(freqs, ret.fit_result, ret.snr) + self.trace_valid.append(invalid_reason is None) + self.trace_invalid_reasons.append(invalid_reason) # Save individual trace analysis with DatasetAnalysis(self.data_loc, folder_name) as trace_ds: @@ -388,21 +444,7 @@ def analyze(self): self.resonance_frequencies = res_f_arr - # Find optimal gain: highest-SNR trace that passes the full quality check - snr_threshold = self.snr_threshold() - max_error = self.max_fit_param_error() - passing_indices = [] - for i, (snr, fit) in enumerate(zip(self.snr_values, self.fit_results)): - if snr < snr_threshold: - continue - bad_fit = any( - (param.stderr is None or param.value == 0 or - abs(param.stderr / param.value) > max_error) - for pname, param in fit.params.items() - if pname not in ["transmission_slope", "phase_slope", "phase_offset"] - ) - if not bad_fit: - passing_indices.append(i) + passing_indices = [i for i, valid in enumerate(self.trace_valid) if valid] if passing_indices: best_idx = max(passing_indices, key=lambda i: self.snr_values[i]) @@ -448,19 +490,8 @@ def _check_low_gain_quality(self) -> CheckResult: failures = [] for i in range(n_low): - snr = self.snr_values[i] - fit = self.fit_results[i] - if snr < threshold: - failures.append(f"trace {i}: SNR={snr:.3f} < {threshold:.3f}") - continue - for pname, param in fit.params.items(): - if pname in ["transmission_slope", "phase_slope", "phase_offset"]: - continue - if param.stderr is None: - failures.append(f"trace {i}/{pname}: no stderr") - elif param.value == 0 or abs(param.stderr / param.value) > max_error: - pct = abs(param.stderr / param.value) * 100 if param.value != 0 else float("inf") - failures.append(f"trace {i}/{pname}: {pct:.0f}%") + if not self.trace_valid[i]: + failures.append(f"trace {i}: {self.trace_invalid_reasons[i]}") passed = len(failures) == 0 desc = (f"first {n_low} traces pass quality check" if passed @@ -468,12 +499,21 @@ def _check_low_gain_quality(self) -> CheckResult: return CheckResult("low_gain_quality_check", passed, desc) def _check_high_snr(self) -> CheckResult: - """At least one trace must exceed the high SNR threshold.""" + """At least one valid trace must exceed the high SNR threshold.""" threshold = self.high_snr_threshold() - best = max(self.snr_values) + valid_indices = [i for i, valid in enumerate(self.trace_valid) if valid] + if not valid_indices: + return CheckResult( + "high_snr_check", + False, + "no valid traces remain after fit-range and fit-quality filtering", + ) + + best_idx = max(valid_indices, key=lambda i: self.snr_values[i]) + best = self.snr_values[best_idx] passed = best >= threshold - desc = (f"best SNR={best:.3f} ≥ {threshold:.3f}" if passed - else f"best SNR={best:.3f} < {threshold:.3f} — no trace meets high-SNR threshold") + desc = (f"best valid trace={best_idx}, SNR={best:.3f} ≥ {threshold:.3f}" if passed + else f"best valid trace={best_idx}, SNR={best:.3f} < {threshold:.3f}") return CheckResult("high_snr_check", passed, desc) def correct(self, result: EvaluateResult) -> EvaluateResult: @@ -499,13 +539,15 @@ def correct(self, result: EvaluateResult) -> EvaluateResult: result = super().correct(result) # check table + success update (writes readout_gain) if result.status == OperationStatus.SUCCESS: - gains = self.independents["gains"][0] + gains = self.independents["gains"] self.report_output.append("\n### Individual Gain Traces\n") for i, (fig_path, g) in enumerate(zip(trace_figures, gains)): + validity = "valid" if self.trace_valid[i] else f"invalid ({self.trace_invalid_reasons[i]})" self.report_output.extend([ f"\n**Trace {i}: Gain = {g:.3f}**\n" f"- SNR: {self.snr_values[i]:.3f}\n" - f"- f_0: {self.resonance_frequencies[i]:.3f} MHz\n", + f"- f_0: {self.resonance_frequencies[i]:.3f} MHz\n" + f"- Status: {validity}\n", fig_path, ]) diff --git a/src/cqedtoolbox/protocols/operations/single_qubit/sat_spec.py b/src/cqedtoolbox/protocols/operations/single_qubit/sat_spec.py index b74cdc1..aa80b0d 100644 --- a/src/cqedtoolbox/protocols/operations/single_qubit/sat_spec.py +++ b/src/cqedtoolbox/protocols/operations/single_qubit/sat_spec.py @@ -27,6 +27,7 @@ ) from cqedtoolbox.measurement_lib.opx.advanced.qubit_tuneup import measure_qubit_ssb_spec_saturation from cqedtoolbox.measurement_lib.qick.single_transmon_v2 import PulseProbeSpectroscopy +from cqedtoolbox.readout.qubit_readout import rotate_complex_qubit_data logger = logging.getLogger(__name__) @@ -552,17 +553,9 @@ def __init__(self, params): self.independents = {"frequencies": []} self.dependents = {"signal": []} - self.fit_result_re = None - self.fit_result_imag = None - self.fit_result_mag = None - self.snr_re = None - self.snr_imag = None - self.snr_mag = None - self.residuals_re = None - self.residuals_imag = None - self.residuals_mag = None - self._best_fit_result = None - self._winner_key: str | None = None + self.fit_result = None + self.snr = None + self.residuals = None def _measure_qick(self) -> Path: logger.info("Starting qick saturation spectroscopy measurement") @@ -581,18 +574,18 @@ def _measure_opx(self) -> Path: return loc def _load_data_qick(self): - path = self.data_loc / "data.ddh5" - if not path.exists(): - raise FileNotFoundError(f"File {path} does not exist") - data = datadict_from_hdf5(path) - - self.independents["frequencies"] = data["freq"]["values"] - self.dependents["signal"] = data["signal"]["values"] + data = load_as_xr(self.data_loc) + rotated = rotate_complex_qubit_data(data)[0] + self.independents["frequencies"] = rotated["freq"].values + self.dependents["signal"] = rotated["signal"].values def _load_data_opx(self): - data = load_as_xr(self.data_loc).mean("repetition") + data = load_as_xr(self.data_loc) + if "repetition" in data.dims: + data = data.mean("repetition") + data, _ = rotate_complex_qubit_data(data) self.independents["frequencies"] = data["ssb_frequency"].values - self.dependents["signal"] = data["signal_Re"].values + 1j * data["signal_Im"].values + self.dependents["signal"] = data["signal"].values def _measure_dummy(self) -> Path: """Create synthetic saturation spectroscopy data using a sweep""" @@ -624,148 +617,58 @@ def _measure_dummy(self) -> Path: def _load_data_dummy(self): """Load dummy data from disk (same as _load_data_qick)""" logger.info("Loading dummy data from disk") - path = self.data_loc / "data.ddh5" - if not path.exists(): - raise FileNotFoundError(f"File {path} does not exist") - data = datadict_from_hdf5(path) - - self.independents["frequencies"] = data["frequencies"]["values"] - self.dependents["signal"] = data["signal"]["values"] - - def _fit_lorentzian_components(self, frequencies, signal, fig_title="") -> tuple: - """ - Fit real, imaginary, and magnitude components with Lorentzian fits. - Returns (fit_result_re, residuals_re, snr_re), (fit_result_imag, ...), (fit_result_mag, ...), - fig_re, fig_imag, fig_mag - """ - signal_re = signal.real - signal_imag = signal.imag - signal_mag = np.abs(signal) - - # Fit real part - fit_re = Lorentzian(frequencies, signal_re) - fit_result_re = fit_re.run(fit_re) - fit_curve_re = fit_result_re.eval() - residuals_re = signal_re - fit_curve_re - amp_re = fit_result_re.params["A"].value - noise_re = np.std(residuals_re) - snr_re = np.abs(amp_re / (4 * noise_re)) - - # Fit imaginary part - fit_imag = Lorentzian(frequencies, signal_imag) - fit_result_imag = fit_imag.run(fit_imag) - fit_curve_imag = fit_result_imag.eval() - residuals_imag = signal_imag - fit_curve_imag - amp_imag = fit_result_imag.params["A"].value - noise_imag = np.std(residuals_imag) - snr_imag = np.abs(amp_imag / (4 * noise_imag)) - - # Fit magnitude - fit_mag = Lorentzian(frequencies, signal_mag) - fit_result_mag = fit_mag.run(fit_mag) - fit_curve_mag = fit_result_mag.eval() - residuals_mag = signal_mag - fit_curve_mag - amp_mag = fit_result_mag.params["A"].value - noise_mag = np.std(residuals_mag) - snr_mag = np.abs(amp_mag / (4 * noise_mag)) - - # Real plot - fig_re, ax_re = plt.subplots() - ax_re.set_title(f"{fig_title} - Real") - ax_re.set_xlabel("Frequency (Hz)") - ax_re.set_ylabel("Signal Real (A.U)") - ax_re.plot(frequencies, signal_re, label="Data") - ax_re.plot(frequencies, fit_curve_re, label="Fit") - ax_re.legend() - - # Imaginary plot - fig_imag, ax_imag = plt.subplots() - ax_imag.set_title(f"{fig_title} - Imaginary") - ax_imag.set_xlabel("Frequency (Hz)") - ax_imag.set_ylabel("Signal Imaginary (A.U)") - ax_imag.plot(frequencies, signal_imag, label="Data") - ax_imag.plot(frequencies, fit_curve_imag, label="Fit") - ax_imag.legend() - - # Magnitude plot - fig_mag, ax_mag = plt.subplots() - ax_mag.set_title(f"{fig_title} - Magnitude") - ax_mag.set_xlabel("Frequency (Hz)") - ax_mag.set_ylabel("Signal Magnitude (A.U)") - ax_mag.plot(frequencies, signal_mag, label="Data") - ax_mag.plot(frequencies, fit_curve_mag, label="Fit") - ax_mag.legend() - - return ( - (fit_result_re, residuals_re, snr_re), - (fit_result_imag, residuals_imag, snr_imag), - (fit_result_mag, residuals_mag, snr_mag), - fig_re, - fig_imag, - fig_mag - ) + data = load_as_xr(self.data_loc) + rotated = rotate_complex_qubit_data(data)[0] + self.independents["frequencies"] = rotated["frequencies"].values + self.dependents["signal"] = rotated["signal"].values + + def _fit_lorentzian(self, frequencies, signal, fig_title="") -> tuple: + fit = Lorentzian(frequencies, signal) + fit_result = fit.run(fit) + fit_curve = fit_result.eval() + residuals = signal - fit_curve + amp = fit_result.params["A"].value + noise = np.std(residuals) + snr = np.abs(amp / (4 * noise)) + + fig, ax = plt.subplots() + ax.set_title(fig_title) + ax.set_xlabel("Frequency (Hz)") + ax.set_ylabel("Rotated Signal (A.U)") + ax.plot(frequencies, signal, label="Data") + ax.plot(frequencies, fit_curve, label="Fit") + ax.legend() + + return fit_result, residuals, snr, fig def analyze(self): with DatasetAnalysis(self.data_loc, self.name) as ds: - result_re, result_imag, result_mag, fig_re, fig_imag, fig_mag = self._fit_lorentzian_components( + self.fit_result, self.residuals, self.snr, fig = self._fit_lorentzian( self.independents["frequencies"], self.dependents["signal"], "Saturation Spectroscopy" ) - - self.fit_result_re, self.residuals_re, self.snr_re = result_re - self.fit_result_imag, self.residuals_imag, self.snr_imag = result_imag - self.fit_result_mag, self.residuals_mag, self.snr_mag = result_mag - - # Determine winner (highest SNR) for use in checks and success update - snr_map = {"re": self.snr_re, "imag": self.snr_imag, "mag": self.snr_mag} - self._winner_key = max(snr_map, key=lambda k: snr_map[k]) - fit_map = { - "re": self.fit_result_re, - "imag": self.fit_result_imag, - "mag": self.fit_result_mag, - } - self._best_fit_result = fit_map[self._winner_key] + self._best_fit_result = self.fit_result ds.add( - fit_result_re=self.fit_result_re, - params_re=serialize_fit_params(self.fit_result_re.params), - snr_re=float(self.snr_re), - fit_result_imag=self.fit_result_imag, - params_imag=serialize_fit_params(self.fit_result_imag.params), - snr_imag=float(self.snr_imag), - fit_result_mag=self.fit_result_mag, - params_mag=serialize_fit_params(self.fit_result_mag.params), - snr_mag=float(self.snr_mag) + fit_result=self.fit_result, + params=serialize_fit_params(self.fit_result.params), + snr=float(self.snr) ) - ds.add_figure(f"{self.name}_real", fig=fig_re) - image_path_re = ds._new_file_path(ds.savefolders[1], f"{self.name}_real", suffix="png") - self.figure_paths.append(image_path_re) - - ds.add_figure(f"{self.name}_imag", fig=fig_imag) - image_path_imag = ds._new_file_path(ds.savefolders[1], f"{self.name}_imag", suffix="png") - self.figure_paths.append(image_path_imag) - - ds.add_figure(f"{self.name}_mag", fig=fig_mag) - image_path_mag = ds._new_file_path(ds.savefolders[1], f"{self.name}_mag", suffix="png") - self.figure_paths.append(image_path_mag) + ds.add_figure(self.name, fig=fig) + image_path = ds._new_file_path(ds.savefolders[1], self.name, suffix="png") + self.figure_paths.append(image_path) # --- checks (pure assessment) --- def _check_fit_quality(self) -> CheckResult: - snr_map = {"re": self.snr_re, "imag": self.snr_imag, "mag": self.snr_mag} - winner_key = max(snr_map, key=lambda k: snr_map[k]) - winner_snr = snr_map[winner_key] - winner_fit = {"re": self.fit_result_re, "imag": self.fit_result_imag, "mag": self.fit_result_mag}[winner_key] - label_map = {"re": "Real", "imag": "Imaginary", "mag": "Magnitude"} - threshold = self.snr_threshold() - snr_passed = winner_snr >= threshold + snr_passed = self.snr >= threshold max_error = self.max_fit_param_error() bad_params = [] - for pname, param in winner_fit.params.items(): + for pname, param in self.fit_result.params.items(): if pname == "of": continue if param.stderr is None: @@ -775,19 +678,15 @@ def _check_fit_quality(self) -> CheckResult: bad_params.append(f"{pname}({pct:.0f}%)") passed = snr_passed and len(bad_params) == 0 - parts = [f"best_SNR={winner_snr:.3f} (threshold={threshold:.3f}, component={label_map[winner_key]})"] + parts = [f"SNR={self.snr:.3f} (threshold={threshold:.3f})"] if bad_params: parts.append(f"high-error params: {', '.join(bad_params)}") return CheckResult("fit_quality", passed, "; ".join(parts)) def _check_single_peak(self) -> CheckResult: - residuals_map = { - "re": self.residuals_re, - "imag": self.residuals_imag, - "mag": self.residuals_mag, - } - residuals = residuals_map[self._winner_key] + # TODO: make sure that the fit frequency is inside the swept range + residuals = self.residuals frequencies = self.independents["frequencies"] fit = Lorentzian(frequencies, residuals) @@ -827,45 +726,21 @@ def _check_single_peak(self) -> CheckResult: return CheckResult("single_peak", passed, description) def correct(self, result: EvaluateResult) -> EvaluateResult: - # Pop all figures before super() to control layout - fig_re = self.figure_paths[0] if len(self.figure_paths) >= 1 else None - fig_imag = self.figure_paths[1] if len(self.figure_paths) >= 2 else None - fig_mag = self.figure_paths[2] if len(self.figure_paths) >= 3 else None + figure = self.figure_paths[0] if self.figure_paths else None self.figure_paths.clear() - # Build component report (winner first) - snr_map = {"re": self.snr_re, "imag": self.snr_imag, "mag": self.snr_mag} - fig_map = {"re": fig_re, "imag": fig_imag, "mag": fig_mag} - fit_map = { - "re": self.fit_result_re, - "imag": self.fit_result_imag, - "mag": self.fit_result_mag, - } - label_map = {"re": "Real", "imag": "Imaginary", "mag": "Magnitude"} - - sorted_keys = sorted(snr_map, key=lambda k: snr_map[k], reverse=True) - winner_key = sorted_keys[0] - header = (f"## Saturation Spectroscopy\n" f"Measured saturation spectroscopy for frequencies: " f"{self.start_freq():.3f}–{self.end_freq():.3f} MHz\n" f"Data Path: `{self.data_loc}`\n\n") self.report_output.append(header) - winner_label = label_map[winner_key] self.report_output.extend([ - f"### **{winner_label} Component (SELECTED, SNR={snr_map[winner_key]:.3f})**\n\n", - fig_map[winner_key], - f"**Fit Report:**\n```\n{fit_map[winner_key].lmfit_result.fit_report()}\n```\n\n", + f"### Rotated Signal Fit (SNR={self.snr:.3f})\n\n", + figure, + f"**Fit Report:**\n```\n{self.fit_result.lmfit_result.fit_report()}\n```\n\n", ]) - for k in sorted_keys[1:]: - self.report_output.extend([ - f"### **{label_map[k]} Component (SNR={snr_map[k]:.3f})**\n\n", - fig_map[k], - f"**Fit Report:**\n```\n{fit_map[k].lmfit_result.fit_report()}\n```\n\n", - ]) - # Let super() add the check table; no auto-figure since figure_paths is cleared result = super().correct(result) diff --git a/src/cqedtoolbox/protocols/operations/single_qubit/t1.py b/src/cqedtoolbox/protocols/operations/single_qubit/t1.py index 253c2c5..8ebd6d4 100644 --- a/src/cqedtoolbox/protocols/operations/single_qubit/t1.py +++ b/src/cqedtoolbox/protocols/operations/single_qubit/t1.py @@ -16,7 +16,7 @@ from labcore.protocols.base import ( ProtocolOperation, serialize_fit_params, - CorrectionParameter, CheckResult, Correction, EvaluateResult, + CorrectionParameter, CheckResult, Correction, EvaluateResult, PlatformTypes ) from cqedtoolbox.protocols.parameters import ( Repetition, @@ -29,6 +29,7 @@ ) from cqedtoolbox.measurement_lib.opx.advanced.qubit_tuneup import measure_t1 from cqedtoolbox.measurement_lib.qick.single_transmon_v2 import T1Program +from cqedtoolbox.readout.qubit_readout import rotate_complex_qubit_data logger = logging.getLogger(__name__) @@ -227,23 +228,15 @@ def __init__(self, params): self._register_check("quality_check", self._check_quality, [self._increase_delay, self._increase_averaging]) - self._register_success_update(self.t1, lambda: self._winner_fit.params["tau"].value) - self._register_success_update(self.delay, lambda: 10 * self._winner_fit.params["tau"].value) + self._register_success_update(self.t1, lambda: self.fit_result.params["tau"].value) + self._register_success_update(self.delay, lambda: 10 * self.fit_result.params["tau"].value) self.independents = {"delays": []} self.dependents = {"signal": []} - self.fit_result_re = None - self.fit_result_imag = None - self.fit_result_mag = None - self.snr_re = None - self.snr_imag = None - self.snr_mag = None - self._winner_fit = None - self._winner_snr = None - self._winner_key = None - self._winner_name = None - self._sorted_components = None + self.fit_result = None + self.residuals = None + self.snr = None def _measure_dummy(self) -> Path: logger.info("Starting dummy T1 measurement") @@ -256,12 +249,10 @@ def _measure_dummy(self) -> Path: return loc def _load_data_dummy(self): - path = self.data_loc / "data.ddh5" - if not path.exists(): - raise FileNotFoundError(f"File {path} does not exist") - data = datadict_from_hdf5(path) - self.independents["delays"] = data["delays"]["values"] - self.dependents["signal"] = data["signal"]["values"] + data = load_as_xr(self.data_loc) + rotated = rotate_complex_qubit_data(data)[0] + self.independents["delays"] = rotated["delays"].values + self.dependents["signal"] = rotated["signal"].values def _measure_qick(self) -> Path: logger.info("Starting qick T1 measurement") @@ -278,178 +269,80 @@ def _measure_opx(self) -> Path: return loc def _load_data_qick(self): - path = self.data_loc / "data.ddh5" - if not path.exists(): - raise FileNotFoundError(f"File {path} does not exist") - data = datadict_from_hdf5(path) - - self.independents["delays"] = data["t"]["values"] - self.dependents["signal"] = data["signal"]["values"] + data = load_as_xr(self.data_loc) + rotated = rotate_complex_qubit_data(data)[0] + self.independents["delays"] = rotated["t"].values + self.dependents["signal"] = rotated["signal"].values def _load_data_opx(self): - data = load_as_xr(self.data_loc).mean("repetition") + data = load_as_xr(self.data_loc) + if "repetition" in data.dims: + data = data.mean("repetition") + data, _ = rotate_complex_qubit_data(data) self.independents["delays"] = data["delay"].values - self.dependents["signal"] = data["signal_Re"].values + 1j * data["signal_Im"].values - - def _fit_exponential_components(self, delays, signal, fig_title="") -> tuple: - """ - Fit real, imaginary, and magnitude components with ExponentialDecay fits. - Returns (fit_result_re, fit_result_imag, fit_result_mag, fig_re, fig_imag, fig_mag) - """ - signal_re = signal.real - signal_imag = signal.imag - signal_mag = np.abs(signal) - - # Fit real part - fit_re = ExponentialDecay(delays, signal_re) - fit_result_re = fit_re.run(fit_re) - fit_curve_re = fit_result_re.eval() - residuals_re = signal_re - fit_curve_re - amp_re = fit_result_re.params["A"].value - noise_re = np.std(residuals_re) - snr_re = np.abs(amp_re / (4 * noise_re)) - - # Fit imaginary part - fit_imag = ExponentialDecay(delays, signal_imag) - fit_result_imag = fit_imag.run(fit_imag) - fit_curve_imag = fit_result_imag.eval() - residuals_imag = signal_imag - fit_curve_imag - amp_imag = fit_result_imag.params["A"].value - noise_imag = np.std(residuals_imag) - snr_imag = np.abs(amp_imag / (4 * noise_imag)) - - # Fit magnitude - fit_mag = ExponentialDecay(delays, signal_mag) - fit_result_mag = fit_mag.run(fit_mag) - fit_curve_mag = fit_result_mag.eval() - residuals_mag = signal_mag - fit_curve_mag - amp_mag = fit_result_mag.params["A"].value - noise_mag = np.std(residuals_mag) - snr_mag = np.abs(amp_mag / (4 * noise_mag)) - - # Create three separate figures - # Real plot - fig_re, ax_re = plt.subplots() - ax_re.set_title(f"{fig_title} - Real") - ax_re.set_xlabel("Delay (μs)") - ax_re.set_ylabel("Signal Real (A.U)") - ax_re.plot(delays, signal_re, label="Data") - ax_re.plot(delays, fit_curve_re, label="Fit") - ax_re.legend() - - # Imaginary plot - fig_imag, ax_imag = plt.subplots() - ax_imag.set_title(f"{fig_title} - Imaginary") - ax_imag.set_xlabel("Delay (μs)") - ax_imag.set_ylabel("Signal Imaginary (A.U)") - ax_imag.plot(delays, signal_imag, label="Data") - ax_imag.plot(delays, fit_curve_imag, label="Fit") - ax_imag.legend() - - # Magnitude plot - fig_mag, ax_mag = plt.subplots() - ax_mag.set_title(f"{fig_title} - Magnitude") - ax_mag.set_xlabel("Delay (μs)") - ax_mag.set_ylabel("Signal Magnitude (A.U)") - ax_mag.plot(delays, signal_mag, label="Data") - ax_mag.plot(delays, fit_curve_mag, label="Fit") - ax_mag.legend() - - return ( - (fit_result_re, residuals_re, snr_re), - (fit_result_imag, residuals_imag, snr_imag), - (fit_result_mag, residuals_mag, snr_mag), - fig_re, - fig_imag, - fig_mag - ) + self.dependents["signal"] = data["signal"].values + + def _fit_exponential(self, delays, signal, fig_title="") -> tuple: + fit = ExponentialDecay(delays, signal) + fit_result = fit.run(fit) + fit_curve = fit_result.eval() + residuals = signal - fit_curve + amp = fit_result.params["A"].value + noise = np.std(residuals) + snr = np.abs(amp / (4 * noise)) + + fig, ax = plt.subplots() + ax.set_title(fig_title) + if self.platform_type == PlatformTypes.OPX: + ax.set_xlabel("Delay (ns)") + else: + ax.set_xlabel("Delay (μs)") + ax.set_ylabel("Rotated Signal (A.U)") + ax.plot(delays, signal, label="Data") + ax.plot(delays, fit_curve, label="Fit") + ax.legend() + + return fit_result, residuals, snr, fig def analyze(self): with DatasetAnalysis(self.data_loc, self.name) as ds: - result_re, result_imag, result_mag, fig_re, fig_imag, fig_mag = self._fit_exponential_components( + self.fit_result, self.residuals, self.snr, fig = self._fit_exponential( self.independents["delays"], self.dependents["signal"], "T1 Measurement" ) - self.fit_result_re, residuals_re, self.snr_re = result_re - self.fit_result_imag, residuals_imag, self.snr_imag = result_imag - self.fit_result_mag, residuals_mag, self.snr_mag = result_mag - # Save all fit results ds.add( - fit_result_re=self.fit_result_re, - params_re=serialize_fit_params(self.fit_result_re.params), - snr_re=float(self.snr_re), - fit_result_imag=self.fit_result_imag, - params_imag=serialize_fit_params(self.fit_result_imag.params), - snr_imag=float(self.snr_imag), - fit_result_mag=self.fit_result_mag, - params_mag=serialize_fit_params(self.fit_result_mag.params), - snr_mag=float(self.snr_mag) + fit_result=self.fit_result, + params=serialize_fit_params(self.fit_result.params), + snr=float(self.snr) ) - # Save all three figures separately - ds.add_figure(f"{self.name}_real", fig=fig_re) - image_path_re = ds._new_file_path(ds.savefolders[1], f"{self.name}_real", suffix="png") - self.figure_paths.append(image_path_re) - - ds.add_figure(f"{self.name}_imag", fig=fig_imag) - image_path_imag = ds._new_file_path(ds.savefolders[1], f"{self.name}_imag", suffix="png") - self.figure_paths.append(image_path_imag) - - ds.add_figure(f"{self.name}_mag", fig=fig_mag) - image_path_mag = ds._new_file_path(ds.savefolders[1], f"{self.name}_mag", suffix="png") - self.figure_paths.append(image_path_mag) - - snr_dict = { - "Real": (self.snr_re, self.fit_result_re, "re"), - "Imaginary": (self.snr_imag, self.fit_result_imag, "imag"), - "Magnitude": (self.snr_mag, self.fit_result_mag, "mag"), - } - self._sorted_components = sorted(snr_dict.items(), key=lambda x: x[1][0], reverse=True) + ds.add_figure(self.name, fig=fig) + image_path = ds._new_file_path(ds.savefolders[1], self.name, suffix="png") + self.figure_paths.append(image_path) def _check_quality(self) -> CheckResult: snr_min = self.snr_min_threshold() - - valid = [ - (name, snr, fit, key) - for name, (snr, fit, key) in self._sorted_components - if snr >= snr_min - ] - - if valid: - self._winner_name, self._winner_snr, self._winner_fit, self._winner_key = valid[0] - max_error = self.max_fit_param_error() - bad_params = [] - for pname, param in self._winner_fit.params.items(): - if param.stderr is None: - bad_params.append(f"{pname}(no stderr)") - elif param.value == 0 or abs(param.stderr / param.value) > max_error: - pct = abs(param.stderr / param.value) * 100 if param.value != 0 else float("inf") - bad_params.append(f"{pname}({pct:.0f}%)") - passed = len(bad_params) == 0 - parts = [f"winner={self._winner_name}, SNR={self._winner_snr:.3f} (threshold={snr_min:.1f})"] - if bad_params: - parts.append(f"high-error params: {', '.join(bad_params)}") - else: - self._winner_name, (self._winner_snr, self._winner_fit, self._winner_key) = self._sorted_components[0] - passed = False - parts = [ - f"no component with SNR >= {snr_min:.1f}", - f"best={self._winner_name}, SNR={self._winner_snr:.3f}", - ] + max_error = self.max_fit_param_error() + bad_params = [] + for pname, param in self.fit_result.params.items(): + if param.stderr is None: + bad_params.append(f"{pname}(no stderr)") + elif param.value == 0 or abs(param.stderr / param.value) > max_error: + pct = abs(param.stderr / param.value) * 100 if param.value != 0 else float("inf") + bad_params.append(f"{pname}({pct:.0f}%)") + passed = self.snr >= snr_min and len(bad_params) == 0 + parts = [f"SNR={self.snr:.3f} (threshold={snr_min:.1f})"] + if bad_params: + parts.append(f"high-error params: {', '.join(bad_params)}") return CheckResult("quality_check", passed, "; ".join(parts)) def correct(self, result: EvaluateResult) -> EvaluateResult: - # Pop all figures before super() auto-appends the last one - fig_re = self.figure_paths.pop(0) if len(self.figure_paths) >= 3 else None - fig_imag = self.figure_paths.pop(0) if self.figure_paths else None - fig_mag = self.figure_paths.pop(0) if self.figure_paths else None - self.figure_paths.clear() # prevent auto-append - - plot_map = {"re": fig_re, "imag": fig_imag, "mag": fig_mag} + figure = self.figure_paths[0] if self.figure_paths else None + self.figure_paths.clear() snr_min = self.snr_min_threshold() self.report_output.append( @@ -458,15 +351,13 @@ def correct(self, result: EvaluateResult) -> EvaluateResult: f"Data Path: `{self.data_loc}`\n\n" ) - for i, (comp_name, (comp_snr, comp_fit, comp_key)) in enumerate(self._sorted_components): - tag = "(SELECTED)" if i == 0 else "(NOT SELECTED)" - self.report_output.append(f"### **{comp_name} Component {tag}**\n") - if plot_map.get(comp_key): - self.report_output.append(plot_map[comp_key]) - self.report_output.append( - f"SNR={comp_snr:.3f}\n\n" - f"**Fit Report:**\n```\n{str(comp_fit.lmfit_result.fit_report())}\n```\n\n" - ) + self.report_output.append("### Rotated Signal Fit\n") + if figure: + self.report_output.append(figure) + self.report_output.append( + f"SNR={self.snr:.3f}\n\n" + f"**Fit Report:**\n```\n{str(self.fit_result.lmfit_result.fit_report())}\n```\n\n" + ) result = super().correct(result) # adds check table + success update lines return result diff --git a/src/cqedtoolbox/protocols/operations/single_qubit/t2e.py b/src/cqedtoolbox/protocols/operations/single_qubit/t2e.py index 7b8d79d..3de02b0 100644 --- a/src/cqedtoolbox/protocols/operations/single_qubit/t2e.py +++ b/src/cqedtoolbox/protocols/operations/single_qubit/t2e.py @@ -16,7 +16,7 @@ from labcore.protocols.base import ( ProtocolOperation, PlatformTypes, serialize_fit_params, - CorrectionParameter, CheckResult, Correction, EvaluateResult, + CorrectionParameter, CheckResult, Correction, EvaluateResult, PlatformTypes ) from cqedtoolbox.protocols.parameters import ( Repetition, @@ -25,10 +25,10 @@ ReadoutGain, ReadoutLength, T2E, - NEchos ) from cqedtoolbox.measurement_lib.opx.advanced.qubit_tuneup import measure_t2 from cqedtoolbox.measurement_lib.qick.single_transmon_v2 import T2nProgram +from cqedtoolbox.readout.qubit_readout import rotate_complex_qubit_data logger = logging.getLogger(__name__) @@ -60,17 +60,6 @@ def _opx_getter(self): return self.params.corrections.t2e.max_fit_param_error() def _opx_setter(self, v): self.params.corrections.t2e.max_fit_param_error(v) -@dataclass -class MaxEchos(CorrectionParameter): - name: str = field(default="t2e_max_echos", init=False) - description: str = field(default="Maximum number of echo pulses to try", init=False) - - def _qick_getter(self): return int(self.params.corrections.t2e.max_echos()) - def _qick_setter(self, v): self.params.corrections.t2e.max_echos(v) - def _opx_getter(self): return int(self.params.corrections.t2e.max_echos()) - def _opx_setter(self, v): self.params.corrections.t2e.max_echos(v) - - @dataclass class AveragingIncreaseFactor(CorrectionParameter): name: str = field(default="t2e_averaging_factor", init=False) @@ -93,51 +82,13 @@ def _opx_getter(self): return int(self.params.corrections.t2e.max_averaging_incr def _opx_setter(self, v): self.params.corrections.t2e.max_averaging_increases(v) -# --------------------------------------------------------------------------- -# Correction subclasses -# --------------------------------------------------------------------------- - -class IncreaseEchosCorrection(Correction): - name = "increase_echos" - description = "Increase number of echo pulses by 1" - triggered_by = "quality_check" - - def __init__(self, n_echos_param, max_echos_param): - self.n_echos_param = n_echos_param - self.max_echos_param = max_echos_param - self._original_echos: int | None = None - self._last_change: str = "" - - def can_apply(self) -> bool: - if self._original_echos is None: - self._original_echos = int(self.n_echos_param()) - return int(self.n_echos_param()) < int(self.max_echos_param()) - - def apply(self) -> None: - if self._original_echos is None: - self._original_echos = int(self.n_echos_param()) - old = int(self.n_echos_param()) - new = old + 1 - self.n_echos_param(new) - self._last_change = f"n_echos: {old} → {new}" - - def report_output(self) -> str: - return self._last_change - - def reset(self) -> None: - if self._original_echos is not None: - self.n_echos_param(self._original_echos) - - class IncreaseAveragingCorrection(Correction): name = "increase_averaging" - description = "Increase number of repetitions and reset echo count" + description = "Increase number of repetitions" triggered_by = "quality_check" - def __init__(self, reps_param, echo_correction: IncreaseEchosCorrection, - factor_param, max_increases_param): + def __init__(self, reps_param, factor_param, max_increases_param): self.reps_param = reps_param - self.echo_correction = echo_correction self.factor_param = factor_param self.max_increases_param = max_increases_param self._original_reps: int | None = None @@ -156,7 +107,6 @@ def apply(self) -> None: self.reps_param(new) self._count += 1 self._last_change = f"reps: {old} → {new}" - self.echo_correction.reset() def report_output(self) -> str: return self._last_change @@ -184,7 +134,6 @@ def __init__(self, params): qubit_gain=QubitGain(params), readout_gain=ReadoutGain(params), readout_length=ReadoutLength(params), - n_echos=NEchos(params), ) self._register_outputs( t2e=T2E(params) @@ -193,14 +142,12 @@ def __init__(self, params): self._register_correction_params( snr_min_threshold=SNRMinThreshold(params), max_fit_param_error=MaxFitParamError(params), - max_echos=MaxEchos(params), averaging_increase_factor=AveragingIncreaseFactor(params), max_averaging_increases=MaxAveragingIncreases(params), ) self._increase_averaging = IncreaseAveragingCorrection( self.repetitions, - self._increase_echos, self.averaging_increase_factor, self.max_averaging_increases, ) @@ -208,22 +155,14 @@ def __init__(self, params): corrections = [self._increase_averaging] self._register_check("quality_check", self._check_quality, corrections) - self._register_success_update(self.t2e, lambda: self._winner_fit.params["tau"].value) + self._register_success_update(self.t2e, lambda: self.fit_result.params["tau"].value) self.independents = {"delays": []} self.dependents = {"signal": []} - self.fit_result_re = None - self.fit_result_imag = None - self.fit_result_mag = None - self.snr_re = None - self.snr_imag = None - self.snr_mag = None - self._winner_fit = None - self._winner_snr = None - self._winner_key = None - self._winner_name = None - self._sorted_components = None + self.fit_result = None + self.residuals = None + self.snr = None def _measure_dummy(self) -> Path: logger.info("Starting dummy T2 Echo measurement") @@ -236,12 +175,10 @@ def _measure_dummy(self) -> Path: return loc def _load_data_dummy(self): - path = self.data_loc / "data.ddh5" - if not path.exists(): - raise FileNotFoundError(f"File {path} does not exist") - data = datadict_from_hdf5(path) - self.independents["delays"] = data["delays"]["values"] - self.dependents["signal"] = data["signal"]["values"] + data = load_as_xr(self.data_loc) + rotated = rotate_complex_qubit_data(data)[0] + self.independents["delays"] = rotated["delays"].values + self.dependents["signal"] = rotated["signal"].values def _measure_qick(self) -> Path: logger.info("Starting qick T2 Echo measurement") @@ -258,175 +195,80 @@ def _measure_opx(self) -> Path: return loc def _load_data_qick(self): - path = self.data_loc / "data.ddh5" - if not path.exists(): - raise FileNotFoundError(f"File {path} does not exist") - data = datadict_from_hdf5(path) - - self.independents["delays"] = data["t"]["values"] - self.dependents["signal"] = data["signal"]["values"] + data = load_as_xr(self.data_loc) + rotated = rotate_complex_qubit_data(data)[0] + self.independents["delays"] = rotated["t"].values + self.dependents["signal"] = rotated["signal"].values def _load_data_opx(self): - data = load_as_xr(self.data_loc).mean("repetition") + data = load_as_xr(self.data_loc) + if "repetition" in data.dims: + data = data.mean("repetition") + data, _ = rotate_complex_qubit_data(data) self.independents["delays"] = data["delay"].values - self.dependents["signal"] = data["signal_Re"].values + 1j * data["signal_Im"].values - - def _fit_exponentially_decaying_sine_components(self, delays, signal, fig_title="") -> tuple: - """ - Fit real, imaginary, and magnitude components with ExponentiallyDecayingSine fits. - Returns (fit_result_re, fit_result_imag, fit_result_mag, fig_re, fig_imag, fig_mag) - """ - signal_re = signal.real - signal_imag = signal.imag - signal_mag = np.abs(signal) - - # Fit real part - fit_re = ExponentiallyDecayingSine(delays, signal_re) - fit_result_re = fit_re.run(fit_re) - fit_curve_re = fit_result_re.eval() - residuals_re = signal_re - fit_curve_re - amp_re = fit_result_re.params["A"].value - noise_re = np.std(residuals_re) - snr_re = np.abs(amp_re / (4 * noise_re)) - - # Fit imaginary part - fit_imag = ExponentiallyDecayingSine(delays, signal_imag) - fit_result_imag = fit_imag.run(fit_imag) - fit_curve_imag = fit_result_imag.eval() - residuals_imag = signal_imag - fit_curve_imag - amp_imag = fit_result_imag.params["A"].value - noise_imag = np.std(residuals_imag) - snr_imag = np.abs(amp_imag / (4 * noise_imag)) - - # Fit magnitude - fit_mag = ExponentiallyDecayingSine(delays, signal_mag) - fit_result_mag = fit_mag.run(fit_mag) - fit_curve_mag = fit_result_mag.eval() - residuals_mag = signal_mag - fit_curve_mag - amp_mag = fit_result_mag.params["A"].value - noise_mag = np.std(residuals_mag) - snr_mag = np.abs(amp_mag / (4 * noise_mag)) - - # Create three separate figures - fig_re, ax_re = plt.subplots() - ax_re.set_title(f"{fig_title} - Real") - ax_re.set_xlabel("Delay (μs)") - ax_re.set_ylabel("Signal Real (A.U)") - ax_re.plot(delays, signal_re, label="Data") - ax_re.plot(delays, fit_curve_re, label="Fit") - ax_re.legend() - - fig_imag, ax_imag = plt.subplots() - ax_imag.set_title(f"{fig_title} - Imaginary") - ax_imag.set_xlabel("Delay (μs)") - ax_imag.set_ylabel("Signal Imaginary (A.U)") - ax_imag.plot(delays, signal_imag, label="Data") - ax_imag.plot(delays, fit_curve_imag, label="Fit") - ax_imag.legend() - - fig_mag, ax_mag = plt.subplots() - ax_mag.set_title(f"{fig_title} - Magnitude") - ax_mag.set_xlabel("Delay (μs)") - ax_mag.set_ylabel("Signal Magnitude (A.U)") - ax_mag.plot(delays, signal_mag, label="Data") - ax_mag.plot(delays, fit_curve_mag, label="Fit") - ax_mag.legend() - - return ( - (fit_result_re, residuals_re, snr_re), - (fit_result_imag, residuals_imag, snr_imag), - (fit_result_mag, residuals_mag, snr_mag), - fig_re, - fig_imag, - fig_mag - ) + self.dependents["signal"] = data["signal"].values + + def _fit_exponentially_decaying_sine(self, delays, signal, fig_title="") -> tuple: + fit = ExponentiallyDecayingSine(delays, signal) + fit_result = fit.run(fit) + fit_curve = fit_result.eval() + residuals = signal - fit_curve + amp = fit_result.params["A"].value + noise = np.std(residuals) + snr = np.abs(amp / (4 * noise)) + + fig, ax = plt.subplots() + ax.set_title(fig_title) + if self.platform_type == PlatformTypes.OPX: + ax.set_xlabel("Delay (ns)") + else: + ax.set_xlabel("Delay (μs)") + ax.set_ylabel("Rotated Signal (A.U)") + ax.plot(delays, signal, label="Data") + ax.plot(delays, fit_curve, label="Fit") + ax.legend() + + return fit_result, residuals, snr, fig def analyze(self): with DatasetAnalysis(self.data_loc, self.name) as ds: - result_re, result_imag, result_mag, fig_re, fig_imag, fig_mag = self._fit_exponentially_decaying_sine_components( + self.fit_result, self.residuals, self.snr, fig = self._fit_exponentially_decaying_sine( self.independents["delays"], self.dependents["signal"], "T2 Echo Measurement" ) - self.fit_result_re, residuals_re, self.snr_re = result_re - self.fit_result_imag, residuals_imag, self.snr_imag = result_imag - self.fit_result_mag, residuals_mag, self.snr_mag = result_mag - # Save all fit results ds.add( - fit_result_re=self.fit_result_re, - params_re=serialize_fit_params(self.fit_result_re.params), - snr_re=float(self.snr_re), - fit_result_imag=self.fit_result_imag, - params_imag=serialize_fit_params(self.fit_result_imag.params), - snr_imag=float(self.snr_imag), - fit_result_mag=self.fit_result_mag, - params_mag=serialize_fit_params(self.fit_result_mag.params), - snr_mag=float(self.snr_mag) + fit_result=self.fit_result, + params=serialize_fit_params(self.fit_result.params), + snr=float(self.snr) ) - # Save all three figures separately - ds.add_figure(f"{self.name}_real", fig=fig_re) - image_path_re = ds._new_file_path(ds.savefolders[1], f"{self.name}_real", suffix="png") - self.figure_paths.append(image_path_re) - - ds.add_figure(f"{self.name}_imag", fig=fig_imag) - image_path_imag = ds._new_file_path(ds.savefolders[1], f"{self.name}_imag", suffix="png") - self.figure_paths.append(image_path_imag) - - ds.add_figure(f"{self.name}_mag", fig=fig_mag) - image_path_mag = ds._new_file_path(ds.savefolders[1], f"{self.name}_mag", suffix="png") - self.figure_paths.append(image_path_mag) - - snr_dict = { - "Real": (self.snr_re, self.fit_result_re, "re"), - "Imaginary": (self.snr_imag, self.fit_result_imag, "imag"), - "Magnitude": (self.snr_mag, self.fit_result_mag, "mag"), - } - self._sorted_components = sorted(snr_dict.items(), key=lambda x: x[1][0], reverse=True) + ds.add_figure(self.name, fig=fig) + image_path = ds._new_file_path(ds.savefolders[1], self.name, suffix="png") + self.figure_paths.append(image_path) def _check_quality(self) -> CheckResult: snr_min = self.snr_min_threshold() - - valid = [ - (name, snr, fit, key) - for name, (snr, fit, key) in self._sorted_components - if snr >= snr_min - ] - - if valid: - self._winner_name, self._winner_snr, self._winner_fit, self._winner_key = valid[0] - max_error = self.max_fit_param_error() - bad_params = [] - for pname, param in self._winner_fit.params.items(): - if param.stderr is None: - bad_params.append(f"{pname}(no stderr)") - elif param.value == 0 or abs(param.stderr / param.value) > max_error: - pct = abs(param.stderr / param.value) * 100 if param.value != 0 else float("inf") - bad_params.append(f"{pname}({pct:.0f}%)") - passed = len(bad_params) == 0 - parts = [f"winner={self._winner_name}, SNR={self._winner_snr:.3f} (threshold={snr_min:.1f})"] - if bad_params: - parts.append(f"high-error params: {', '.join(bad_params)}") - else: - self._winner_name, (self._winner_snr, self._winner_fit, self._winner_key) = self._sorted_components[0] - passed = False - parts = [ - f"no component with SNR >= {snr_min:.1f}", - f"best={self._winner_name}, SNR={self._winner_snr:.3f}", - ] + max_error = self.max_fit_param_error() + bad_params = [] + for pname, param in self.fit_result.params.items(): + if param.stderr is None: + bad_params.append(f"{pname}(no stderr)") + elif param.value == 0 or abs(param.stderr / param.value) > max_error: + pct = abs(param.stderr / param.value) * 100 if param.value != 0 else float("inf") + bad_params.append(f"{pname}({pct:.0f}%)") + passed = self.snr >= snr_min and len(bad_params) == 0 + parts = [f"SNR={self.snr:.3f} (threshold={snr_min:.1f})"] + if bad_params: + parts.append(f"high-error params: {', '.join(bad_params)}") return CheckResult("quality_check", passed, "; ".join(parts)) def correct(self, result: EvaluateResult) -> EvaluateResult: - # Pop all figures before super() auto-appends the last one - fig_re = self.figure_paths.pop(0) if len(self.figure_paths) >= 3 else None - fig_imag = self.figure_paths.pop(0) if self.figure_paths else None - fig_mag = self.figure_paths.pop(0) if self.figure_paths else None - self.figure_paths.clear() # prevent auto-append - - plot_map = {"re": fig_re, "imag": fig_imag, "mag": fig_mag} + figure = self.figure_paths[0] if self.figure_paths else None + self.figure_paths.clear() snr_min = self.snr_min_threshold() self.report_output.append( @@ -435,15 +277,13 @@ def correct(self, result: EvaluateResult) -> EvaluateResult: f"Data Path: `{self.data_loc}`\n\n" ) - for i, (comp_name, (comp_snr, comp_fit, comp_key)) in enumerate(self._sorted_components): - tag = "(SELECTED)" if i == 0 else "(NOT SELECTED)" - self.report_output.append(f"### **{comp_name} Component {tag}**\n") - if plot_map.get(comp_key): - self.report_output.append(plot_map[comp_key]) - self.report_output.append( - f"SNR={comp_snr:.3f}\n\n" - f"**Fit Report:**\n```\n{str(comp_fit.lmfit_result.fit_report())}\n```\n\n" - ) + self.report_output.append("### Rotated Signal Fit\n") + if figure: + self.report_output.append(figure) + self.report_output.append( + f"SNR={self.snr:.3f}\n\n" + f"**Fit Report:**\n```\n{str(self.fit_result.lmfit_result.fit_report())}\n```\n\n" + ) result = super().correct(result) # adds check table + success update line return result diff --git a/src/cqedtoolbox/protocols/operations/single_qubit/t2r.py b/src/cqedtoolbox/protocols/operations/single_qubit/t2r.py index e830a01..2e3f3bb 100644 --- a/src/cqedtoolbox/protocols/operations/single_qubit/t2r.py +++ b/src/cqedtoolbox/protocols/operations/single_qubit/t2r.py @@ -16,7 +16,7 @@ from labcore.protocols.base import ( ProtocolOperation, serialize_fit_params, - CorrectionParameter, CheckResult, Correction, EvaluateResult, + CorrectionParameter, CheckResult, Correction, EvaluateResult, PlatformTypes ) from cqedtoolbox.protocols.parameters import ( Repetition, @@ -25,10 +25,10 @@ ReadoutGain, ReadoutLength, T2R, - NEchos ) from cqedtoolbox.measurement_lib.opx.advanced.qubit_tuneup import measure_t2 from cqedtoolbox.measurement_lib.qick.single_transmon_v2 import T2RProgram +from cqedtoolbox.readout.qubit_readout import rotate_complex_qubit_data logger = logging.getLogger(__name__) @@ -138,7 +138,6 @@ def __init__(self, params): qubit_gain=QubitGain(params), readout_gain=ReadoutGain(params), readout_length=ReadoutLength(params), - n_echos=NEchos(params) ) self._register_outputs( t2r=T2R(params) @@ -159,22 +158,14 @@ def __init__(self, params): self._register_check("quality_check", self._check_quality, self._increase_averaging) - self._register_success_update(self.t2r, lambda: self._winner_fit.params["tau"].value) + self._register_success_update(self.t2r, lambda: self.fit_result.params["tau"].value) self.independents = {"delays": []} self.dependents = {"signal": []} - self.fit_result_re = None - self.fit_result_imag = None - self.fit_result_mag = None - self.snr_re = None - self.snr_imag = None - self.snr_mag = None - self._winner_fit = None - self._winner_snr = None - self._winner_key = None - self._winner_name = None - self._sorted_components = None + self.fit_result = None + self.residuals = None + self.snr = None def _measure_dummy(self) -> Path: logger.info("Starting dummy T2 Ramsey measurement") @@ -187,12 +178,10 @@ def _measure_dummy(self) -> Path: return loc def _load_data_dummy(self): - path = self.data_loc / "data.ddh5" - if not path.exists(): - raise FileNotFoundError(f"File {path} does not exist") - data = datadict_from_hdf5(path) - self.independents["delays"] = data["delays"]["values"] - self.dependents["signal"] = data["signal"]["values"] + data = load_as_xr(self.data_loc) + rotated = rotate_complex_qubit_data(data)[0] + self.independents["delays"] = rotated["delays"].values + self.dependents["signal"] = rotated["signal"].values def _measure_qick(self) -> Path: logger.info("Starting qick T2 Ramsey measurement") @@ -209,175 +198,80 @@ def _measure_opx(self) -> Path: return loc def _load_data_qick(self): - path = self.data_loc / "data.ddh5" - if not path.exists(): - raise FileNotFoundError(f"File {path} does not exist") - data = datadict_from_hdf5(path) - - self.independents["delays"] = data["t"]["values"] - self.dependents["signal"] = data["signal"]["values"] + data = load_as_xr(self.data_loc) + rotated = rotate_complex_qubit_data(data)[0] + self.independents["delays"] = rotated["t"].values + self.dependents["signal"] = rotated["signal"].values def _load_data_opx(self): - data = load_as_xr(self.data_loc).mean("repetition") + data = load_as_xr(self.data_loc) + if "repetition" in data.dims: + data = data.mean("repetition") + data, _ = rotate_complex_qubit_data(data) self.independents["delays"] = data["delay"].values - self.dependents["signal"] = data["signal_Re"].values + 1j * data["signal_Im"].values - - def _fit_exponentially_decaying_sine_components(self, delays, signal, fig_title="") -> tuple: - """ - Fit real, imaginary, and magnitude components with ExponentiallyDecayingSine fits. - Returns (fit_result_re, fit_result_imag, fit_result_mag, fig_re, fig_imag, fig_mag) - """ - signal_re = signal.real - signal_imag = signal.imag - signal_mag = np.abs(signal) - - # Fit real part - fit_re = ExponentiallyDecayingSine(delays, signal_re) - fit_result_re = fit_re.run(fit_re) - fit_curve_re = fit_result_re.eval() - residuals_re = signal_re - fit_curve_re - amp_re = fit_result_re.params["A"].value - noise_re = np.std(residuals_re) - snr_re = np.abs(amp_re / (4 * noise_re)) - - # Fit imaginary part - fit_imag = ExponentiallyDecayingSine(delays, signal_imag) - fit_result_imag = fit_imag.run(fit_imag) - fit_curve_imag = fit_result_imag.eval() - residuals_imag = signal_imag - fit_curve_imag - amp_imag = fit_result_imag.params["A"].value - noise_imag = np.std(residuals_imag) - snr_imag = np.abs(amp_imag / (4 * noise_imag)) - - # Fit magnitude - fit_mag = ExponentiallyDecayingSine(delays, signal_mag) - fit_result_mag = fit_mag.run(fit_mag) - fit_curve_mag = fit_result_mag.eval() - residuals_mag = signal_mag - fit_curve_mag - amp_mag = fit_result_mag.params["A"].value - noise_mag = np.std(residuals_mag) - snr_mag = np.abs(amp_mag / (4 * noise_mag)) - - # Create three separate figures - fig_re, ax_re = plt.subplots() - ax_re.set_title(f"{fig_title} - Real") - ax_re.set_xlabel("Delay (μs)") - ax_re.set_ylabel("Signal Real (A.U)") - ax_re.plot(delays, signal_re, label="Data") - ax_re.plot(delays, fit_curve_re, label="Fit") - ax_re.legend() - - fig_imag, ax_imag = plt.subplots() - ax_imag.set_title(f"{fig_title} - Imaginary") - ax_imag.set_xlabel("Delay (μs)") - ax_imag.set_ylabel("Signal Imaginary (A.U)") - ax_imag.plot(delays, signal_imag, label="Data") - ax_imag.plot(delays, fit_curve_imag, label="Fit") - ax_imag.legend() - - fig_mag, ax_mag = plt.subplots() - ax_mag.set_title(f"{fig_title} - Magnitude") - ax_mag.set_xlabel("Delay (μs)") - ax_mag.set_ylabel("Signal Magnitude (A.U)") - ax_mag.plot(delays, signal_mag, label="Data") - ax_mag.plot(delays, fit_curve_mag, label="Fit") - ax_mag.legend() - - return ( - (fit_result_re, residuals_re, snr_re), - (fit_result_imag, residuals_imag, snr_imag), - (fit_result_mag, residuals_mag, snr_mag), - fig_re, - fig_imag, - fig_mag - ) + self.dependents["signal"] = data["signal"].values + + def _fit_exponentially_decaying_sine(self, delays, signal, fig_title="") -> tuple: + fit = ExponentiallyDecayingSine(delays, signal) + fit_result = fit.run(fit) + fit_curve = fit_result.eval() + residuals = signal - fit_curve + amp = fit_result.params["A"].value + noise = np.std(residuals) + snr = np.abs(amp / (4 * noise)) + + fig, ax = plt.subplots() + ax.set_title(fig_title) + if self.platform_type == PlatformTypes.OPX: + ax.set_xlabel("Delay (ns)") + else: + ax.set_xlabel("Delay (μs)") + ax.set_ylabel("Rotated Signal (A.U)") + ax.plot(delays, signal, label="Data") + ax.plot(delays, fit_curve, label="Fit") + ax.legend() + + return fit_result, residuals, snr, fig def analyze(self): with DatasetAnalysis(self.data_loc, self.name) as ds: - result_re, result_imag, result_mag, fig_re, fig_imag, fig_mag = self._fit_exponentially_decaying_sine_components( + self.fit_result, self.residuals, self.snr, fig = self._fit_exponentially_decaying_sine( self.independents["delays"], self.dependents["signal"], "T2 Ramsey Measurement" ) - self.fit_result_re, residuals_re, self.snr_re = result_re - self.fit_result_imag, residuals_imag, self.snr_imag = result_imag - self.fit_result_mag, residuals_mag, self.snr_mag = result_mag - # Save all fit results ds.add( - fit_result_re=self.fit_result_re, - params_re=serialize_fit_params(self.fit_result_re.params), - snr_re=float(self.snr_re), - fit_result_imag=self.fit_result_imag, - params_imag=serialize_fit_params(self.fit_result_imag.params), - snr_imag=float(self.snr_imag), - fit_result_mag=self.fit_result_mag, - params_mag=serialize_fit_params(self.fit_result_mag.params), - snr_mag=float(self.snr_mag) + fit_result=self.fit_result, + params=serialize_fit_params(self.fit_result.params), + snr=float(self.snr) ) - # Save all three figures separately - ds.add_figure(f"{self.name}_real", fig=fig_re) - image_path_re = ds._new_file_path(ds.savefolders[1], f"{self.name}_real", suffix="png") - self.figure_paths.append(image_path_re) - - ds.add_figure(f"{self.name}_imag", fig=fig_imag) - image_path_imag = ds._new_file_path(ds.savefolders[1], f"{self.name}_imag", suffix="png") - self.figure_paths.append(image_path_imag) - - ds.add_figure(f"{self.name}_mag", fig=fig_mag) - image_path_mag = ds._new_file_path(ds.savefolders[1], f"{self.name}_mag", suffix="png") - self.figure_paths.append(image_path_mag) - - snr_dict = { - "Real": (self.snr_re, self.fit_result_re, "re"), - "Imaginary": (self.snr_imag, self.fit_result_imag, "imag"), - "Magnitude": (self.snr_mag, self.fit_result_mag, "mag"), - } - self._sorted_components = sorted(snr_dict.items(), key=lambda x: x[1][0], reverse=True) + ds.add_figure(self.name, fig=fig) + image_path = ds._new_file_path(ds.savefolders[1], self.name, suffix="png") + self.figure_paths.append(image_path) def _check_quality(self) -> CheckResult: snr_min = self.snr_min_threshold() - - valid = [ - (name, snr, fit, key) - for name, (snr, fit, key) in self._sorted_components - if snr >= snr_min - ] - - if valid: - self._winner_name, self._winner_snr, self._winner_fit, self._winner_key = valid[0] - max_error = self.max_fit_param_error() - bad_params = [] - for pname, param in self._winner_fit.params.items(): - if param.stderr is None: - bad_params.append(f"{pname}(no stderr)") - elif param.value == 0 or abs(param.stderr / param.value) > max_error: - pct = abs(param.stderr / param.value) * 100 if param.value != 0 else float("inf") - bad_params.append(f"{pname}({pct:.0f}%)") - passed = len(bad_params) == 0 - parts = [f"winner={self._winner_name}, SNR={self._winner_snr:.3f} (threshold={snr_min:.1f})"] - if bad_params: - parts.append(f"high-error params: {', '.join(bad_params)}") - else: - self._winner_name, (self._winner_snr, self._winner_fit, self._winner_key) = self._sorted_components[0] - passed = False - parts = [ - f"no component with SNR >= {snr_min:.1f}", - f"best={self._winner_name}, SNR={self._winner_snr:.3f}", - ] + max_error = self.max_fit_param_error() + bad_params = [] + for pname, param in self.fit_result.params.items(): + if param.stderr is None: + bad_params.append(f"{pname}(no stderr)") + elif param.value == 0 or abs(param.stderr / param.value) > max_error: + pct = abs(param.stderr / param.value) * 100 if param.value != 0 else float("inf") + bad_params.append(f"{pname}({pct:.0f}%)") + passed = self.snr >= snr_min and len(bad_params) == 0 + parts = [f"SNR={self.snr:.3f} (threshold={snr_min:.1f})"] + if bad_params: + parts.append(f"high-error params: {', '.join(bad_params)}") return CheckResult("quality_check", passed, "; ".join(parts)) def correct(self, result: EvaluateResult) -> EvaluateResult: - # Pop all figures before super() auto-appends the last one - fig_re = self.figure_paths.pop(0) if len(self.figure_paths) >= 3 else None - fig_imag = self.figure_paths.pop(0) if self.figure_paths else None - fig_mag = self.figure_paths.pop(0) if self.figure_paths else None - self.figure_paths.clear() # prevent auto-append - - plot_map = {"re": fig_re, "imag": fig_imag, "mag": fig_mag} + figure = self.figure_paths[0] if self.figure_paths else None + self.figure_paths.clear() snr_min = self.snr_min_threshold() self.report_output.append( @@ -386,15 +280,13 @@ def correct(self, result: EvaluateResult) -> EvaluateResult: f"Data Path: `{self.data_loc}`\n\n" ) - for i, (comp_name, (comp_snr, comp_fit, comp_key)) in enumerate(self._sorted_components): - tag = "(SELECTED)" if i == 0 else "(NOT SELECTED)" - self.report_output.append(f"### **{comp_name} Component {tag}**\n") - if plot_map.get(comp_key): - self.report_output.append(plot_map[comp_key]) - self.report_output.append( - f"SNR={comp_snr:.3f}\n\n" - f"**Fit Report:**\n```\n{str(comp_fit.lmfit_result.fit_report())}\n```\n\n" - ) + self.report_output.append("### Rotated Signal Fit\n") + if figure: + self.report_output.append(figure) + self.report_output.append( + f"SNR={self.snr:.3f}\n\n" + f"**Fit Report:**\n```\n{str(self.fit_result.lmfit_result.fit_report())}\n```\n\n" + ) result = super().correct(result) # adds check table + success update line return result diff --git a/src/cqedtoolbox/protocols/parameters.py b/src/cqedtoolbox/protocols/parameters.py index 2f92fe6..ee8b90b 100644 --- a/src/cqedtoolbox/protocols/parameters.py +++ b/src/cqedtoolbox/protocols/parameters.py @@ -283,6 +283,14 @@ def _qick_getter(self): def _qick_setter(self, value): active_qubit = nestedAttributeFromString(self.params, "active.qubit")() return nestedAttributeFromString(self.params, f"{active_qubit}.qubit.chi")(value) + + def _opx_getter(self): + active_qubit = nestedAttributeFromString(self.params, "active.qubit")() + return nestedAttributeFromString(self.params, f"{active_qubit}.chi")() + + def _opx_setter(self, value): + active_qubit = nestedAttributeFromString(self.params, "active.qubit")() + return nestedAttributeFromString(self.params, f"{active_qubit}.chi")(value) @dataclass diff --git a/src/cqedtoolbox/protocols/qubit_tuneup.py b/src/cqedtoolbox/protocols/qubit_tuneup.py index 82d01d2..c48b6a8 100644 --- a/src/cqedtoolbox/protocols/qubit_tuneup.py +++ b/src/cqedtoolbox/protocols/qubit_tuneup.py @@ -6,17 +6,17 @@ class QubitTuneup(ProtocolBase): - def __init__(self, params, report_path: Path = Path(".")): + def __init__(self, params, geometry, report_path: Path = Path(".")): super().__init__(report_path) self.root_branch = BranchBase("QubitTuneup") self.root_branch.extend([ - ResonatorSpectroscopy(params), - ResonatorSpectroscopyVsGain(params), + ResonatorSpectroscopy(params, geometry), + ResonatorSpectroscopyVsGain(params, geometry), SaturationSpectroscopy(params), PowerRabi(params), PiSpectroscopy(params), - ResonatorSpectroscopyAfterPi(params), + ResonatorSpectroscopyAfterPi(params, geometry), T1Operation(params), T2ROperation(params), T2EOperation(params),