Python lmfit.Minimizer类代码示例

def test_multidimensional_fit_GH205():
    # test that you don't need to flatten the output from the objective
    # function. Tests regression for GH205.
    pos = np.linspace(0, 99, 100)
    xv, yv = np.meshgrid(pos, pos)
    f = lambda xv, yv, lambda1, lambda2: (np.sin(xv * lambda1)
                                             + np.cos(yv * lambda2))

    data = f(xv, yv, 0.3, 3)
    assert_(data.ndim, 2)

    def fcn2min(params, xv, yv, data):
        """ model decaying sine wave, subtract data"""
        lambda1 = params['lambda1'].value
        lambda2 = params['lambda2'].value
        model = f(xv, yv, lambda1, lambda2)
        return model - data

    # create a set of Parameters
    params = Parameters()
    params.add('lambda1', value=0.4)
    params.add('lambda2', value=3.2)

    mini = Minimizer(fcn2min, params, fcn_args=(xv, yv, data))
    res = mini.minimize()

def autobk(energy, mu, rbkg=1, nknots=None, group=None, e0=None,
           kmin=0, kmax=None, kw=1, dk=0, win=None, vary_e0=True,
           chi_std=None, nfft=2048, kstep=0.05, _larch=None):
    if _larch is None:
        raise Warning("cannot calculate autobk spline -- larch broken?")

    # get array indices for rkbg and e0: irbkg, ie0
    rgrid = np.pi/(kstep*nfft)
    if rbkg < 2*rgrid: rbkg = 2*rgrid
    irbkg = int(1.01 + rbkg/rgrid)
    if e0 is None:
        e0 = find_e0(energy, mu, group=group, _larch=_larch)
    ie0 = _index_nearest(energy, e0)

    # save ungridded k (kraw) and grided k (kout)
    # and ftwin (*k-weighting) for FT in residual
    kraw = np.sqrt(ETOK*(energy[ie0:] - e0))
    if kmax is None:
        kmax = max(kraw)
    kout  = kstep * np.arange(int(1.01+kmax/kstep))
    ftwin = kout**kw * ftwindow(kout, xmin=kmin, xmax=kmax,
                                window=win, dx=dk)

    # calc k-value and initial guess for y-values of spline params
    nspline = max(4, min(60, 2*int(rbkg*(kmax-kmin)/np.pi) + 1))
    spl_y  = np.zeros(nspline)
    spl_k  = np.zeros(nspline)
    for i in range(nspline):
        q = kmin + i*(kmax-kmin)/(nspline - 1)
        ik = _index_nearest(kraw, q)
        i1 = min(len(kraw)-1, ik + 5)
        i2 = max(0, ik - 5)
        spl_k[i] = kraw[ik]
        spl_y[i] = (2*mu[ik] + mu[i1] + mu[i2] ) / 4.0
    # get spline represention: knots, coefs, order=3
    # coefs will be varied in fit.
    knots, coefs, order = splrep(spl_k, spl_y)

    # set fit parameters from initial coefficients
    fparams = Parameters()
    for i, v in enumerate(coefs):
        fparams.add("c%i" % i, value=v, vary=i<len(spl_y))

    fitkws = dict(knots=knots, order=order, kraw=kraw, mu=mu[ie0:],
                  irbkg=irbkg, kout=kout, ftwin=ftwin, nfft=nfft)
    # do fit
    fit = Minimizer(__resid, fparams, fcn_kws=fitkws)
    fit.leastsq()

    # write final results
    coefs = [p.value for p in fparams.values()]
    bkg, chi = spline_eval(kraw, mu[ie0:], knots, coefs, order, kout)
    obkg  = np.zeros(len(mu))
    obkg[:ie0] = mu[:ie0]
    obkg[ie0:] = bkg
    if _larch.symtable.isgroup(group):
        setattr(group, 'bkg',  obkg)
        setattr(group, 'chie', mu-obkg)
        setattr(group, 'k',    kout)
        setattr(group, 'chi',  chi)

    def fit_single_line(self, x, y, zero_lev, err_continuum, fitting_parameters, bootstrap_iterations = 1000):

        #Simple fit
        if self.fit_dict['MC_iterations'] == 1:
            fit_output = lmfit_minimize(residual_gauss, fitting_parameters, args=(x, y, zero_lev, err_continuum))
            self.fit_dict['area_intg'] = simps(y, x) - simps(zero_lev, x)
            self.fit_dict['area_intg_err'] = 0.0
             
        #Bootstrap
        else:
            mini_posterior  = Minimizer(lnprob_gaussCurve, fitting_parameters, fcn_args = ([x, y, zero_lev, err_continuum]))
            fit_output      = mini_posterior.emcee(steps=200, params = fitting_parameters)
            
            #Bootstrap for the area of the lines
            area_array = empty(bootstrap_iterations) 
            len_x_array = len(x)
            for i in range(bootstrap_iterations):
                y_new =  y + np_normal_dist(0.0, err_continuum, len_x_array)
                area_array[i] = simps(y_new, x) - simps(zero_lev, x)
            self.fit_dict['area_intg'] = mean(area_array)
            self.fit_dict['area_intg_err'] = std(area_array)           
        
        #Store the fitting parameters
        output_params = fit_output.params
        for key in self.fit_dict['parameters_list']:
            self.fit_dict[key + '_norm'] = output_params[key].value
            self.fit_dict[key + '_norm_er'] = output_params[key].stderr
            
        return

    def setup(self):
        self.x = np.linspace(1, 10, 250)
        np.random.seed(0)
        self.y = (3.0 * np.exp(-self.x / 2)
                  - 5.0 * np.exp(-(self.x - 0.1) / 10.)
                  + 0.1 * np.random.randn(len(self.x)))

        self.p = Parameters()
        self.p.add_many(('a1', 4., True, 0., 10.),
                        ('a2', 4., True, -10., 10.),
                        ('t1', 3., True, 0.01, 10.),
                        ('t2', 3., True, 0.01, 20.))

        self.p_emcee = deepcopy(self.p)
        self.p_emcee.add('noise', 0.2, True, 0.001, 1.)

        self.mini_de = Minimizer(Minimizer_Residual,
                                 self.p,
                                 fcn_args=(self.x, self.y),
                                 kws={'seed': 1,
                                      'polish': False,
                                      'maxiter': 100})

        self.mini_emcee = Minimizer(Minimizer_lnprob,
                                    self.p_emcee,
                                    fcn_args=(self.x, self.y))

def NIST_Test(DataSet, start='start2', plot=True):

    NISTdata = ReadNistData(DataSet)
    resid, npar, dimx = Models[DataSet]
    y = NISTdata['y']
    x = NISTdata['x']

    params = Parameters()
    for i in range(npar):
        pname = 'b%i' % (i+1)
        cval  = NISTdata['cert_values'][i]
        cerr  = NISTdata['cert_stderr'][i]
        pval1 = NISTdata[start][i]
        params.add(pname, value=pval1)


    myfit = Minimizer(resid, params, fcn_args=(x,), fcn_kws={'y':y},
                      scale_covar=True)

    myfit.prepare_fit()
    myfit.leastsq()

    digs = Compare_NIST_Results(DataSet, myfit, params, NISTdata)

    if plot and HASPYLAB:
        fit = -resid(params, x, )
        pylab.plot(x, y, 'r+')
        pylab.plot(x, fit, 'ko--')
        pylab.show()

    return digs > 2

def polyfit(x,y,u=None):
    '''Determine the weighted least-squares fit for 2nd order polynomial'''

    x = np.asarray(x)    
    y = np.asarray(y)

    if u is not None: 
        u[u==0]=1
        weight = 1./np.asarray(u) 
    else:
        weight = np.ones_like(x)

    params = Parameters()
    params.add('a', value=0)
    params.add('b', value=(y.max()-y.min())/(x.max()-x.min()))
    params.add('c', value=0.0)

    def residual(pars, x, data=None,w=None):
        model = pars['a'].value + pars['b'].value*x + pars['c'].value*x**2
        if data is None:
            return model
        return (model - data) #* w
    
    myfit = Minimizer(residual, params,fcn_args=(x,), fcn_kws={'data':y,'w':weight})
    myfit.leastsq()     
    
    return [params['c'].value,params['b'].value,params['a'].value]

class LmModel(object):
    """ 
    Base class for all models. 
    
    Models take x and y and return 
    """
    def __init__(self,x,y):
        self.x, self.y=x,y
        self.parameters=Parameters()
        self.min=Minimizer(self.residual, self.parameters)
    
    def print_para(self):
        for i in self.parameters.values:
            print i
            
    def func(self,paras):
        raise NotImplementedError
    
    def est_startvals(self):
        raise NotImplementedError
    
    def residual(self,paras):
        return self.func(paras)-self.y
        
    def fit(self):
        self.min.leastsq()
        self.y_model=self.func(self.parameters)

def test_constraints1():
    def residual(pars, x, sigma=None, data=None):
        yg = gaussian(x, pars['amp_g'], pars['cen_g'], pars['wid_g'])
        yl = lorentzian(x, pars['amp_l'], pars['cen_l'], pars['wid_l'])

        model =  yg +  yl + pars['line_off'] + x * pars['line_slope']
        if data is None:
            return model
        if sigma is None:
            return (model - data)
        return (model - data)/sigma


    n = 601
    xmin = 0.
    xmax = 20.0
    x = linspace(xmin, xmax, n)

    data = (gaussian(x, 21, 8.1, 1.2) +
            lorentzian(x, 10, 9.6, 2.4) +
            random.normal(scale=0.23,  size=n) +
            x*0.5)


    pfit = Parameters()
    pfit.add(name='amp_g',  value=10)
    pfit.add(name='cen_g',  value=9)
    pfit.add(name='wid_g',  value=1)

    pfit.add(name='amp_tot',  value=20)
    pfit.add(name='amp_l',  expr='amp_tot - amp_g')
    pfit.add(name='cen_l',  expr='1.5+cen_g')
    pfit.add(name='wid_l',  expr='2*wid_g')

    pfit.add(name='line_slope', value=0.0)
    pfit.add(name='line_off', value=0.0)

    sigma = 0.021  # estimate of data error (for all data points)

    myfit = Minimizer(residual, pfit,
                      fcn_args=(x,), fcn_kws={'sigma':sigma, 'data':data},
                      scale_covar=True)

    myfit.prepare_fit()
    init = residual(myfit.params, x)

    result = myfit.leastsq()

    print(' Nfev = ', result.nfev)
    print( result.chisqr, result.redchi, result.nfree)

    report_fit(result.params)
    pfit= result.params
    fit = residual(result.params, x)
    assert(pfit['cen_l'].value == 1.5 + pfit['cen_g'].value)
    assert(pfit['amp_l'].value == pfit['amp_tot'].value - pfit['amp_g'].value)
    assert(pfit['wid_l'].value == 2 * pfit['wid_g'].value)

def test_peakfit():
    from lmfit.utilfuncs import gaussian
    def residual(pars, x, data=None):
        g1 = gaussian(x, pars['a1'].value, pars['c1'].value, pars['w1'].value)
        g2 = gaussian(x, pars['a2'].value, pars['c2'].value, pars['w2'].value)
        model = g1 + g2
        if data is None:
            return model
        return (model - data)

    n    = 601
    xmin = 0.
    xmax = 15.0
    noise = np.random.normal(scale=.65, size=n)
    x = np.linspace(xmin, xmax, n)

    org_params = Parameters()
    org_params.add_many(('a1', 12.0, True, None, None, None),
                        ('c1',  5.3, True, None, None, None),
                        ('w1',  1.0, True, None, None, None),
                        ('a2',  9.1, True, None, None, None),
                        ('c2',  8.1, True, None, None, None),
                        ('w2',  2.5, True, None, None, None))

    data  = residual(org_params, x) + noise


    fit_params = Parameters()
    fit_params.add_many(('a1',  8.0, True, None, 14., None),
                        ('c1',  5.0, True, None, None, None),
                        ('w1',  0.7, True, None, None, None),
                        ('a2',  3.1, True, None, None, None),
                        ('c2',  8.8, True, None, None, None))

    fit_params.add('w2', expr='2.5*w1')

    myfit = Minimizer(residual, fit_params,
                      fcn_args=(x,), fcn_kws={'data':data})

    myfit.prepare_fit()

    init = residual(fit_params, x)


    myfit.leastsq()

    print(' N fev = ', myfit.nfev)
    print(myfit.chisqr, myfit.redchi, myfit.nfree)

    report_fit(fit_params)

    fit = residual(fit_params, x)
    check_paras(fit_params, org_params)

def test_minizers():
    """
    test scale minimizers except newton-cg (needs jacobian) and
    anneal (doesn't work out of the box).
    """
    methods = ['Nelder-Mead', 'Powell', 'CG', 'BFGS',
               'L-BFGS-B', 'TNC', 'COBYLA', 'SLSQP']
    p_true = Parameters()
    p_true.add('amp', value=14.0)
    p_true.add('period', value=5.33)
    p_true.add('shift', value=0.123)
    p_true.add('decay', value=0.010)

    def residual(pars, x, data=None):
        amp = pars['amp'].value
        per = pars['period'].value
        shift = pars['shift'].value
        decay = pars['decay'].value

        if abs(shift) > pi/2:
            shift = shift - np.sign(shift)*pi
        model = amp*np.sin(shift + x/per) * np.exp(-x*x*decay*decay)
        if data is None:
            return model
        return (model - data)

    n = 2500
    xmin = 0.
    xmax = 250.0
    noise = np.random.normal(scale=0.7215, size=n)
    x     = np.linspace(xmin, xmax, n)
    data  = residual(p_true, x) + noise

    fit_params = Parameters()
    fit_params.add('amp', value=11.0, min=5, max=20)
    fit_params.add('period', value=5., min=1., max=7)
    fit_params.add('shift', value=.10,  min=0.0, max=0.2)
    fit_params.add('decay', value=6.e-3, min=0, max=0.1)

    init = residual(fit_params, x)
    mini = Minimizer(residual, fit_params, [x, data])
    for m in methods:
        print(m)
        mini.scalar_minimize(m, x)

        fit = residual(fit_params, x)

        for name, par in fit_params.items():
            nout = "%s:%s" % (name, ' '*(20-len(name)))
            print("%s: %s (%s) " % (nout, par.value, p_true[name].value))

        for para, true_para in zip(fit_params.values(), p_true.values()):
            check_wo_stderr(para, true_para.value)

def fit_axis(image_nparray2D,axis,minim_method="nelder"):
	axis_data = np.sum(pic_data,axis = 1) if (axis == 0) else np.sum(pic_data,axis = 0)
	axis_points = np.linspace(1,len(axis_data),len(axis_data))
	param_estimates = startparams_estimate(axis_data)
	params_for_fit = Parameters()
	params_for_fit.add('I_zero',value=param_estimates[0],min=0,max=np.amax(axis_data))
	params_for_fit.add('r_zero',value=param_estimates[1],min=1,max=len(axis_data))
	params_for_fit.add('omega_zero',value=param_estimates[2],min=1,max=len(axis_data))
	params_for_fit.add('backgr',value=param_estimates[3])
	fit = Minimizer(residual,params_for_fit,fcn_args=(axis_points,),\
		fcn_kws={"data":axis_data})
	fit_res = fit.minimize(minim_method)
	return (axis_points,axis_data,fit_res)

	def fitevent(self, edat, initguess):
		try:
			dt = 1000./self.Fs 	# time-step in ms.

			# control numpy error reporting
			np.seterr(invalid='ignore', over='ignore', under='ignore')

			ts = np.array([ t*dt for t in range(0,len(edat)) ], dtype='float64')

			self.nStates=len(initguess)
			initRCConst=dt*5.

			# setup fit params
			params=Parameters()

			for i in range(0, len(initguess)):
				params.add('a'+str(i), value=initguess[i][0]) 
				params.add('mu'+str(i), value=initguess[i][1]) 
				if self.LinkRCConst:				
					if i==0:
						params.add('tau'+str(i), value=initRCConst)
					else:
						params.add('tau'+str(i), value=initRCConst, expr='tau0')
				else:
					params.add('tau'+str(i), value=initRCConst)

			params.add('b', value=self.baseMean )
			

			igdict=params.valuesdict()

			optfit=Minimizer(self._objfunc, params, fcn_args=(ts,edat,))
			optfit.prepare_fit()
			result=optfit.leastsq(xtol=self.FitTol,ftol=self.FitTol,maxfev=self.FitIters)

			if result.success:
				tt=[init[0] for init, final in zip(igdict.items(), (result.params.valuesdict()).items()) if init==final]
				if len(tt) > 0:
					self.flagEvent('wInitGuessUnchanged')

				self._recordevent(result)
			else:
				#print optfit.message, optfit.lmdif_message
				self.rejectEvent('eFitConvergence')
		except KeyboardInterrupt:
			self.rejectEvent('eFitUserStop')
			raise
		except InvalidEvent:
			self.rejectEvent('eInvalidEvent')
		except:
	 		self.rejectEvent('eFitFailure')

    def __fit2D(self,minim_method="nelder",rotation=False):
        
        self.__fit_axis(0,minim_method)
        self.__fit_axis(1,minim_method)

        # we first take all the initial parameters from 1D fits
        bgr2D_est = self.axis0fitparams.valuesdict()["backgr"]/len(self.axis0pts)
        x2D_est = self.axis0fitparams.valuesdict()["r_zero"]
        omegaX2D_est = self.axis0fitparams.valuesdict()["omega_zero"]
        y2D_est = self.axis1fitparams.valuesdict()["r_zero"]
        omegaY2D_est = self.axis1fitparams.valuesdict()["omega_zero"]

        smoothened_image = gaussian_filter(self.image_array,50)
        peakheight2D_est = np.amax(smoothened_image)
        #now we need to programatically cut out the region of interest out of the
        #whole picture so that fitting takes way less time

        # NOTE! In this implementation, if the beam is small compared to picture size
        # and is very close to the edge, the fitting will fail, because the x and y
        # center position estimates will be off

        self.__format_picture(x2D_est,omegaX2D_est,y2D_est,omegaY2D_est)
        cropped_data = self.formatted_array
        xvals = np.linspace(1,cropped_data.shape[0],cropped_data.shape[0])
        yvals = np.linspace(1,cropped_data.shape[1],cropped_data.shape[1])
        x, y = np.meshgrid(yvals,xvals)
        # NOTE! there's apparently some weird convention, this has to do with
        # Cartesian vs. matrix indexing, which is explain in numpy.meshgrid manual

        estimates_2D = Parameters()
        estimates_2D.add("I_zero",value=peakheight2D_est,min=bgr2D_est)
        estimates_2D.add("x_zero",value=0.5*len(yvals),min=0,max=len(yvals)) # NOTE! weird indexing conventions
        estimates_2D.add("y_zero",value=0.5*len(xvals),min=0,max=len(xvals)) # NOTE! weird indexing conventions
        estimates_2D.add("omegaX_zero",value=omegaX2D_est)
        estimates_2D.add("omegaY_zero",value=omegaY2D_est)
        estimates_2D.add("theta_rot",value=0*np.pi,min = 0,max = np.pi) #just starting with 0
        estimates_2D.add("backgr",value=bgr2D_est)


        if rotation:
            fit2D = Minimizer(residual_G2D,estimates_2D,fcn_args=(x,y),fcn_kws={"data":cropped_data})
            print("Including rotation")
        else:
            fit2D = Minimizer(residual_G2D_norotation,estimates_2D,fcn_args=(x,y),fcn_kws={"data":cropped_data})
            print("Not including rotation")

        fit_res2D = fit2D.minimize(minim_method)

        self.x2Dgrid = x
        self.y2Dgrid = y
        self.fit2Dparams = fit_res2D.params

    def fit(self, image):
        """Fit a image of a hologram with the current attribute 
        parameters.

        Example:
        >>> p = {'x':0, 'y':0, 'z':100, 'a_p':0.5, 'n_p':1.5, 'n_m':1.337, 
        ...      'mpp':0.135, 'lamb':0.447}
        >>> mie_fit = Mie_Fitter(p)
        >>> mit_fit.result(image)
        """
        dim = image.shape
        minner = Minimizer(mie_loss, self.p, fcn_args=(image, dim))
        self.result = minner.minimize()
        return self.result

 def __call__(self):
     #out = minimize(self.residual, 
     #               self.params,
     #               scale_covar = False
     #              #method = 'cg'
     #               )
     mini = Minimizer(self.residual, self.params)
     out = mini.emcee(burn = 10000, steps  = 60000, thin = 1, workers = 1, params = self.params)
     self.H0 = 10**(out.params['a_nu'].value + 5 + 
                    0.2 * (out.params['m04258'].value - 5*log10(out.params['mu_geometric'].value) - 25))
     #print  5*log10(out.params['mu_geometric'].value) + 25
     self.e_H0 = model.H0 * sqrt((out.params['a_nu'].stderr * log(10))**2 
                         + (log(10)/5 *out.params['m04258'].stderr )**2
                         + (out.params['mu_geometric'].stderr/out.params['mu_geometric'].value)**2)
     return out