ResonatorModel giving poor guesses

Issue

I have been noticing that the fit provided by ResonatorSpectroscopyAnalysis is not very robust. In fact, I faced it many times in the lab, even with clean and well-behaved datasets that should have resulted in good fits.

The root of the issue appears to lie (at least partly) in the way the phase guesses are generated.

Code example with a failed fit

import numpy as np
import matplotlib.pyplot as plt
np.set_printoptions(precision=3)
from quantify_core.analysis import fitting_models as fm

sorted_frequencies= np.array([75000000.  , 75625000.  , 76250000.  , 76875000.  , 77500000.  ,
       78125000.  , 78281250.  , 78359375.  , 78437500.  , 78476562.5 ,
       78515625.  , 78535156.25, 78554687.5 , 78574218.75, 78593750.  ,
       78632812.5 , 78652343.75, 78671875.  , 78710937.5 , 78750000.  ,
       78828125.  , 78906250.  , 79062500.  , 79218750.  , 79296875.  ,
       79375000.  , 79453125.  , 79531250.  , 79687500.  , 79843750.  ,
       80000000.  , 80312500.  , 80625000.  , 81250000.  , 81875000.  ,
       82500000.  , 83125000.  , 83750000.  , 84375000.  , 85000000.  ])
sorted_IQ= np.array([ 0.00519+0.01065j,  0.01118+0.00364j,  0.01015-0.0048j ,
        0.00444-0.0111j , -0.0053 -0.01126j, -0.01218-0.00298j,
       -0.01221+0.00051j, -0.0116 +0.00252j, -0.01046+0.00397j,
       -0.00938+0.00468j, -0.00833+0.0053j , -0.00763+0.00525j,
       -0.00703+0.00518j, -0.00667+0.00487j, -0.00581+0.00491j,
       -0.00541+0.00438j, -0.00468+0.00402j, -0.00426+0.00399j,
       -0.00376+0.00334j, -0.0033 +0.003j  , -0.00293+0.00217j,
       -0.00262+0.00153j, -0.00297+0.00054j, -0.00411+0.00053j,
       -0.00469+0.00101j, -0.00536+0.0015j , -0.00571+0.00235j,
       -0.00604+0.00323j, -0.00548+0.00521j, -0.0046 +0.00685j,
       -0.00308+0.00832j,  0.00105+0.0093j ,  0.00503+0.00833j,
        0.01011+0.00196j,  0.0089 -0.00607j,  0.00222-0.01101j,
       -0.00596-0.00991j, -0.01104-0.00348j, -0.01026+0.00503j,
       -0.00403+0.01062j])


model = fm.ResonatorModel()
guess = model.guess(sorted_IQ, f=sorted_frequencies,fit_kws={'xtol':1e-15, 'ftol':1e-15})
#******** Working Delay Corrctions ***********
# guess['fr'].value = 7.90625e+07
# guess['Ql'].value = 59.01213
# guess['Qe'].value = 57.84357
# guess['A'].value  = 0.00865632
# guess['theta'].value = 0
# guess['phi_0'].value = 66.55737
# guess['phi_v'].value = -8.725337e-07

#########

fit_result = model.fit(sorted_IQ, params=guess, f=sorted_frequencies,
		fit_kws={'xtol':1e-15, 'ftol':1e-15})

n_fit_frequencies = 2000
independent_var = model.independent_vars[0]
freq_array = fit_result.userkws[independent_var]
fit_freqs = np.linspace(np.min(freq_array), np.max(freq_array), n_fit_frequencies)
print( fit_result.fit_report() )

fit_y = model.eval(fit_result.params, **{independent_var: fit_freqs})
range_casting = 'abs'
range_cast_func = getattr(np, range_casting)
fit_y = range_cast_func(fit_y)

plt.plot( sorted_frequencies, np.abs(sorted_IQ),'bo',ms=3.0)
plt.plot( fit_freqs, fit_y,'r-',ms=3.0)
plt.show()

Which yields

image

Test fix

I tried generating the guesses for phi_0 and phi_v the same way as in PycQED. Namely, I looked at the following code: https://github.com/DiCarloLab-Delft/PycQED_py3/blob/develop/pycqed/analysis/measurement_analysis.py#L6225-L6235

I tried replicating it in the resonator_phase_guess function in quantify_core.analysis.fitting_models. The new code is:

def resonator_phase_guess(s21: np.ndarray, freq: np.ndarray) -> Tuple[float, float]:
    """
    Guesses the phase velocity in resonator spectroscopy,
    based on the median of all the differences between consecutive phases.

    Parameters
    ----------
    s21:
        Resonator S21 data
    freq:
        Frequency of the spectroscopy pulse

    Returns
    -------
    phi_0:
        Guess for the phase offset
    phi_v:
        Guess for the phase velocity
    """
    phase = np.angle(s21)
    argmin_s21 = np.abs(s21).argmin()
    fr_guess = freq[argmin_s21]

    nbr_phase_jumps = (np.diff(phase) > 4).sum()
    guess_phi_v = (2 * np.pi * nbr_phase_jumps + (phase[0] - phase[-1])) / (
    min(freq) - max(freq))
    # phi_0
    angle_resonance = phase[int(len(s21) / 2)]
    phase_evolution_resonance = np.exp(1j * guess_phi_v * fr_guess)
    angle_phase_evolution = np.arctan2(
    np.imag(phase_evolution_resonance), np.real(phase_evolution_resonance))
    guess_phi_0 = angle_resonance - angle_phase_evolution

    return guess_phi_0, guess_phi_v

The result:

image