Python lmfit.Model类代码示例

def gauss_step_const(signal, guess):
    """
    Fits high contrast data very well
    """
    if guess == False:
        return [0, 0]
    else:
        amp, centre, stdev, offset = guess
        
        data = np.array([range(len(signal)), signal]).T
        X = data[:,0]
        Y = data[:,1]

#         gauss_mod = Model(gaussian)
        gauss_mod = Model(gaussian)
        const_mod = ConstantModel()
        step_mod = StepModel(prefix='step')
        
        pars = gauss_mod.make_params(height=amp, center=centre, width=stdev / 3., offset=offset)
#         pars = gauss_mod.make_params(amplitude=amp, center=centre, sigma=stdev / 3.)
        gauss_mod.set_param_hint('sigma', value = stdev / 3., min=stdev / 2., max=stdev)
        pars += step_mod.guess(Y, x=X, center=centre)

        pars += const_mod.guess(Y, x=X)
    
        
        mod = const_mod + gauss_mod + step_mod
        result = mod.fit(Y, pars, x=X)
        # write error report
        #print result.fit_report()
        print "contrast fit", result.redchi
    
    return X, result.best_fit, result.redchi

def lmDDOPhaseFit(xdata, ydata, params, f0_range = 1.2, Q_range = 3):    
    f0= params[0]
    Q=1/(2*params[2])
    
    x = xdata
    y = ydata
#Define a linear model and a Damped Oscillator Model    
#    ddophase_mod = ExpressionModel('off + m*x- arctan(1/Q * 1/(f0/x - x/f0))-')
    ddophase_mod = Model(phase)
#Initial Pars for Linear Model
    pars =  ddophase_mod.make_params(off=0, m=0, f0=f0, Q=Q)

#Add fit parameters, Center, Amplitude, and Sigma
    pars['f0'].set(min=f0/f0_range, max=f0*f0_range)
    pars['Q'].set(min=Q/Q_range, max=Q*Q_range)
#Create full model. Add linear model and all peaks
#Initialize fit
    init = ddophase_mod.eval(pars, x=x)
#Do the fit. The weight exponential can weight the points porportional to the
#amplitude of y point. In this way, points on peak can be given more weight.     
    out=ddophase_mod.fit(y, pars,x=x)
#Get the fit parameters
    fittedf0= out.params['f0'].value
    fittedQ = out.params['Q'].value
#Returns the output fit as well as an array of the fit parameters
    """Returns output fit as will as list of important fitting parameters"""
    return out, [fittedf0, np.nan, np.nan, fittedQ]

def rescale_reframe(scisub, refsub, verbose=False):
    """Example function with types documented in the docstring.

        `PEP 484`_ type annotations are supported. If attribute, parameter, and
        return types are annotated according to `PEP 484`_, they do not need to be
        included in the docstring:

        Args:
            param1 (int): The first parameter.
            param2 (str): The second parameter.

        Returns:
            bool: The return value. True for success, False otherwise.

        .. _PEP 484:
            https://www.python.org/dev/peps/pep-0484/

    """
    
    params = Parameters()
    params.add('sigma', 1.0, True, 0.0, inf)
    
    image_sub         = Model(res_sigma_image_flat, independent_vars=['scisub', 'refsub'])
    image_sub_results = image_sub.fit(data=scisub.ravel(), params=params, scisub=scisub, refsub=refsub)
    
    if verbose: print('Sigma of Rescale: {}'.format(image_sub_results.params['sigma'].value))
    
    return partial_res_sigma_image(image_sub_results.params['sigma'].value)

def make_models(model,peaks):
    """ Make composite models for multiple peaks 
        
        Arguments:
            -- models
            -- peaks: dict of {prefix:"peak_name_", model_params:()}
    """
    if len(peaks)<2:
        mod = Model(model,prefix=peaks[0][0])
        p_guess = mod.make_params()
        values = peaks[0][1:]
        update_params(p_guess,values)

    elif len(peaks)>1:

        mod = Model(model,prefix=peaks[0][0])
        
        values = peaks[0][1:]
        print(values)
        for i in peaks[1:]:
            mod += Model(model,prefix=i[0])
            values.extend(i[1:])

        p_guess = mod.make_params()
        
        update_params(p_guess,values)

        return mod, p_guess

def get_resvar(x, y, mod='linear', zero_intercept=False):
    """
    If `mod == 'linear'` then uses linear regression; otherwise
    `mod == 'sigmoid'` and uses sigmoid regression.
    """
    assert len(x) == len(y)
    if len(x) < 3 or mod == 'linear':
        model = Model(linear)
        params = model.make_params()
        params['a'].value = (y.max() - y.min()) / (x.max() - x.min())
        params['b'].value = y.min()
        if zero_intercept:
            params['b'].value = 0
            params['b'].vary = False
        modfit = model.fit(y, x=x, params=params)
    else:
        model = Model(sigmoid)
        modfit = model.fit(y, x=x, a=y.max(), b=1 / (x.max() - x.min()),
                           c=x.max())
    finex = np.linspace(x.min(), x.max(), 50)
    result = {
        'mod': mod,
        'params': modfit.best_values,
        'resvar': modfit.residual.var(),
        'y': modfit.best_fit,
        'finex': finex,
        'finey': modfit.eval(x=finex),
        'r2': 1 - modfit.residual.var() / np.var(y),
        'modfit': modfit}
    return result

def fit_data_bg(x, y, peak_pos, peak_type='LO', width=None, bg_ord=0):
    """ 
    Builds a lmfit model of peaks in listed by index in `peak_pos`

    Parameters
    ----------
    peak_type : string (default='lorentizian')
        Peaks can be of the following types:

        - 'LO' : symmetric lorentzian
        - 'GA' : symmetric gaussain
        - 'VO' : symmetric pseudo voigt

    max_width : int (default = total points/10)
        max width (in data points) that peak fitted can be

    bg_ord: int
        order of the background polynomial
        0: constant, 1: linear, ...

    Returns
    -------
    out: 
        fitted model

    """
    # need to define peak width finding
    if width is None:
        width = guess_peak_width(x, y)
    
    # start with polynomial background
    model = PolynomialModel(bg_ord, prefix='bg_')
    pars = model.make_params()

    if peak_type == 'LO':
        peak_function = lorentzian
    elif peak_type == 'GA':
        peak_function = gaussian
    elif peak_type == 'VO':
        peak_function = voigt

    # add peak type for all peaks
    for i, peak in enumerate(peak_pos):
        temp_model = Model(peak_function, prefix='p%s_' % i)
        pars.update(temp_model.make_params())
        model += temp_model

    # set initial background as flat line at zeros
    for i in range(bg_ord + 1):
        pars['bg_c%i' % i].set(0)

    # give values for other peaks, keeping width and height positive
    for i, peak in enumerate(peak_pos):
        pars['p%s_x0' % i].set(x[peak])
        pars['p%s_fwhm' % i].set(width, min=0)
        pars['p%s_amp' % i].set(y[peak], min=0)

    out = model.fit(y, pars, x=x)
    return out

def multi_disk_fit(arcmin_sma, intens, intens_err, bounds, psf):
    h_mat_0 = numpy.ones(len(bounds - 1))
    mod = Model(multi_disk_func)
    pars = mod.make_params(I0 = 1.0, h_mat = h_mat_0, bounds = bounds, PSF = PSF)
    pars['bounds'].vary = False
    pars['PSF'].vary = False
    results = mod.fit(intens, params = pars, x = arcmin_sma, weights = 1.0/intens_err)
    pdb.set_trace()

def fit_D_phi_powerlaw(phi_, D):
    """fit stuff"""
    Dmod = Model(Dhop_powerlaw)
    params = Dmod.make_params() 
    params['p'].set(0.3)
    result = Dmod.fit(D, params, phi=phi_)
    print result.fit_report(min_correl=0.5)
    p = result.params['p'].value
    return p

def lmfit(mjd,flux,fluxerr):
    t0_guess = mjd[np.argmax(flux)]
    tau_fall_guess = 40.
    tau_rise_guess = -5.
    A = 150.
    B = 20.
    # nflux = np.zeros(2+len(np.array(flux)))
    # nfluxerr = np.ones(2+len(np.array(flux)))/10.
    # nmjd = np.zeros(2+len(np.array(flux)))
    #
    # nflux[1:-1] = flux
    # nfluxerr[1:-1] = fluxerr
    # nmjd[1:-1] = mjd
    # nmjd[1] = mjd[0]-100.
    # nmjd[-1] = mjd[-1]+150
    #
    # flux = nflux
    # fluxerr = nfluxerr
    # mjd = nmjd

    bmod = Model(bazinfunc)
    bmod.set_param_hint('t0', value=t0_guess, min=t0_guess-20, max=t0_guess+20)
    bmod.set_param_hint('tau_fall', value=tau_fall_guess)
    bmod.set_param_hint('tau_rise', value=tau_rise_guess)
    bmod.set_param_hint('A',value=A)
    bmod.set_param_hint('B',value=B)

    pars = bmod.make_params()
    #print(bmod.param_names)
    #print(bmod.independent_vars)
    # print(np.array(flux))
    # print(np.array(1./np.array(fluxerr)))
    # print(np.array(mjd))
    result = bmod.fit(np.array(flux),method='leastsq',weights=1./np.array(fluxerr), t=np.array(mjd))

    #print(result.fit_report())
    # plt.clf()
    # plt.errorbar(np.array(mjd), np.array(flux), yerr=fluxerr,fmt='o')
    # plt.plot(np.array(mjd), result.init_fit, 'k--')
    # plt.plot(np.array(mjd), result.best_fit, 'r-')
    # #plt.xlim(mjd[1],mjd[-2])
    # plt.savefig('bazinfit.png')



    chisq = result.redchi
    ps = result.best_values
    popt = [ps['t0'],ps['tau_fall'],ps['tau_rise'],ps['A'],ps['B']]
    #print('popt',popt)
    #sys.exit()
    # if chisq < 2.:
    #     input('good chisq!')

    # popt, pcov, infodict, errmsg, ier = curve_fit(bazinfunc, mjd, flux,
    #                                               sigma=fluxerr, p0=p0, maxfev=2000000, full_output=True)
    #
    # chisq = (infodict['fvec'] ** 2).sum() / (len(infodict['fvec']) - len(popt))
    return chisq,popt

    def test_extra_param_issues_warning(self):
        # The function accepts extra params, Model will warn but not raise.
        def flexible_func(x, amplitude, center, sigma, **kwargs):
            return gaussian(x, amplitude, center, sigma)

        flexible_model = Model(flexible_func)
        pars = flexible_model.make_params(**self.guess())
        with warnings.catch_warnings(record=True) as w:
            warnings.simplefilter("always")
            flexible_model.fit(self.data, pars, x=self.x, extra=5)
        self.assertTrue(len(w) == 1)
        self.assertTrue(issubclass(w[-1].category, UserWarning))

def linefit(fitam, fitcps, fitweights):
    try:
        model = Model(linefunc)
        params = Parameters()
        params.add('linem', value=5)
        params.add('linek', value=0.25, min=0, max=1)
        fit = model.fit(asarray(fitcps), params, fitweights, linex=asarray(fitam))
        fitvalues = fit.best_values
        fiterr = fit.covar
    except RuntimeError:
        return float('nan'), float('nan')
    return [[fitvalues['linem'], fitvalues['linek']], fiterr]

def Gaussian_Fit(xData,yData, cen, ifPlot=False):

    gmod = Model(gaussian)
    result = gmod.fit(yData, x=xData, amp=5, cen=cen, wid=1)
    
    if ifPlot == True : 
        plt.plot(xData, yData,         'bo')
        plt.plot(xData, result.best_fit, 'r-')
        plt.ylim(min(min(result.best_fit),min(xData)),max(max(result.best_fit),max(yData)))
        plt.show()
    
    return result.best_values, result.best_fit

def fit_D_phi_alt(phi_, D, c1, c2, c1o, c2o):
    """fit stuff"""
    Dmod = Model(Dhopalt)
    params = Dmod.make_params()
    params['c1'].set(c1, vary=False)
    params['c2'].set(c2, vary=False)
    params['c1o'].set(c1o, vary=False)
    params['c2o'].set(c2o, vary=False)
    params['c3'].set(0.25)
    result = Dmod.fit(D, params, phi=phi_)
    print result.fit_report(min_correl=0.5)
    c3 = result.params['c3'].value
    return c3

def fit_Dinf_model(Input, diffusion, startindex, endindex):
    time_ = diffusion[startindex:endindex,0]
    D = diffusion[startindex:endindex,1]
    Dmod = Model(Dinf_model)
    params = Dmod.make_params()
    params['var1'].set(0.1) 
    params['Dinf'].set(0.0001)#, min = 0, max=0.1)
    result = Dmod.fit(D, params, time=time_)
    print result.fit_report(min_correl=0.5)
    var1 = result.params['var1'].value
    Dinf = result.params['Dinf'].value
    Dinfstderr = result.params['Dinf'].stderr
    return var1, Dinf, Dinfstderr

def fit_D_phi(phi_, D, c1, c2):
    """fit stuff"""
    Dmod = Model(Dhop)
    params = Dmod.make_params()
    params['c1'].set(c1, vary=False)
    params['c2'].set(c2, vary=False)
    params['c3'].set(1, vary=False)
    params['beta'].set(0.70, min=0)
    result = Dmod.fit(D, params, phi=phi_)
    print result.fit_report(min_correl=0.5)
    beta = result.params['beta'].value
    c3 = result.params['c3'].value
    return c3, beta