From ab40c5a4fe7a4b2867e5468d80787327545e4534 Mon Sep 17 00:00:00 2001 From: Naoki Kanazawa Date: Thu, 11 Jan 2024 07:20:08 +0900 Subject: [PATCH] Replace the fitter of StarkRamseyXYAmpScan experiment (#1243) ### Summary This PR replaces the fitter of `StarkRamseyXYAmpScan` experiment. Current fitter puts strong assumption on model parameter and it makes the fitting quite unstable to the change in dephasing rate. This PR doesn't introduce any breaking API change and a separate release note is not necessary because the experiment is not released yet. ### Details and comments The current fitter simultaneously fits the P1 values of the Ramsey X and Y experiment. Given the Stark shift is the third order polynominal of tone amplitude x, the P1 values for X and Y are ``` Px(x) = A cos(2 pi t_S (c1 x + c2 x**2 + c3 x**3 + f_err)) + offset Px(x) = A sin(2 pi t_S (c1 x + c2 x**2 + c3 x**3 + f_err)) + offset ``` where t_S is the fixed Stark delay and others are the fit parameters. Note that current fitter assumes that A is independent of x, however, since this is a variant of Ramsey experiment, the P1 amplitude actually depends on x. This is because qubit T2 may depend on the frequency, or at larger Stark amplitude the qubit heating may cause faster dephasing. In the new fitter, we directly extract the phase polynominal and perform fit on this synthesized data. Namely, ``` poly = unwrap(arctan2(Py, Px)) ~ 2 pi t_S (c1 x + c2 x**2 + c3 x**3 + f_err) ``` Because A is canceled out in this form, the fitting becomes robust to the amplitude dependence of the dephasing rate. Usually we cannot obtain the reasonable model for A(x). For example, the experiment data below shows the difference of fit on the poly data and raw Px, Py data (actually the parameters are fit on the poly data and Px, Py data are just visualized with the fitted parameters). ![image](https://github.com/Qiskit-Extensions/qiskit-experiments/assets/39517270/d14b8537-dcc9-489c-b277-fe7df7af8938) As you can see, the Px, Py data are damped quickly on positive amplitudes, and the envelope of the trigonometric function is asymmetric with respect to 0. In this situation the current fitter doesn't work well. In contrast, new phase fitter works nicely. --------- Co-authored-by: Will Shanks Co-authored-by: Will Shanks --- .../curve_analysis/scatter_table.py | 2 +- .../analysis/ramsey_xy_analysis.py | 410 +++++++++++------- .../characterization/test_stark_ramsey_xy.py | 62 ++- 3 files changed, 297 insertions(+), 177 deletions(-) diff --git a/qiskit_experiments/curve_analysis/scatter_table.py b/qiskit_experiments/curve_analysis/scatter_table.py index 8f573226a7..4361274b9e 100644 --- a/qiskit_experiments/curve_analysis/scatter_table.py +++ b/qiskit_experiments/curve_analysis/scatter_table.py @@ -127,7 +127,7 @@ def y_err(self) -> np.ndarray: ) def shots(self): """Shot number of data points.""" - return self.shots.to_numpy() + return self["shots"].to_numpy() @property @deprecate_func( diff --git a/qiskit_experiments/library/characterization/analysis/ramsey_xy_analysis.py b/qiskit_experiments/library/characterization/analysis/ramsey_xy_analysis.py index cf8d74515b..1a0ac92837 100644 --- a/qiskit_experiments/library/characterization/analysis/ramsey_xy_analysis.py +++ b/qiskit_experiments/library/characterization/analysis/ramsey_xy_analysis.py @@ -16,6 +16,7 @@ import lmfit import numpy as np +from uncertainties import unumpy as unp import qiskit_experiments.curve_analysis as curve import qiskit_experiments.visualization as vis @@ -215,19 +216,48 @@ class StarkRamseyXYAmpScanAnalysis(curve.CurveAnalysis): # section: overview - This analysis is a variant of :class:`RamseyXYAnalysis` in which - the data is fit for a trigonometric function model with a linear phase. - By contrast, in this model, the phase is assumed to be a polynomial of the x-data, - and techniques to compute a good initial guess for these polynomial coefficients - are not trivial. For example, when the phase is a linear function of the x-data, - one may apply a Fourier transform to the data to estimate the coefficient, - but this technique can not be used for a higher order polynomial. + This analysis is a variant of :class:`RamseyXYAnalysis`. In both cases, the X and Y + data are treated as the real and imaginary parts of a complex oscillating signal. + In :class:`RamseyXYAnalysis`, the data are fit assuming a phase varying linearly with + the x-data corresponding to a constant frequency and assuming an exponentially + decaying amplitude. By contrast, in this model, the phase is assumed to be + a third order polynomial :math:`\theta(x)` of the x-data. + Additionally, the amplitude is not assumed to follow a specific form. + Techniques to compute a good initial guess for the polynomial coefficients inside + a trigonometric function like this are not trivial. Instead, this analysis extracts the + raw phase and runs fits on the extracted data to a polynomial :math:`\theta(x)` directly. - This analysis assumes the following polynomial for the phase imparted by the Stark shift. + The measured P1 values for a Ramsey X and Y experiment can be written in the form of + a trignometric function taking the phase polynomial :math:`\theta(x)`: .. math:: - \theta_{\text Stark}(x) = 2 \pi t_S f_S(x), + P_X = \text{amp}(x) \cdot \cos \theta(x) + \text{offset},\\ + P_Y = \text{amp}(x) \cdot \sin \theta(x) + \text{offset}. + + Hence the phase polynomial can be extracted as follows + + .. math:: + + \theta(x) = \tan^{-1} \frac{P_Y}{P_X}. + + Because the arctangent is implemented by the ``atan2`` function + defined in :math:`[-\pi, \pi]`, the computed :math:`\theta(x)` is unwrapped to + ensure continuous phase evolution. + + We call attention to the fact that :math:`\text{amp}(x)` is also Stark tone amplitude + dependent because of the qubit frequency dependence of the dephasing rate. + In general :math:`\text{amp}(x)` is unpredictable due to dephasing from + two-level systems distributed randomly in frequency + or potentially due to qubit heating. This prevents us from precisely fitting + the raw :math:`P_X`, :math:`P_Y` data. Fitting only the phase data makes the + analysis robust to amplitude dependent dephasing. + + In this analysis, the phase polynomial is defined as + + .. math:: + + \theta(x) = 2 \pi t_S f_S(x) where @@ -245,73 +275,21 @@ class StarkRamseyXYAmpScanAnalysis(curve.CurveAnalysis): .. math:: - F_{X+} = \text{amp} \cdot \cos \left( 2 \pi t_S f_S^+(x) \right) + \text{offset}, \\ - F_{Y+} = \text{amp} \cdot \sin \left( 2 \pi t_S f_S^+(x) \right) + \text{offset}, \\ - F_{X-} = \text{amp} \cdot \cos \left( 2 \pi t_S f_S^-(x) \right) + \text{offset}, \\ - F_{Y-} = \text{amp} \cdot \sin \left( 2 \pi t_S f_S^-(x) \right) + \text{offset}, - - where - - .. math :: - - f_S^\nu(x) = c_1^\nu x + c_2^\nu x^2 + c_3^\nu x^3 + f_{\rm err}. + \theta^\nu(x) = c_1^\nu x + c_2^\nu x^2 + c_3^\nu x^3 + f_{\rm err}, + where :math:`\nu \in \{+, -\}`. The Stark shift is asymmetric with respect to :math:`x=0`, because of the anti-crossings of higher energy levels. In a typical transmon qubit, these levels appear only in :math:`f_S < 0` because of the negative anharmonicity. To precisely fit the results, this analysis uses different model parameters for positive (:math:`x > 0`) and negative (:math:`x < 0`) shift domains. - To obtain the initial guess, the following calculation is employed in this analysis. - First, oscillations in each quadrature are normalized and the offset is subtracted. - The amplitude and offset can be accurately estimated from the experiment data - when the oscillation involves multiple cycles. - - .. math :: - - {\cal F}_{X} = \cos \left( 2 \pi t_S f_S(x) \right), \\ - {\cal F}_{Y} = \sin \left( 2 \pi t_S f_S(x) \right). - - Next, these normalized oscillations are differentiated - - .. math :: - - \dot{{\cal F}}_X = - 2 \pi t_S \frac{d f_S}{dx} {\cal F}_Y, \\ - \dot{{\cal F}}_Y = 2 \pi t_S \frac{d f_S}{dx} {\cal F}_X. \\ - - The square root of the sum of the squares of the above quantities yields - - .. math :: - - \sqrt{\dot{{\cal F}}_X^2 + \dot{{\cal F}}_Y^2} - = 2 \pi t_S \frac{d}{dx} f_S = 2 \pi t_S (c_1 + 2 c_2 x + 3 c_3 x^2). - - By computing this synthesized data on the left hand side, one can estimate - the initial guess of the polynomial coefficients by quadratic regression. - This fit protocol is independently conducted for the experiment data on the - positive and negative shift domain. - # section: fit_parameters - defpar \rm amp: - desc: Amplitude of both series. - init_guess: Median of root sum square of Ramsey X and Y oscillation. - bounds: [0, 1] - - defpar \rm offset: - desc: Base line of all series. - init_guess: The average of the data. - bounds: [-1, 1] - - defpar t_S: - desc: Fixed parameter from the ``stark_length`` experiment option. - init_guess: Automatically set from metadata when this analysis is run. - bounds: None - defpar c_1^+: desc: The linear term coefficient of the positive Stark shift (fit parameter: ``stark_pos_coef_o1``). - init_guess: See the fit model description. + init_guess: 0. bounds: None defpar c_2^+: @@ -320,19 +298,19 @@ class StarkRamseyXYAmpScanAnalysis(curve.CurveAnalysis): induce blue shift when its sign is positive. Note that the quadratic term is the primary term (fit parameter: ``stark_pos_coef_o2``). - init_guess: See the fit model description. + init_guess: 1e6. bounds: [0, inf] defpar c_3^+: desc: The cubic term coefficient of the positive Stark shift (fit parameter: ``stark_pos_coef_o3``). - init_guess: See the fit model description. + init_guess: 0. bounds: None defpar c_1^-: desc: The linear term coefficient of the negative Stark shift. (fit parameter: ``stark_neg_coef_o1``). - init_guess: See the fit model description. + init_guess: 0. bounds: None defpar c_2^-: @@ -341,13 +319,13 @@ class StarkRamseyXYAmpScanAnalysis(curve.CurveAnalysis): induce red shift when its sign is negative. Note that the quadratic term is the primary term (fit parameter: ``stark_neg_coef_o2``). - init_guess: See the fit model description. + init_guess: -1e6. bounds: [-inf, 0] defpar c_3^-: desc: The cubic term coefficient of the negative Stark shift (fit parameter: ``stark_neg_coef_o3``). - init_guess: See the fit model description. + init_guess: 0. bounds: None defpar f_{\rm err}: @@ -364,44 +342,80 @@ class StarkRamseyXYAmpScanAnalysis(curve.CurveAnalysis): def __init__(self): - models = [] - for direction in ("pos", "neg"): - # Ramsey phase := 2π ts Δf(x); Δf(x) = c1 x + c2 x^2 + c3 x^3 + f_err - fs = f"(c1_{direction} * x + c2_{direction} * x**2 + c3_{direction} * x**3 + f_err)" - models.extend( - [ - lmfit.models.ExpressionModel( - expr=f"amp * cos(2 * pi * ts * {fs}) + offset", - name=f"X{direction}", - ), - lmfit.models.ExpressionModel( - expr=f"amp * sin(2 * pi * ts * {fs}) + offset", - name=f"Y{direction}", - ), - ] - ) - + models = [ + lmfit.models.ExpressionModel( + expr="c1_pos * x + c2_pos * x**2 + c3_pos * x**3 + f_err", + name="FREQpos", + ), + lmfit.models.ExpressionModel( + expr="c1_neg * x + c2_neg * x**2 + c3_neg * x**3 + f_err", + name="FREQneg", + ), + ] super().__init__(models=models) @classmethod def _default_options(cls): - """Default analysis options.""" + """Default analysis options. + + Analysis Options: + pulse_len (float): Duration of effective Stark pulse in units of sec. + """ ramsey_plotter = vis.CurvePlotter(vis.MplDrawer()) ramsey_plotter.set_figure_options( xlabel="Stark tone amplitude", - ylabel="P(1)", - ylim=(0, 1), + ylabel=["Stark shift", "P1"], + yval_unit=["Hz", None], series_params={ - "Xpos": {"color": "#123FE8", "symbol": "o", "label": "Ramsey X(+)"}, - "Ypos": {"color": "#6312E8", "symbol": "^", "label": "Ramsey Y(+)"}, - "Xneg": {"color": "#E83812", "symbol": "o", "label": "Ramsey X(-)"}, - "Yneg": {"color": "#E89012", "symbol": "^", "label": "Ramsey Y(-)"}, + "Fpos": { + "color": "#123FE8", + "symbol": "^", + "label": "", + "canvas": 0, + }, + "Fneg": { + "color": "#123FE8", + "symbol": "v", + "label": "", + "canvas": 0, + }, + "Xpos": { + "color": "#123FE8", + "symbol": "o", + "label": "Ramsey X", + "canvas": 1, + }, + "Ypos": { + "color": "#6312E8", + "symbol": "^", + "label": "Ramsey Y", + "canvas": 1, + }, + "Xneg": { + "color": "#E83812", + "symbol": "o", + "label": "Ramsey X", + "canvas": 1, + }, + "Yneg": { + "color": "#E89012", + "symbol": "^", + "label": "Ramsey Y", + "canvas": 1, + }, }, + sharey=False, ) - ramsey_plotter.set_options(style=vis.PlotStyle({"figsize": (12, 5)})) + ramsey_plotter.set_options(subplots=(2, 1), style=vis.PlotStyle({"figsize": (10, 8)})) options = super()._default_options() options.update_options( + data_subfit_map={ + "Xpos": {"series": "X", "direction": "pos"}, + "Ypos": {"series": "Y", "direction": "pos"}, + "Xneg": {"series": "X", "direction": "neg"}, + "Yneg": {"series": "Y", "direction": "neg"}, + }, result_parameters=[ curve.ParameterRepr("c1_pos", "stark_pos_coef_o1", "Hz"), curve.ParameterRepr("c2_pos", "stark_pos_coef_o2", "Hz"), @@ -411,17 +425,76 @@ def _default_options(cls): curve.ParameterRepr("c3_neg", "stark_neg_coef_o3", "Hz"), curve.ParameterRepr("f_err", "stark_ferr", "Hz"), ], - data_subfit_map={ - "Xpos": {"series": "X", "direction": "pos"}, - "Ypos": {"series": "Y", "direction": "pos"}, - "Xneg": {"series": "X", "direction": "neg"}, - "Yneg": {"series": "Y", "direction": "neg"}, - }, plotter=ramsey_plotter, + fit_category="freq", + pulse_len=None, ) return options + def _freq_phase_coef(self) -> float: + """Return a coefficient to convert frequency into phase value.""" + try: + return 2 * np.pi * self.options.pulse_len + except TypeError as ex: + raise TypeError( + "A float-value duration in units of sec of the Stark pulse must be provided. " + f"The pulse_len option value {self.options.pulse_len} is not valid." + ) from ex + + def _format_data( + self, + curve_data: curve.ScatterTable, + category: str = "freq", + ) -> curve.ScatterTable: + + curve_data = super()._format_data(curve_data, category="ramsey_xy") + ramsey_xy = curve_data[curve_data.category == "ramsey_xy"] + + # Create phase data by arctan(Y/X) + columns = list(curve_data.columns) + phase_data = np.empty((0, len(columns))) + y_mean = ramsey_xy.yval.mean() + + grouped = ramsey_xy.groupby("name") + for m_id, direction in enumerate(("pos", "neg")): + x_quadrature = grouped.get_group(f"X{direction}") + y_quadrature = grouped.get_group(f"Y{direction}") + if not np.array_equal(x_quadrature.xval, y_quadrature.xval): + raise ValueError( + "Amplitude values of X and Y quadrature are different. " + "Same values must be used." + ) + x_uarray = unp.uarray(x_quadrature.yval, x_quadrature.yerr) + y_uarray = unp.uarray(y_quadrature.yval, y_quadrature.yerr) + + amplitudes = x_quadrature.xval.to_numpy() + + # pylint: disable=no-member + phase = unp.arctan2(y_uarray - y_mean, x_uarray - y_mean) + phase_n = unp.nominal_values(phase) + phase_s = unp.std_devs(phase) + + # Unwrap phase + # We assume a smooth slope and correct 2pi phase jump to minimize the change of the slope. + unwrapped_phase = np.unwrap(phase_n) + if amplitudes[0] < 0: + # Preserve phase value closest to 0 amplitude + unwrapped_phase = unwrapped_phase + (phase_n[-1] - unwrapped_phase[-1]) + + # Store new data + tmp = np.empty((len(amplitudes), len(columns)), dtype=object) + tmp[:, columns.index("xval")] = amplitudes + tmp[:, columns.index("yval")] = unwrapped_phase / self._freq_phase_coef() + tmp[:, columns.index("yerr")] = phase_s / self._freq_phase_coef() + tmp[:, columns.index("name")] = f"FREQ{direction}" + tmp[:, columns.index("class_id")] = m_id + tmp[:, columns.index("shots")] = x_quadrature.shots + y_quadrature.shots + tmp[:, columns.index("category")] = category + phase_data = np.r_[phase_data, tmp] + + return curve_data.append_list_values(other=phase_data) + def _generate_fit_guesses( self, user_opt: curve.FitOptions, @@ -436,65 +509,91 @@ def _generate_fit_guesses( Returns: List of fit options that are passed to the fitter function. """ - # Compute offset guess + user_opt.bounds.set_if_empty(c2_pos=(0, np.inf), c2_neg=(-np.inf, 0)) user_opt.p0.set_if_empty( - offset=np.mean(curve_data.y), - f_err=0.0, - ) - user_opt.bounds.set_if_empty( - offset=(-1, 1), - amp=(0, 1), - c2_pos=(0, np.inf), - c2_neg=(-np.inf, 0), + c1_pos=0, c2_pos=1e6, c3_pos=0, c1_neg=0, c2_neg=-1e6, c3_neg=0, f_err=0 ) - est_offs = user_opt.p0["offset"] + return user_opt - # Compute amplitude guess - amps = np.zeros(0) - for direction in ("pos", "neg"): - ram_x_off = curve_data.get_subset_of(f"X{direction}").y - est_offs - ram_y_off = curve_data.get_subset_of(f"Y{direction}").y - est_offs - amps = np.concatenate([amps, np.sqrt(ram_x_off**2 + ram_y_off**2)]) - user_opt.p0.set_if_empty(amp=np.median(amps)) - est_a = user_opt.p0["amp"] - const = 2 * np.pi * user_opt.p0["ts"] - - # Compute polynomial coefficients + def _create_figures( + self, + curve_data: curve.ScatterTable, + ) -> List["matplotlib.figure.Figure"]: + + # plot unwrapped phase on first axis + for d in ("pos", "neg"): + sub_data = curve_data[(curve_data.name == f"FREQ{d}") & (curve_data.category == "freq")] + self.plotter.set_series_data( + series_name=f"F{d}", + x_formatted=sub_data.xval.to_numpy(), + y_formatted=sub_data.yval.to_numpy(), + y_formatted_err=sub_data.yerr.to_numpy(), + ) + + # plot raw RamseyXY plot on second axis + for name in ("Xpos", "Ypos", "Xneg", "Yneg"): + sub_data = curve_data[(curve_data.name == name) & (curve_data.category == "ramsey_xy")] + self.plotter.set_series_data( + series_name=name, + x_formatted=sub_data.xval.to_numpy(), + y_formatted=sub_data.yval.to_numpy(), + y_formatted_err=sub_data.yerr.to_numpy(), + ) + + # find base and amplitude guess + ramsey_xy = curve_data[curve_data.category == "ramsey_xy"] + offset_guess = 0.5 * (ramsey_xy.yval.min() + ramsey_xy.yval.max()) + amp_guess = 0.5 * np.ptp(ramsey_xy.yval) + + # plot frequency and Ramsey fit lines + line_data = curve_data[curve_data.category == "fitted"] for direction in ("pos", "neg"): - ram_x_data = curve_data.get_subset_of(f"X{direction}") - ram_y_data = curve_data.get_subset_of(f"Y{direction}") + sub_data = line_data[line_data.name == f"FREQ{direction}"] + if len(sub_data) == 0: + continue + xval = sub_data.xval.to_numpy() + yn = sub_data.yval.to_numpy() + ys = sub_data.yerr.to_numpy() + yval = unp.uarray(yn, ys) * self._freq_phase_coef() + + # Ramsey fit lines are predicted from the phase fit line. + # Note that this line doesn't need to match with the expeirment data + # because Ramsey P1 data may fluctuate due to phase damping. + + # pylint: disable=no-member + ramsey_cos = amp_guess * unp.cos(yval) + offset_guess + ramsey_sin = amp_guess * unp.sin(yval) + offset_guess + + self.plotter.set_series_data( + series_name=f"F{direction}", + x_interp=xval, + y_interp=yn, + ) + self.plotter.set_series_data( + series_name=f"X{direction}", + x_interp=xval, + y_interp=unp.nominal_values(ramsey_cos), + ) + self.plotter.set_series_data( + series_name=f"Y{direction}", + x_interp=xval, + y_interp=unp.nominal_values(ramsey_sin), + ) - if not np.array_equal(ram_x_data.x, ram_y_data.x): - raise ValueError( - f"The x values for {direction} direction are different in the scan " - "in X and Y quadrature. The scan must be identical in both quadrature " - "to compute initial guess." + if np.isfinite(ys).all(): + self.plotter.set_series_data( + series_name=f"F{direction}", + y_interp_err=ys, ) - xvals = ram_x_data.x - - # Get normalized sinusoidals - xnorm = (ram_x_data.y - est_offs) / est_a - ynorm = (ram_y_data.y - est_offs) / est_a - - # Compute derivative to extract polynomials from sinusoidal - dx = np.diff(xnorm) / np.diff(xvals) - dy = np.diff(ynorm) / np.diff(xvals) - - # Eliminate sinusoidal - phase_poly = np.sqrt(dx**2 + dy**2) - - # Do polyfit up to 2rd order. - # This must correspond to the 3rd order in the original function. - vmat_xpoly = np.vstack((xvals[1:] ** 2, xvals[1:], np.ones(xvals.size - 1))).T - coeffs = np.linalg.lstsq(vmat_xpoly, phase_poly, rcond=-1)[0] - poly_guess = { - f"c1_{direction}": coeffs[2] / 1 / const, - f"c2_{direction}": coeffs[1] / 2 / const, - f"c3_{direction}": coeffs[0] / 3 / const, - } - user_opt.p0.set_if_empty(**poly_guess) - - return user_opt + self.plotter.set_series_data( + series_name=f"X{direction}", + y_interp_err=unp.std_devs(ramsey_cos), + ) + self.plotter.set_series_data( + series_name=f"Y{direction}", + y_interp_err=unp.std_devs(ramsey_sin), + ) + return [self.plotter.figure()] def _initialize( self, @@ -503,6 +602,5 @@ def _initialize( super()._initialize(experiment_data) # Set scaling factor to convert phase to frequency - fixed_params = self.options.fixed_parameters.copy() - fixed_params["ts"] = experiment_data.metadata["stark_length"] - self.set_options(fixed_parameters=fixed_params) + if "stark_length" in experiment_data.metadata: + self.set_options(pulse_len=experiment_data.metadata["stark_length"]) diff --git a/test/library/characterization/test_stark_ramsey_xy.py b/test/library/characterization/test_stark_ramsey_xy.py index 30ac1bdb2f..795fde5297 100644 --- a/test/library/characterization/test_stark_ramsey_xy.py +++ b/test/library/characterization/test_stark_ramsey_xy.py @@ -224,15 +224,16 @@ def test_circuit_valid_delays(self): self.assertEqual(stark_u.duration % backend.target.granularity, 0) @named_data( - ["ideal_quadratic", 0.5, 0.0, 30e6, 0.0, 0.0, -30e6, 0.0, 0.0, 0.5], - ["with_all_terms", 0.5, 15e6, 200e6, -100e6, 15e6, -200e6, -100e6, 300e3, 0.5], - ["with_spam_like", 0.3, 0.0, 30e6, 0.0, 0.0, -30e6, 0.0, 0.0, 0.4], - ["asymmetric_shift", 0.4, -20e6, 200e6, -100e6, -15e6, -180e6, -90e6, 200e3, 0.5], - ["large_cubic_term", 0.5, 10e6, 15e6, 30e6, 5e6, -10e6, 40e6, 0.0, 0.5], + ["ideal_quadratic", 0.0, 30e6, 0.0, 0.0, -30e6, 0.0, 0.0], + ["with_all_terms", 15e6, 200e6, -100e6, 15e6, -200e6, -100e6, 300e3], + ["asymmetric_shift", -20e6, 200e6, -100e6, -15e6, -180e6, -90e6, 200e3], + ["large_cubic_term", 10e6, 15e6, 30e6, 5e6, -10e6, 40e6, 0.0], ) @unpack - def test_ramsey_fast_analysis(self, amp, c1p, c2p, c3p, c1n, c2n, c3n, ferr, off): + def test_ramsey_fast_analysis(self, c1p, c2p, c3p, c1n, c2n, c3n, ferr): """End-to-end test for Ramsey fast analysis with artificial data.""" + amp = 0.5 + off = 0.5 rng = np.random.default_rng(seed=123) shots = 1000 @@ -270,17 +271,38 @@ def test_ramsey_fast_analysis(self, amp, c1p, c2p, c3p, c1n, c2n, c3n, ferr, off analysis.run(exp_data, replace_results=True) self.assertExperimentDone(exp_data) - # Check if fit parameters are consistent with input values - to_test = { - "stark_pos_coef_o1": c1p, - "stark_pos_coef_o2": c2p, - "stark_pos_coef_o3": c3p, - "stark_neg_coef_o1": c1n, - "stark_neg_coef_o2": c2n, - "stark_neg_coef_o3": c3n, - "stark_ferr": ferr, - } - for name, refv in to_test.items(): - # Error must be within 1 percent or 1 MHz - val = exp_data.analysis_results(name).value.n - self.assertAlmostEqual(val, refv, delta=3e6) + # Check the fitted parameter can approximate the same polynominal + x_pos = np.linspace(0, 1, 51) + x_neg = np.linspace(-1, 0, 51) + ref_yvals_pos = c1p * x_pos + c2p * x_pos**2 + c3p * x_pos**3 + ferr + ref_yvals_neg = c1n * x_neg + c2n * x_neg**2 + c3n * x_neg**3 + ferr + + # Note that these parameter values are not necessary the same with input values + # as long as they can approximate the original phase polynominal. + c1p_est = exp_data.analysis_results("stark_pos_coef_o1").value.n + c2p_est = exp_data.analysis_results("stark_pos_coef_o2").value.n + c3p_est = exp_data.analysis_results("stark_pos_coef_o3").value.n + c1n_est = exp_data.analysis_results("stark_neg_coef_o1").value.n + c2n_est = exp_data.analysis_results("stark_neg_coef_o2").value.n + c3n_est = exp_data.analysis_results("stark_neg_coef_o3").value.n + ferr_est = exp_data.analysis_results("stark_ferr").value.n + + test_yvals_pos = c1p_est * x_pos + c2p_est * x_pos**2 + c3p_est * x_pos**3 + ferr_est + test_yvals_neg = c1n_est * x_neg + c2n_est * x_neg**2 + c3n_est * x_neg**3 + ferr_est + + # Check similality of reconstructed polynominals + # Curves must be agree within the torelance of 1.5 * 1 MHz. + np.testing.assert_array_almost_equal( + test_yvals_pos, + ref_yvals_pos, + decimal=-6, + err_msg="Reconstructed phase polynominal on positive frequency shift side " + "doesn't match with the original curve.", + ) + np.testing.assert_array_almost_equal( + test_yvals_neg, + ref_yvals_neg, + decimal=-6, + err_msg="Reconstructed phase polynominal on negative frequency shift side " + "doesn't match with the original curve.", + )