samedi 5 décembre 2020

curve fitting with condition

I want to get the best fit of multiple files at a time. I have more than 700 files with two different feature some of them can be fitted with a single frequency equation [(amp * sin(x2piomega+phase)) * exp(-xdecay)] and the remaining files can be with the (first harmonics) equation [(amp * np.sin(x2np.piomega+phase)+ ampbnp.sin(x4np.piomega+phaseb)) * np.exp(-xdecay)]. So I want to pass all files to the first equation and other files that cannot be nicely fitted with the first equation(this means the frequency is very high or low compared to a neighbor) automatically pass to the second equation to get the best fit. My aim is to get a very small standard deviation of the mean frequency. Here I have applied the if-else condition but I could not get the expected result. So I think there is some error with the conditional command. Here is my code (I want to fit those green peaks(in figure)with second formula)

import numpy as np
from numpy import loadtxt
import lmfit
import glob
import os

def resid(params, x, ydata):
   decay = params['decay'].value
   phase = params['phase'].value
   omega = params['omega'].value
   amp = params['amp'].value
   ampb=params['ampb'].value
   phaseb=params['phaseb'].value
   
   if params['omega']> 9.70 and params['omega'].value< 10.40:
       y_model =  (amp * np.sin(x*2*np.pi*omega+phase)) * np.exp(-x*decay)
       return y_model - ydata
            
   else:
       y_model =  (amp * np.sin(x*2*np.pi*omega+phase)+ ampb*np.sin(x*4*np.pi*omega+phaseb)) * np.exp(-x*decay)
       return y_model - ydata 
       
       
       
filelist = glob.glob("C:/Users/USER/Desktop/more_data/Lower_frequencies/*.txt") 
#print(filelist)
for file in filelist:
   data = np.loadtxt(file)
   x = data[:,0]
   y = data[:,1]
   file_name=file
   temp = os.path.splitext(file_name)[0] 
   params = lmfit.Parameters()
   params.add('phase',    0.0,      min=-np.pi, max=np.pi)
   params.add('omega',  10.32,      min=9,      max=12)
   params.add('amp',      0.2,      min=0,      max=5.0)
   params.add('decay',      4,      min=0,      max=10.0)
   params.add('ampb',    0.03,      min=0.0,    max=5)
   params.add('phaseb',   0.0,      min=-np.pi, max=np.pi)
 
   
   fit = lmfit.minimize(resid, params, args=(x, y), method='least_squares')
   params = {}
   for name, param in fit.params.items():
       elif name == 'omega':
           params['omega'] = param.value
           params['omega_stderr'] = param.stderr
       
       e
   f = open('c:/Users/USER/Desktop/outputA2.txt','a')
   print(
       os.path.basename(temp),
       params['omega'],
       params['decay'],
       fit.success,
       fit.nfev,
       sep="\t", 
     file=f
   )

This figure is the comparison between all files with the first equation(green line) and then all the files fitted with the second equation(red line). My aim is to pass some files automatically to the second equation which have high and low peak values in the first fit (green lines).

Aucun commentaire:

Enregistrer un commentaire