我编写了一个程序来拟合一些拉曼光谱峰。
我需要返回已拟合参数(位置,幅度,半高全宽)。
我使用了lmfit模块创建了一个带约束的洛伦兹峰。
从我的图形绘制中可以看出,我得到了良好的拟合峰和原始数据之间的一致性。但是在提取拟合后的参数时,我遇到了问题,该程序仅返回初始值。
我尝试过“report_fit模块”和更改初始参数,但没有成功。参数值没有发生变化。
让我困扰的是,这个程序在我的同事电脑上可以正常运行,但是在我的电脑上不行。因此,问题可能来自于我的Python版本。
我正在使用Spyder 2.3.9,并在Windows 10下使用Anaconda安装的Python 3.4。lmfit模块0.9.3似乎部分工作,因为我可以得到很好的拟合结果(来自图形plt.plot)。但是我无法返回拟合后的参数值。
这是我的代码:
#!/usr/bin/python3
# -*- coding:utf-8 -*-
import os
import numpy as np
import matplotlib.pyplot as plt
from math import factorial
from scipy.interpolate import interp1d
from lmfit import minimize, Parameters #,report_fit
##############################################################################
# Fit of Raman peaks
def fmin150(pars,x,y):
amp= pars['Amp_150'].value
cent=pars['Cent_150'].value
hwhm=pars['Wid_150'].value
a=pars['a_150'].value
b=pars['b_150'].value
peak = (amp*hwhm)/(((x-cent)**2)+(hwhm**2)) + ((a*x)+b)
return peak - y
def fmin220(pars,x,y):
amp= pars['Amp_220'].value
cent=pars['Cent_220'].value
hwhm=pars['Wid_220'].value
a=pars['a_220'].value
b=pars['b_220'].value
peak = (amp*hwhm)/(((x-cent)**2)+(hwhm**2)) + ((a*x)+b)
return peak - y
def fmin2d(pars,x,y):
amp= pars['Amp_2D'].value
cent=pars['Cent_2D'].value
hwhm=pars['Wid_2D'].value
a=pars['a_2D'].value
b=pars['b_2D'].value
peak = (amp*hwhm)/(((x-cent)**2)+(hwhm**2)) + ((a*x)+b)
return peak - y
##############################################################################
def fit(filename):
"""
Fit Raman spectrum file
Return filename, for each peak (*f*whm, position, height)
PL (Position, Intensity, Area)
"""
print ("----------------------------")
print("Treating file : ")
print(filename)
try:
data = np.loadtxt(filename)
except:
print("Unable to load file")
xx=data[:,0]
yy=data[:,1]
#### Define fitting interval #####
# Cu oxides (unités en cm-1)
xminim150 = 120
xmaxim150 = 170
xminim220 = 175
xmaxim220 = 275
xminim300 = 280
xmaxim300 = 340
xminim640 = 345
xmaxim640 = 800
# Graphene
xminimG = 1540
xmaximG = 1680
xminim2D = 2600
xmaxim2D = 2820
# Definition Baground = place without the fitted peaks
zone1 = (xx > min(xx)) & (xx < xminim150)
zone2 = (xx > xmaxim150) & (xx < xminim220)
zone3 = (xx > xmaxim220) & (xx < xminim300)
zone4 = (xx > xmaxim300) & (xx < xminim640)
zone5 = (xx > xmaxim640) & (xx < xminimG)
zone6 = (xx > xmaximG) & (xx < xminim2D)
zone7 = (xx > xmaxim2D) & (xx < max(xx))
x_BG = np.concatenate((xx[zone1],xx[zone2],xx[zone3],xx[zone4],xx[zone5],xx[zone6],xx[zone7]))
y_BG = np.concatenate((yy[zone1],yy[zone2],yy[zone3],yy[zone4],yy[zone5],yy[zone6],yy[zone7]))
#Creation de l'interpolation lineaire
f1 = interp1d(x_BG, y_BG, kind='linear')
xinterpmin= x_BG[0] # valeur de x_BG min
xinterpmax= x_BG[len(x_BG)-1] # valeur de x_BG max
nbxinterp = len(xx) * 4 #(nb de point en x)* 4 pour une interpolation correcte
xnew = np.linspace(xinterpmin, xinterpmax, num=nbxinterp, endpoint=True)
ynew= f1(xnew)
########################## Fit 2D peaks ###############################
# create a set of Parameters
pars150 = Parameters()
pars220 = Parameters()
pars300 = Parameters()
pars640 = Parameters()
parsg = Parameters()
pars2d = Parameters()
##### Cu2O pic 150 cm-1 #####
pars150.add('Amp_150', value=10, min=0.0001, max=100000) # Amplitude ~intensity
pars150.add('Cent_150', value=150, min=140, max=160) # Center ~position
pars150.add('Wid_150', value=5, min=4, max=15 ) # Width is the HWHM
pars150.add('a_150', value=1, min=-100000, max=100000)
pars150.add('b_150', value=10, min=-100000, max=100000)
##### Cu2O pic 220 cm-1 #####
pars220.add('Amp_220', value=10, min=0.0001, max=100000)
pars220.add('Cent_220', value=220, min=200, max=230)
pars220.add('Wid_220', value=5, min=4, max=15 )
pars220.add('a_220', value=1, min=-100000, max=100000)
pars220.add('b_220', value=10, min=-100000, max=100000)
##### Graphene 2D #####
pars2d.add('Amp_2D', value=15, min=0.0001, max=100000)
pars2d.add('Cent_2D', value=2711, min=2670, max=2730)
pars2d.add('Wid_2D', value=15, min=4, max=25 )
pars2d.add('a_2D', value=1, min=-100000, max=100000)
pars2d.add('b_2D', value=10, min=-100000, max=100000)
# define x for each peaks
interv_150 = (xx > xminim150) & (xx < xmaxim150)
x_150 = xx[interv_150]
y_150 = yy[interv_150]
interv_220 = (xx > xminim220) & (xx < xmaxim220)
x_220 = xx[interv_220]
y_220 = yy[interv_220]
interv_2D = (xx > xminim2D) & (xx < xmaxim2D)
x_2D = xx[interv_2D]
y_2D = yy[interv_2D]
###########################################################
# Performe FIT with leastsq model ###########
result_150 = minimize(fmin150, pars150, args=(x_150, y_150))
result_220 = minimize(fmin220, pars220, args=(x_220, y_220))
result_2D = minimize(fmin2d, pars2d, args=(x_2D, y_2D))
# calculate final result
final_150 = y_150 + result_150.residual
final_220 = y_220 + result_220.residual
final_2D = y_2D + result_2D.residual
###########################################################
# Parameter after fit #
amp_150= pars150['Amp_150'].value
cent_150=pars150['Cent_150'].value
fwhm_150=2*pars150['Wid_150'].value
amp_220= pars220['Amp_220'].value
cent_220=pars220['Cent_220'].value
fwhm_220=2*pars220['Wid_220'].value
amp_2D= pars2d['Amp_2D'].value
cent_2D=pars2d['Cent_2D'].value
fwhm_2D=2*pars2d['Wid_2D'].value
###########
#Plot data#
plt.plot(xx, yy, 'k+' ,x_150, final_150, 'r', x_220, final_220,'r', x_2D, final_2D,'b')
plt.xlabel(r'Raman shift (cm$^{-1}$)', fontsize=14)
plt.ylabel('Intensity (a.u.)', fontsize=14)
plt.xlim(0,3000)
plt.title(filename, fontsize=16)
savename=filename[:-4]+".png"
print(savename)
plt.savefig(savename)
plt.clf()
return filename, amp_150, cent_150, fwhm_150, amp_220, cent_220, fwhm_220, amp_2D, cent_2D, fwhm_2D
def main():
"""main program loop"""
print("----------------------------")
liste = []
for filename in os.listdir(os.getcwd()):
if filename.endswith(".txt"):
liste.append(filename)
f1 = open("TestFit_all.dat","w")
header = "Filename\tI_150\tCentre_150\tFWHM_150\tI_220\tCentre_220\tFWHM_220\tI_300\tCentre_300\tFWHM_300\tI_640\tCentre_640\tFWHM_640\tI_G\tCentre_G\tFWHM_G\tI_2D\tCentre_2D\tFWHM_2D\tI_PL\tI_PL_err\tCent_PL\tCent_PL_err\tArea1000_PL\n"
f1.write(header)
for i in liste:
txt=fit(i)
print(txt)
#text = str(txt[0])+"\t"+str(txt[1])+"\t"+str(txt[2])+"\t"+str(txt[3])+"\t"+str(txt[4])+"\t"+"\n"
text = str(txt[0])+"\t"+str(txt[1])+"\t"+str(txt[2])+"\t"+str(txt[3])+"\t"+str(txt[4])+"\t"+str(txt[5])+"\t"+str(txt[6])+"\t"+str(txt[7])+"\t"+str(txt[8])+"\t"+str(txt[9])+"\t"+"\t"+"\n"
f1.write(text)
f1.close()
print("----------------------------")
print("Done")
###################################################################
# start
if __name__ == "__main__":
main()
感谢您的帮助 :)
一个有效的适配示例
请将您的回复发送至deniz.cakir@etu.umontpellier.fr