I wrote a program to fit some Raman spectra peaks. I need to return the fitted parameters (position, amplitude, HWHM).
I used the modul lmfit to create a lorentzian peak with constraints.
I have a good agreement between the fitted peak and the raw data according to my figure plot. But when it comes to exracte the parameters after fitting, I have a problem, the programm only return the initial values.
I tied the 'report_fit module' and changing the inital parameter without sucess. The parameters values does not evolve.
What bothers me is that this program works on my colleague's PC but not on mine. So the problem may comes from my python version.
I am using spyder 2.3.9, and python 3.4 installed with anaconda under windows 10. The lmfit module 0.9.3, seems to work partially, since I can get a good fitting agreament (from the figure plt.plot). But I cannot return the parameter values after fitting.
Here is my code:
#!/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()
Thank you for you help :)
Please send your reply to deniz.cakir@etu.umontpellier.fr
Your script could be a lot shorter to illustrate the question. But, the final Parameters are held in result_150.params and so forth. The parameters passed into minimize() are not altered in the fit.
Well, that's true for lmfit version 0.9.0 and later. Perhaps your colleague has an older version on lmfit.
