# Repository: https://gitlab.com/quantify-os/quantify-core
# Licensed according to the LICENCE file on the main branch
"""Module containing an analysis class for two-state readout calibration."""
from __future__ import annotations
import math
import lmfit
import matplotlib.pyplot as plt
import numpy as np
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
import quantify_core.data.handling as dh
from quantify_core.analysis import base_analysis as ba
from quantify_core.visualization import mpl_plotting as qpl
from quantify_core.visualization.SI_utilities import format_value_string
[docs]
class ReadoutCalibrationAnalysis(ba.BaseAnalysis):
# pylint: disable=line-too-long
"""
Find threshold and angle which discriminates qubit state.
.. admonition:: Example
.. jupyter-execute::
import os
import quantify_core.data.handling as dh
from quantify_core.analysis.readout_calibration_analysis import (
ReadoutCalibrationAnalysis,
)
# load example data
test_data_dir = "../tests/test_data"
dh.set_datadir(test_data_dir)
ReadoutCalibrationAnalysis(tuid="20230509-152441-841-faef49").run().display_figs_mpl()
"""
# pylint: disable=no-member
# pylint: disable=attribute-defined-outside-init
# pylint: disable=broad-exception-caught
[docs]
def process_data(self) -> None:
"""Process the data so that the analysis can make assumptions on the format."""
self.dataset_processed = dh.to_gridded_dataset(self.dataset)
[docs]
def run_fitting(self) -> None:
"""Fit a state discriminator to the readout calibration data."""
points = np.zeros((self.dataset_processed.x0.data.size, 2), dtype=float)
points[:, 0] = self.dataset_processed.y0.data # real (I)
points[:, 1] = self.dataset_processed.y1.data # imag (Q)
# make a dummy model just to store the params with analysis infrastructure
lin_model = lmfit.models.LinearModel()
result = lin_model.fit([0, 1], x=[0, 1])
# overwrite parameters
result.params = lmfit.Parameters()
try:
self._lda = LinearDiscriminantAnalysis(solver="svd", store_covariance=True)
self._lda.fit(points, self.dataset_processed.x0.data)
intercept = -self._lda.intercept_[0] / self._lda.coef_[0][1]
slope = -self._lda.coef_[0][0] / self._lda.coef_[0][1]
result.params.add("intercept", value=intercept, vary=False)
result.params.add("slope", value=slope, vary=False)
result.params.add(
"acq_threshold",
# compute distance from origin to line
value=abs(intercept) / math.sqrt(slope**2 + 1),
vary=False,
)
result.params.add(
"acq_rotation_rad",
# compute CW angle from y-axis
value=-math.pi / 2 - math.atan(slope),
vary=False,
)
result.success = True
result.errorbars = { # type: ignore
"intercept": False,
"slope": False,
"acq_threshold": False,
"acq_rotation_rad": False,
}
except Exception as lda_error:
result.success = False
result.message = "Error during fit:\n" + str(lda_error) # type: ignore
self.fit_results.update({"linear_discriminator": result})
def _get_points(self) -> tuple:
points = np.zeros((self.dataset_processed.x0.data.size, 2), dtype=float)
points[:, 0] = self.dataset_processed.y0.data # real (I)
points[:, 1] = self.dataset_processed.y1.data # imag (Q)
y_pred = self._lda.predict(points)
true_pos = self.dataset_processed.x0.data == y_pred
true_pos0, true_pos1 = (
true_pos[self.dataset_processed.x0.data == 0],
true_pos[self.dataset_processed.x0.data == 1],
)
x0, x1 = (
points[self.dataset_processed.x0.data == 0],
points[self.dataset_processed.x0.data == 1],
)
return (x0[true_pos0], x0[~true_pos0], x1[true_pos1], x1[~true_pos1])
[docs]
def analyze_fit_results(self) -> None:
"""Check the fit success and populate :code:`.quantities_of_interest`."""
self.quantities_of_interest = {
"fit_success": self.fit_results["linear_discriminator"].success
}
if not self.quantities_of_interest["fit_success"]:
self.quantities_of_interest["fit_msg"] = self.fit_results[
"linear_discriminator"
].message
return
# Set fitted parameters as quantities of interest
for parameter, value in self.fit_results["linear_discriminator"].params.items():
self.quantities_of_interest[parameter] = ba.lmfit_par_to_ufloat(value)
# Get centroid from LDA model
self.quantities_of_interest["avg_ket0"] = (
self._lda.means_[0][0] + self._lda.means_[0][1] * 1j
)
self.quantities_of_interest["avg_ket1"] = (
self._lda.means_[1][0] + self._lda.means_[1][1] * 1j
)
# Compute fidelity estimates
(x0_tp, x0_fp, x1_tp, x1_fp) = self._get_points()
self.quantities_of_interest["fid_est_0"] = x0_tp.shape[0] / (
x0_tp.shape[0] + x0_fp.shape[0]
)
self.quantities_of_interest["fid_est_1"] = x1_tp.shape[0] / (
x1_tp.shape[0] + x1_fp.shape[0]
)
text_msg = "Sample means:\n"
text_msg += format_value_string(
r"$|0\rangle$",
self.quantities_of_interest["avg_ket0"],
unit=self.dataset_processed.y0.units,
end_char="\n",
)
text_msg += format_value_string(
r"$|1\rangle$",
self.quantities_of_interest["avg_ket1"],
unit=self.dataset_processed.y0.units,
end_char="\n\n",
)
text_msg += "Discriminator:\n"
text_msg += format_value_string(
"Threshold",
self.quantities_of_interest["acq_threshold"],
unit=self.dataset_processed.y0.units,
end_char="\n",
)
text_msg += format_value_string(
r"$\theta$ from y-axis",
self.quantities_of_interest["acq_rotation_rad"],
unit="rad",
end_char="\n\n",
)
text_msg += "Fidelity estimates:\n"
text_msg += format_value_string(
r"$|0\rangle$",
self.quantities_of_interest["fid_est_0"] * 100,
unit="%",
end_char="\n",
)
text_msg += format_value_string(
r"$|1\rangle$",
self.quantities_of_interest["fid_est_1"] * 100,
unit="%",
end_char="\n",
)
self.quantities_of_interest["fit_msg"] = text_msg
