Tests are working with these base class changes

pycqed/analysis/fitting_models.py

@@ -334,9 +334,8 @@ def HangerFuncComplex(f, pars):
 
     return S21
 
-def hanger_func_complex_SI(f: float, f0: float, Ql: float, Qe: float,
-                           A: float, theta: float, phi_v: float, phi_0: float,
-                           alpha:float =1):
+def hanger_func_complex_SI(f, f0, Q, Qe,A, theta, phi_v, phi_0,slope =1):
+
     '''
     This is the complex function for a hanger (lamda/4 resonator).
     See equation 3.1 of the Asaad master thesis.
@@ -349,23 +348,37 @@ def hanger_func_complex_SI(f: float, f0: float, Ql: float, Qe: float,
         f   : frequency
         f0  : resonance frequency
         A   : background transmission amplitude
-        Ql  : loaded quality factor
+        Q  : loaded quality factor
         Qe  : extrinsic quality factor
         theta:  phase of Qe (in rad)
         phi_v:  phase to account for propagation delay to sample
         phi_0:  phase to account for propagation delay from sample
-        alpha:  slope of signal around the resonance
+        slope:  slope of signal around the resonance
 
     '''
-    slope_corr = (1+alpha*(f-f0)/f0)
+    slope_corr = (1+slope*(f-f0)/f0)
     propagation_delay_corr = np.exp(1j * (phi_v * f + phi_0))
-    hanger_contribution = (1 - Ql / Qe * np.exp(1j * theta)/
-                               (1 + 2.j * Ql * (f  - f0) / f0))
+    hanger_contribution = (1 - Q / Qe * np.exp(1j * theta)/
+                               (1 + 2.j * Q * (f  - f0) / f0))
     S21 = A *  slope_corr * hanger_contribution * propagation_delay_corr
 
     return S21
 
-
+def hanger_func_complex_SI_pars(f,pars):
+     '''
+     This function is used in the minimization fitting which requires parameters
+     '''
+
+     f0 = pars['f0']
+     Ql = pars['Ql']
+     Qe = pars['Qe']
+     A = pars['A']
+     theta = pars['theta']
+     phi_v = pars['phi_v']
+     phi_0 = pars['phi_0']
+     alpha = pars['alpha']
+     return hanger_func_complex_SI(f, f0, Ql, Qe,
+                            A, theta, phi_v, phi_0, alpha)
 
 def PolyBgHangerFuncAmplitude(f, f0, Q, Qe, A, theta, poly_coeffs):
     # This is the function for a hanger (lambda/4 resonator) which takes into
@@ -991,6 +1004,110 @@ def ro_double_gauss_guess(model, data, x, fixed_p01 = False, fixed_p10 = False):
     return model.make_params()
 
 
+def SlopedHangerFuncAmplitudeGuess(data, f, fit_window=None):
+     xvals = f
+     peaks = a_tools.peak_finder(xvals, data)
+       # Search for peak
+     if peaks['dip'] is not None:    # look for dips first
+         f0 = peaks['dip']
+         amplitude_factor = -1.
+     elif peaks['peak'] is not None:  # then look for peaks
+         f0 = peaks['peak']
+         amplitude_factor = 1.
+     else:                                 # Otherwise take center of range
+         f0 = np.median(xvals)
+         amplitude_factor = -1.
+
+     min_index = np.argmin(data)
+     max_index = np.argmax(data)
+     min_frequency = xvals[min_index]
+     max_frequency = xvals[max_index]
+
+     amplitude_guess = max(dm_tools.reject_outliers(data))
+
+     # Creating parameters and estimations
+     S21min = (min(dm_tools.reject_outliers(data)) /
+               max(dm_tools.reject_outliers(data)))
+
+     Q = f0 / abs(min_frequency - max_frequency)
+
+     Qe = abs(Q / abs(1 - S21min))
+     guess_dict = {'f0': {'value': f0*1e-9,
+                          'min': min(xvals)*1e-9,
+                          'max': max(xvals)*1e-9},
+                   'A': {'value': amplitude_guess},
+                   'Q': {'value': Q,
+                         'min': 1,
+                         'max': 50e6},
+                   'Qe': {'value': Qe,
+                          'min': 1,
+                          'max': 50e6},
+                   'Qi': {'expr': 'abs(1./(1./Q-1./Qe*cos(theta)))',
+                          'vary': False},
+                   'Qc': {'expr': 'Qe/cos(theta)',
+                          'vary': False},
+                   'theta': {'value': 0,
+                             'min': -np.pi/2,
+                             'max': np.pi/2},
+                   'slope': {'value':0,
+                             'vary':True}}
+     return guess_dict
+
+
+def hanger_func_complex_SI_Guess(data, f, fit_window=None):
+     ## This is complete garbage, just to get some return value
+     xvals = f
+     abs_data = np.abs(data)
+     peaks = a_tools.peak_finder(xvals, abs_data)
+       # Search for peak
+     if peaks['dip'] is not None:    # look for dips first
+         f0 = peaks['dip']
+         amplitude_factor = -1.
+     elif peaks['peak'] is not None:  # then look for peaks
+         f0 = peaks['peak']
+         amplitude_factor = 1.
+     else:                                 # Otherwise take center of range
+         f0 = np.median(xvals)
+         amplitude_factor = -1.
+
+     min_index = np.argmin(abs_data)
+     max_index = np.argmax(abs_data)
+     min_frequency = xvals[min_index]
+     max_frequency = xvals[max_index]
+
+     amplitude_guess = max(dm_tools.reject_outliers(abs_data))
+
+     # Creating parameters and estimations
+     S21min = (min(dm_tools.reject_outliers(abs_data)) /
+               max(dm_tools.reject_outliers(abs_data)))
+
+     Q = f0 / abs(min_frequency - max_frequency)
+
+     Qe = abs(Q / abs(1 - S21min))
+     guess_dict = {'f0': {'value': f0*1e-9,
+                          'min': min(xvals)*1e-9,
+                          'max': max(xvals)*1e-9},
+                   'A': {'value': amplitude_guess},
+                   'Qe': {'value': Qe,
+                          'min': 1,
+                          'max': 50e6},
+                   'Ql': {'value': Q,
+                          'min': 1,
+                          'max': 50e6},
+                   'theta': {'value': 0,
+                             'min': -np.pi/2,
+                             'max': np.pi/2},
+                   'alpha': {'value':0,
+                             'vary':True},
+                   'phi_0': {'value':0,
+                             'vary':True},
+                   'phi_v': {'value':0,
+                             'vary':True}}
+
+     return guess_dict
+
+
+
 def sum_int(x,y):
     return np.cumsum(y[:-1]*(x[1:]-x[:-1]))
 

pycqed/analysis_v2/base_analysis.py

@@ -224,7 +224,7 @@ def run_analysis(self):
             self.plot(key_list='auto')  # make the plots
 
         if self.options_dict.get('save_figs', False):
-            self.save_figures(close_figs=self.options_dict.get('close_figs', False))
+            self.save_figures(close_figs=self.options_dict.get('close_figs', True))
 
     def get_timestamps(self):
         """
@@ -485,40 +485,70 @@ def run_fitting(self):
             fit_yvals = fit_dict['fit_yvals']
             fit_xvals = fit_dict['fit_xvals']
 
+            fitting_type = fit_dict.get('fitting_type', 'model')
+
             model = fit_dict.get('model', None)
             if model is None:
                 fit_fn = fit_dict.get('fit_fn', None)
                 model = fit_dict.get('model', lmfit.Model(fit_fn))
             fit_guess_fn = fit_dict.get('fit_guess_fn', None)
-            if fit_guess_fn is None and fit_dict.get('fit_guess', True):
-                fit_guess_fn = model.guess
+            if fit_guess_fn is None:
+                if  fitting_type == 'model' and fit_dict.get('fit_guess', True):
+                    fit_guess_fn = model.guess
 
             if guess_pars is None:
                 if fit_guess_fn is not None:
+                    if fitting_type is 'minimize':
+                        guess_pars = fit_guess_fn(**fit_yvals, **fit_xvals, **guessfn_pars)
+                        params = lmfit.Parameters()
+                        for gd_key, val in guess_pars.items():
+                            params.add(gd_key)
+                            for attr, attr_val in val.items():
+                                setattr(params[gd_key], attr, attr_val)
                     # a fit function should return lmfit parameter objects
                     # but can also work by returning a dictionary of guesses
-                    guess_pars = fit_guess_fn(**fit_yvals, **fit_xvals, **guessfn_pars)
-                    if not isinstance(guess_pars, lmfit.Parameters):
-                        for gd_key, val in list(guess_pars.items()):
-                            model.set_param_hint(gd_key, **val)
-                        guess_pars = model.make_params()
-
-                    if guess_dict is not None:
-                        for gd_key, val in guess_dict.items():
-                            for attr, attr_val in val.items():
-                                # e.g. setattr(guess_pars['frequency'], 'value', 20e6)
-                                setattr(guess_pars[gd_key], attr, attr_val)
-                    # A guess can also be specified as a dictionary.
-                    # additionally this can be used to overwrite values
-                    # from the guess functions.
+                    elif fitting_type is 'model':
+                        guess_pars = fit_guess_fn(**fit_yvals, **fit_xvals, **guessfn_pars)
+                        if not isinstance(guess_pars, lmfit.Parameters):
+                            for gd_key, val in list(guess_pars.items()):
+                                model.set_param_hint(gd_key, **val)
+                            guess_pars = model.make_params()
+                                                ####
+                        if guess_dict is not None:
+                             for gd_key, val in guess_dict.items():
+                                 for attr, attr_val in val.items():
+                                     # e.g. setattr(guess_pars['frequency'], 'value', 20e6)
+                                     setattr(guess_pars[gd_key], attr, attr_val)
+                        ####
                 elif guess_dict is not None:
-                    for key, val in list(guess_dict.items()):
-                        model.set_param_hint(key, **val)
-                    guess_pars = model.make_params()
-            fit_dict['fit_res'] = model.fit(**fit_xvals, **fit_yvals,
-                                            params=guess_pars)
+                    if fitting_type is 'minimize':
+                        params = lmfit.Parameters()
+                        for key, val in list(guess_dict.items()):
+                         # for key, val in guess_dict.items():
+                            params.add(key)
+                            for attr, attr_val in val.items():
+                                setattr(params[key], attr, attr_val)
+                    elif fitting_type is 'model':
+                        for key, val in list(guess_dict.items()):
+                            model.set_param_hint(key, **val)
+                        guess_pars = model.make_params()
+            else:
+                if fitting_type is 'minimize':
+                    raise NotImplementedError('Conversion from guess_pars to params with lmfit.Parameters() needs to be implemented')
+            if fitting_type is 'model':
+                fit_dict['fit_res'] = model.fit(**fit_xvals, **fit_yvals,
+                                                params=guess_pars)
+                self.fit_res[key] = fit_dict['fit_res']
+
+            elif fitting_type is 'minimize':
+                fit_dict['fit_res'] = lmfit.minimize(fcn=_complex_residual_function,
+                        params=params,
+                        args=(fit_fn, fit_xvals, fit_yvals))
+                fit_dict['fit_res'].initial_params = params
+                fit_dict['fit_res'].userkws = fit_xvals
+                fit_dict['fit_res'].fit_fn = fit_fn
+                self.fit_res[key] = fit_dict['fit_res']
 
-            self.fit_res[key] = fit_dict['fit_res']
 
     def save_fit_results(self):
         """
@@ -1173,29 +1203,51 @@ def plot_fit(self, pdict, axs):
             logging.warning('fit_res is an empty dictionary, cannot plot.')
             return
 
-        model = pdict['fit_res'].model
-        #TODO: make this work for type(model) is lmfit.minimizer.MinimizerResult
-        assert (isinstance(model, lmfit.model.Model) or
-                    isinstance(model, lmfit.model.ModelResult),
-            'The passed item in "fit_res" needs to be a fitting model, but is {}'.format(type(model)))
-
         plot_init = pdict.get('plot_init', False)  # plot the initial guess
         plot_normed = pdict.get('plot_normed', False)
         pdict['marker'] = pdict.get('marker', '')  # different default
         plot_linestyle_init = pdict.get('init_linestyle', '--')
         plot_numpoints = pdict.get('num_points', 1000)
 
-        if len(model.independent_vars) == 1:
-            independent_var = model.independent_vars[0]
+        if hasattr(pdict['fit_res'],'model'):
+            model = pdict['fit_res'].model
+            assert (isinstance(model, lmfit.model.Model) or
+                        isinstance(model, lmfit.model.ModelResult),
+                'The passed item in "fit_res" needs to be a fitting model, but is {}'.format(type(model)))
+
+            if len(model.independent_vars) == 1:
+                independent_var = model.independent_vars[0]
+            else:
+                raise ValueError('Fit can only be plotted if the model function'
+                                 ' has one independent variable.')
+            x_arr = pdict['fit_res'].userkws[independent_var]
+            pdict['xvals'] = np.linspace(np.min(x_arr), np.max(x_arr),
+                                         plot_numpoints)
+            pdict['yvals'] = model.eval(pdict['fit_res'].params,
+                                        **{independent_var: pdict['xvals']})
         else:
-            raise ValueError('Fit can only be plotted if the model function'
-                             ' has one independent variable.')
-
-        x_arr = pdict['fit_res'].userkws[independent_var]
-        pdict['xvals'] = np.linspace(np.min(x_arr), np.max(x_arr),
-                                     plot_numpoints)
-        pdict['yvals'] = model.eval(pdict['fit_res'].params,
-                                    **{independent_var: pdict['xvals']})
+
+            '''
+            This is the case for the minimizier fit
+            '''
+            fit_xvals = pdict['fit_res'].userkws
+            if len(fit_xvals.keys()) == 1:
+                independent_var = list(fit_xvals.keys())[0]
+            else:
+                raise ValueError('Fit can only be plotted if the model function'
+                                  ' has one independent variable.')
+            x_arr = pdict['fit_res'].userkws[independent_var]
+            pdict['xvals'] = np.linspace(np.min(x_arr), np.max(x_arr),
+                                         plot_numpoints)
+            fit_fn = pdict['fit_res'].fit_fn
+            output = fit_fn(**pdict['fit_res'].params,
+                            **{independent_var: pdict['xvals']})
+            output_mod_fn = pdict.get('output_mod_fn', None)
+            if output_mod_fn is not None:
+                pdict['yvals'] = output_mod_fn(output)
+            else:
+                pdict['yvals'] = output
+
         if plot_normed:
             pdict['yvals']=pdict['yvals']/pdict['yvals'][0]
 
@@ -1205,9 +1257,20 @@ def plot_fit(self, pdict, axs):
             # The initial guess
             pdict_init = copy.copy(pdict)
             pdict_init['linestyle'] = plot_linestyle_init
-            pdict_init['yvals'] = model.eval(
-                **pdict['fit_res'].init_values,
-                **{independent_var: pdict['xvals']})
+            if hasattr(pdict_init['fit_res'],'model'):
+                # The initial guess
+                pdict_init['yvals'] = model.eval(
+                    **pdict_init['fit_res'].init_values,
+                    **{independent_var: pdict_init['xvals']})
+            else:
+                output = fit_fn(**pdict_init['fit_res'].initial_params,
+                            **{independent_var: pdict_init['xvals']})
+                output_mod_fn = pdict_init.get('output_mod_fn', None)
+                if output_mod_fn is not None:
+                    pdict_init['yvals'] = output_mod_fn(output)
+                else:
+                    pdict_init['yvals'] = output
+
             pdict_init['setlabel'] += ' init'
             self.plot_line(pdict_init, axs)
 
@@ -1402,3 +1465,18 @@ def _merge_dict_rec(dict_a: dict, dict_b: dict):
         if k not in dict_a:
             dict_a[k] = dict_b[k]
     return dict_a
+
+
+def _complex_residual_function(pars, fit_function, x, data):
+    '''
+    Residual of a complex function with complex results dictionary data and
+    real input values dictionary 'x'
+    data should be of the format data = {'data':[np.array]}
+    x should be of the format x = {'x':[np.array}} where x is the variable
+    name in the function fit_function
+    For resonators 'x' is the the frequency, 'data' the complex transmission
+    '''
+    cmp_values = fit_function(**x, **pars)
+    res = cmp_values-data['data']
+    res = np.append(res.real, res.imag)
+    return res