Python models.GaussianModel类代码示例

def gaussFitFunction(fitOut,minimum,maximum,ptNumber=100,prefix=None):
	x = np.linspace(minimum,maximum,ptNumber)
	gaussian = GaussianModel()
	if prefix==None: params = fitOut.best_values
	else: params = getKeysWithoutPrefix(fitOut.best_values,prefix)
	gauss = gaussian.func(x=x,center=params["center"],sigma=params["sigma"],amplitude=params["amplitude"])
	return x, gauss

def test_2gaussians():

  x = np.linspace(0.0, 10.0, num=1000)
  y = gaussian(x, -1, 3, 0.75) + gaussian(x, -0.5, 5, 0.8) + np.random.normal(0, 0.01, x.shape[0])

  gauss1  = GaussianModel(prefix='g1_')
  
  pars = gauss1.guess(y, x=x)
  pars['g1_amplitude'].set(-0.9)
  pars['g1_center'].set(2.5)
  pars['g1_sigma'].set(0.5)

  gauss2  = GaussianModel(prefix='g2_')
  pars.update(gauss2.make_params())
  pars['g2_amplitude'].set(-0.4)
  pars['g2_center'].set(5)
  pars['g2_sigma'].set(0.5)

  mod = gauss1 + gauss2

  init = mod.eval(pars, x=x)

  plt.plot(x, y)
  plt.plot(x, init, 'k--')

  out = mod.fit(y, pars, x=x)

  print(out.fit_report(min_correl=0.5))

  plt.plot(x, out.best_fit, 'r-')
  plt.show()

 def peakFit(self, x, y, model = 'gau', pi = None, di = None):
     ti = time.time()
     if pi:
         NumPeaks = len(pi)
         center = []
         fwhm = []
         amp = []
         numVal = len(x)
         for i in range(NumPeaks):
             pImin = pi[i]-di
             if pImin < 0:
                 pImin = 0
             pImax = pi[i] + di
             if pImax > (numVal-1):
                 pImax = numVal-1
             __y = y[pImin:pImax]
             __x = x[pImin:pImax]
             
             __y = np.power(10,__y/10) #np.array(y)- np.min(y)
             
             mod = GaussianModel()
             pars = mod.guess(__y, x=__x)
             out  = mod.fit(__y, pars, x=__x)
             center.append(out.best_values['center'])
             fwhm.append(out.best_values['sigma']*2.3548)
             amp.append(out.best_values['amplitude'])
         #print 'fit:', time.time()-ti
         return center, fwhm ,amp

def lmfit_ngauss(x,y, *params):
  params = params[0]
  mods = []
  prefixes = []
  for i in range(0, len(params), 3):
    pref = "g%02i_" % (i/3)
    gauss_i = GaussianModel(prefix=pref)

    if i == 0:
      pars = gauss_i.guess(y, x=x)
    else:
      pars.update(gauss_i.make_params())

    A = params[i]
    l_cen = params[i+1]
    sigma = params[i+2]

    pars[pref+'amplitude'].set(A)
    pars[pref+'center'].set(l_cen)
    pars[pref+'sigma'].set(sigma)

    mods.append(gauss_i)
    prefixes.append(pref)

  mod = mods[0]

  if len(mods) > 1:
    for m in mods[1:]:
      mod += m

  print mod

  init = mod.eval(pars, x=x)
  out = mod.fit(y, pars, x=x)
  return mod, out, init

def GaussConst(signal, guess):
    
    amp, centre, stdev, offset = guess
    
    data = np.array([range(len(signal)), signal]).T
    X = data[:,0]
    Y = data[:,1]

    gauss_mod = GaussianModel(prefix='gauss_')
    const_mod = ConstantModel(prefix='const_')
    
    pars = gauss_mod.make_params(center=centre, sigma=stdev, amplitude=amp)
    pars += const_mod.guess(Y, x=X)
    pars['gauss_center'].min = centre - 5.
    pars['gauss_center'].max = centre + 5.
    pars['gauss_sigma'].max = stdev + 5.
    
    mod = gauss_mod + const_mod
    result = mod.fit(Y, pars, x=X)
    
    fwhm = result.best_values['gauss_sigma'] #* 2.3548
    centr = result.best_values['gauss_center']
    
    # Values within two stdevs i.e. 95%
    pl.plot(np.repeat(centr - fwhm * 2, len(Y)),
            np.arange(len(Y)), 'b-')
    pl.plot(np.repeat(centr + fwhm * 2, len(Y)),
            np.arange(len(Y)), 'b-')
    
    return X, result.best_fit, result.best_values['gauss_sigma'] * 4

def fit_profile(profile, guess):
    "Fit a profile to a Gaussian + Contant"
    x = np.arange(len(profile))

    model = GaussianModel(missing='drop') + ConstantModel(missing='drop')
    result = model.fit(profile, x=x, verbose=False, **guess)
    return result

def GaussConst(signal, guess):
    
    if guess == False:
        return [0, 0, 0]
    else:
        amp, centre, stdev, offset = guess
        
        data = np.array([range(len(signal)), signal]).T
        X = data[:,0]
        Y = data[:,1]

        gauss_mod = GaussianModel(prefix='gauss_')
        const_mod = ConstantModel(prefix='const_')
        
        #pars = lorentz_mod.make_params(amplitude=amp, center=centre, sigma=stdev / 3.)
        #lorentz_mod.set_param_hint('sigma', value = stdev / 3., min=0., max=stdev)
        
        pars = gauss_mod.guess(Y, x=X, center=centre, sigma=stdev / 3., amplitude=amp)
        #pars += step_mod.guess(Y, x=X, center=centre)
        pars += const_mod.guess(Y, x=X)
        
        pars['gauss_sigma'].vary = False
        mod = gauss_mod + const_mod
        result = mod.fit(Y, pars, x=X)
        # write error report
        #print result.fit_report()
        fwhm = result.best_values['gauss_sigma'] * 2.3548

        
    return X, result.best_fit, result.redchi, fwhm

def call_gauss(x, y, cen, count, pars):
	label='g'+str(count)+'_'
	gauss = GaussianModel(prefix=label)
	pars.update(gauss.make_params())
	pars[label+'center'].set(cen, min=cen-0.01, max=cen+0.01)
	pars[label+'amplitude'].set(-0.5, min=-10., max=0.0001)
	pars[label+'sigma'].set(0.1, min=0.005, max=0.25)
	return gauss

def gaussian_fit(x, y, title_name):
	mod = GaussianModel()
	pars = mod.guess(y, x=x)
	out = mod.fit(y, pars, x=x)
	plt.figure()
	plt.plot(x, y)
	plt.plot(x, out.best_fit, 'r-')
	plt.title(title_name)
	print(out.fit_report(min_correl = 0.25))
	print('Center at ' + str(out.best_values['center']) + ' Angstrom')
	plt.show()

def get_gaussianmodel(amplitude=1.0, center=5.0, sigma=1.0, noise=0.1):
    # create data to be fitted
    np.random.seed(7392)
    x = np.linspace(-20, 20, 201)
    y = gaussian(x, amplitude, center=center, sigma=sigma)
    y = y + np.random.normal(size=len(x), scale=noise)

    model = GaussianModel()
    params = model.make_params(amplitude=amplitude/5.0,
                               center=center-1.0,
                               sigma=sigma*2.0)
    return x, y, model, params

def xrf_calib_init_roi(mca, roiname):
    """initial calibration step for MCA:
    find energy locations for one ROI
    """
    if not isLarchMCAGroup(mca):
        print( 'Not a valid MCA')
        return
    energy = 1.0*mca.energy
    chans = 1.0*np.arange(len(energy))
    counts = mca.counts
    bgr = getattr(mca, 'bgr', None)
    if bgr is not None:
        counts = counts - bgr
    if not hasattr(mca, 'init_calib'):
        mca.init_calib = OrderedDict()

    roi = None
    for xroi in mca.rois:
        if xroi.name == roiname:
            roi = xroi
            break
    if roi is None:
        return
    words = roiname.split()
    elem = words[0].title()
    family = 'Ka'
    if len(words) > 1:
        family = words[1].title()
    if family == 'Lb':
        family = 'Lb1'
    eknown = xray_line(elem, family)[0]/1000.0
    llim = max(0, roi.left - roi.bgr_width)
    hlim = min(len(chans)-1, roi.right + roi.bgr_width)
    segcounts = counts[llim:hlim]
    maxcounts = max(segcounts)
    ccen = llim + np.where(segcounts==maxcounts)[0]
    ecen = ccen * mca.slope + mca.offset
    bkgcounts = counts[llim] + counts[hlim]
    if maxcounts < 2*bkgcounts:
        mca.init_calib[roiname] = (eknown, ecen, 0.0, ccen, None)
    else:
        model = GaussianModel() + ConstantModel()
        params = model.make_params(amplitude=maxcounts,
                                   sigma=(chans[hlim]-chans[llim])/2.0,
                                   center=ccen-llim, c=0.00)
        params['center'].min = -10
        params['center'].max = hlim - llim + 10
        params['c'].min = -10
        out = model.fit(counts[llim:hlim], params, x=chans[llim:hlim])
        ccen = llim + out.params['center'].value
        ecen = ccen * mca.slope + mca.offset
        fwhm = out.params['fwhm'].value * mca.slope
        mca.init_calib[roiname] = (eknown, ecen, fwhm, ccen, out)

def fluxError(counts, wavelength, error, continuum):
	flux_vector = []
	E_W_vector = []
	cont_avg = np.mean(continuum)
	#plt.close('all')

	for i in range(100):
		#plt.errorbar(wavelength, counts, yerr=error)
		new_counts=[]
		j = 0
		for point in counts:
			new_counts.append(np.random.normal(point, error[j]))
			j = j + 1
		new_counts = np.array(new_counts)
		#So for each N in 1000 a new counts vector is generated randomly
		#Take this data against the wavelength values and fit a gaussian 
		#each time to compute the flux. Append this to a vector and then 
		#find the standard deviation to get the flux error for that emission line
		#Note this is to be encorporated in the fitLines module so that each of the emission 
		#lines is fit in turn. Next step here is to construct the model with lmfit, 
		#guess the initial parameters and then fit the gaussian and take 
		#out.best_values('amplitude') as the flux and store in flux_vector	

		#Now use the lmfit package to perform gaussian fits to the data	
		#Construct the gaussian model
		mod = GaussianModel()

	#Take an initial guess at what the model parameters are 
	#In this case the gaussian model has three parameters, 
	#Which are amplitude, center and sigma
		pars = mod.guess(new_counts, x=wavelength)

	#We know from the redshift what the center of the gaussian is, set this
	#And choose the option not to vary this parameter 
	#Leave the guessed values of the other parameters
		pars['center'].set(value = np.mean(wavelength))


	#Now perform the fit to the data using the set and guessed parameters 
	#And the inverse variance weights form the fits file 
		out  = mod.fit(new_counts, pars, x=wavelength)
		flux = out.best_values['amplitude']
		E_W = out.best_values['amplitude'] / cont_avg
		flux_vector.append(flux)
		E_W_vector.append(E_W)
		#plt.scatter(wavelength, new_counts)
		#plt.plot(wavelength, out.best_fit)
		

	print 'Hello', flux_vector
	#Now return the standard deviation of the flux_vector as the flux error 
	return {'flux_error' : np.std(flux_vector), 'E_W_error' : np.std(E_W_vector)}	

def fit_gaussian(y,x):    
    x=array(x)
    y=array(y)
    mod=GaussianModel()
    pars=mod.guess(y,x=x)
    result=mod.fit(y,pars,x=x)
    a=result.params['amplitude'].value
    b=result.params['center'].value
    c=result.params['sigma'].value
    best=result.best_fit
    chsqred=result.redchi
    chisq=result.chisqr
    fwhm=result.params['fwhm'].value
    
    return a,b,c,best,fwhm,chisq,chsqred

 def peakFit(self, x, y):
     if len(x) == 0:
         y = self.getdBmSpec()
         y = y[self.__scalePos]
         x = self.__scaledWavelength
     y = np.power(10,y/10)
     mod = GaussianModel()
     pars = mod.guess(y, x=x)
     out  = mod.fit(y, pars, x=x)
     
     print(out.fit_report(min_correl=0.25))
     center = out.best_values['center']
     fwhm = out.best_values['sigma']*2.3548
     
     return center, fwhm#, amp

def lmfit_ngauss_constrains(x,y, params, constrains):
  """
  INPUT:
  x - is the wavelength array
  y - is the normalized flux
  params - is a list/array of initial guess values for the parameters
  		   (this controls the number of gaussians to be fitted
  		   	number of gaussians: len(params)/3 - 3 parameters per Gaussian)
  contrains - the limits of the constrains for the fit of the parameters
  OUTPUT:
  mod - the lmfit model object used for the fit
  out - the lmfit fit object that contains all the results of the fit
  init- array with the initial guess model (usefull to see the initial guess when plotting)
  """

  mods = []
  prefixes = []
  for i in range(0, len(params), 3):
    pref = "g%02i_" % (i/3)
    gauss_i = GaussianModel(prefix=pref)

    if i == 0:
      pars = gauss_i.guess(y, x=x)
    else:
      pars.update(gauss_i.make_params())
    A = params[i]
    limA = constrains[i]
    l_cen = params[i+1]
    limL = constrains[i+1]
    sigma = params[i+2]
    limS = constrains[i+2]

    pars[pref+'amplitude'].set(A, min=limA[0], max=limA[1])
    pars[pref+'center'].set(l_cen, min=limL[0], max=limL[1])
    pars[pref+'sigma'].set(sigma, min=limS[0], max=limS[1])

    mods.append(gauss_i)
    prefixes.append(pref)

  mod = mods[0]

  if len(mods) > 1:
    for m in mods[1:]:
      mod += m
  init = mod.eval(pars, x=x)
  out = mod.fit(y, pars, x=x)

  return mod, out, init