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.", + )