Errors when using lmfit library

61 Views Asked by At

I am have experimental data for the double slit experiment that I wish to fit a numerical curve to. I have opted to using the lmfit library because I like their fit result output a lot more than other curve fitting packages available to python. To start, I will post my code, then go on to explain the errors I am seeing and what I have read to try and better explain my problem

So here is my code

from PIL import Image
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy.signal import argrelextrema
from lmfit.models import Model

#Determine the observed intensity distribution

img_data = pd.read_csv('Image31 intensity profile.csv')

pixels = np.array(img_data['Distance_(_)'])
counts = np.array(img_data['Gray_Value'])

norm_intensity = counts / np.linalg.norm(counts)


# Double slit interferometry coherent beam

def keV2m(keV):
    """Converts Photon Energy [keV] to wavelength [m].
    .. note:: Calculation in Vacuum.
    """
    wl = 1./(keV*1000)*4.1356*(10**(-7))*2.9998
    return wl

ener = 10#kev
I_0 = 1
a = 2.0e-6 # slit width in  m
d = 10e-6 # slit seperation in m
L = 6.0 # slit to detec dist in m
detec_pixel_size = 3.1e-6 #pixel size

x = np.arange(-500,500)*detec_pixel_size + 1e-9

def coh_inter_patt(x):
    lam = keV2m(ener)

    sinq_part = ((np.sin((np.pi*x*a)/(lam*L))/((np.pi*x*a)/(lam*L)))**2)

    sinq_part[sinq_part == np.nan] = 1

    inter_pat = I_0*(np.cos((np.pi*x*d)/(lam*L))**2)*sinq_part

    return inter_pat

def coh_inter_patt_sft(x,detec_pixel_size,sft=10):
    lam = keV2m(ener)

    sinq_part = ((np.sin((np.pi*(x+(sft*detec_pixel_size))*a)/(lam*L))/((np.pi*(x+(sft*detec_pixel_size))*a)/(lam*L)))**2)

    sinq_part[sinq_part == np.nan] = 1

    inter_pat = I_0*(np.cos((np.pi*(x+(sft*detec_pixel_size))*d)/(lam*L))**2)*sinq_part

    return inter_pat


inter_pat = coh_inter_patt(x)
inter_pat_sft = coh_inter_patt_sft(x,detec_pixel_size,sft=10)

max_vals = np.sort(argrelextrema(inter_pat, np.greater)[0])

min_vals = np.sort(argrelextrema(inter_pat, np.less)[0])

visibility = (inter_pat[max_vals] - inter_pat[min_vals[1:]])/(inter_pat[max_vals] + inter_pat[min_vals[1:]])


# 'Blurred' patterns
extension = 40 # 10*detec_pixel_size long extended source
blur_patt  = np.zeros_like(coh_inter_patt(x))

for i in range(extension):

    blur_patt += coh_inter_patt_sft(x,detec_pixel_size,sft=i)
    blur_patt += coh_inter_patt_sft(x,detec_pixel_size,sft=-i)

blur_patt = blur_patt/np.linalg.norm(blur_patt)

max_vals = np.sort(argrelextrema(inter_pat, np.greater)[0])
min_vals = np.sort(argrelextrema(inter_pat, np.less)[0])

visibility = (inter_pat[max_vals] - inter_pat[min_vals[1:]])/(inter_pat[max_vals] + inter_pat[min_vals[1:]])

blur_max_vals = np.sort(argrelextrema(blur_patt, np.greater)[0])
blur_min_vals = np.sort(argrelextrema(blur_patt, np.less)[0])
blur_visibility = (blur_patt[max_vals] - blur_patt[min_vals[1:]])/(blur_patt[max_vals] + blur_patt[min_vals[1:]])

#Look at plot output
fig = plt.figure(figsize=(10,6))
ax1 = fig.add_subplot(111)
ax1.plot(pixels, norm_intensity, color='red', label='Observed')
ax1.plot(blur_patt,label='Analytical, vis %f'%blur_visibility.mean(), color='k')
ax1.legend()


#iterative curve fitting

def Blurred_fringes(N):
  extension = N # 10*detec_pixel_size long extended source
  blur_patt  = np.zeros_like(coh_inter_patt(x))

  for i in range(extension):

    blur_patt += coh_inter_patt_sft(x,detec_pixel_size,sft=i)
    blur_patt += coh_inter_patt_sft(x,detec_pixel_size,sft=-i)

  #inter_pat = inter_pat/np.linalg.norm(inter_pat)
  blur_patt = blur_patt/np.linalg.norm(blur_patt)
  return blur_patt

model = Model(Blurred_fringes)

params1 = model.make_params(N=35)
results1 = model.fit(counts, params1, x=pixels)

Any variable or calculation involving maxima, minima, and visibility are a seperate calculation and not involved in my problem as far I can tell. So my issue with the code that I posted is that I run the error, "Key Error: N". I know in general that the lmfit fit() function wants dictionary objects made from the make_params() function. And I know from this article here that the first argument to the model function is considered an independent variable by default

To try and work around that, I changed my curve fitting block to the following:

X = np.arange(0,1000,1)

def Blurred_fringes(x, N):
  extension = N # 10*detec_pixel_size long extended source
  blur_patt  = np.zeros_like(coh_inter_patt(x))

  for i in range(int(extension)):

    blur_patt += coh_inter_patt_sft(x,detec_pixel_size,sft=i)
    blur_patt += coh_inter_patt_sft(x,detec_pixel_size,sft=-i)

  #inter_pat = inter_pat/np.linalg.norm(inter_pat)
  blur_patt = blur_patt/np.linalg.norm(blur_patt)
  return blur_patt



test_y = Blurred_fringes(x, 20)

plt.plot(X,test_y)
plt.show()


model = Model(Blurred_fringes)

params1 = model.make_params(N=35)
results1 = model.fit(counts, params1, x=pixels)

That is, I defined a numpy arange that covers the same x range as before, and fed that into my function and fitting commands. And I am instead met with this error "ValueError: The model function generated NaN values and the fit aborted! Please check your model function and/or set boundaries on parameters where applicable.". I believe these NaN values are associated with these output warnings:

:30: RuntimeWarning: invalid value encountered in divide sinq_part = ((np.sin((np.pixa)/(lamL))/((np.pixa)/(lamL)))2) :41: RuntimeWarning: invalid value encountered in divide sinq_part = ((np.sin((np.pi(x+(sftdetec_pixel_size))a)/(lamL))/((np.pi*(x+(sft*detec_pixel_size))a)/(lamL)))**2)

In general, you may run into divide by 0 errors with sinc functions. But I sorted that out with sinq_part[sinq_part == np.nan] = 1.

All in all, my end goal is to get a fit working on this function before I add an iteration step where I look at the confidence levels for a range of N values and find the best fit that way. N should be the only fitting parameter in my case. All of the other variables defined much ealier are experimental parameters and will be fixed (may change slightly due to associated error within those values, for example the sampe-to-detector distance might be 6.2 instead of 6.0). If anyone has any insight on how to proceed, please let me know. Thanks!

0

There are 0 best solutions below