IPKM/mdevaluate
6
1
Fork 0
Code Issues 1 Pull Requests 1 Activity

simplified quick non gaussian fit function so that it is more robust

Browse Source
This commit is contained in:
robrobo
2025-07-17 18:44:20 +02:00
parent dae2d6ed95
commit 0eff84910b
1 changed files with 10 additions and 11 deletions
Download Patch File Download Diff File Expand all files Collapse all files

21
src/mdevaluate/utils.py
View File

@@ -357,14 +357,15 @@ def quick1etau(t: ArrayLike, C: ArrayLike, n: int = 7) -> float:
return tau_est return tau_est
def quicknongaussfit(t, C, width=4): def quicknongaussfit(t, C, width=2):
""" """
Estimates the time and height of the peak in the non-Gaussian function. Estimates the time and height of the peak in the non-Gaussian function.
C is C(t) the correlation function C is C(t) the correlation function
""" """
# TODO this is a very experimental interpolation, can fail def ffunc(t,y0,A_main,log_tau_main,sig_main):
def ppoly(x,a,b,c,d,e,A,mu,sig): main_peak = A_main*np.exp(-(t - log_tau_main)**2 / (2 * sig_main**2))
return A*np.exp(-(x - mu)**2 / (2 * sig**2))+a+(b*x+e)*1/(1+np.exp(c*(x-d))) return y0 + main_peak
# first rough estimate, the closest time. This is returned if the interpolation fails! # first rough estimate, the closest time. This is returned if the interpolation fails!
tau_est = t[np.argmax(C)] tau_est = t[np.argmax(C)]
nG_max = np.amax(C) nG_max = np.amax(C)
@@ -375,13 +376,11 @@ def quicknongaussfit(t, C, width=4):
tau = time[np.argmax(corr)] tau = time[np.argmax(corr)]
mask = (time>tau-width/2) & (time<tau+width/2) mask = (time>tau-width/2) & (time<tau+width/2)
time = time[mask] ; corr = corr[mask] time = time[mask] ; corr = corr[mask]
guess = [0.001,-0.001,5,tau-0.5,1.0,nG_max, tau, 1.5] nG_min = C[t > 0].min()
popt = curve_fit(ppoly, time, corr, p0=guess, maxfev=10000)[0] guess = [nG_min, nG_max-nG_min, tau, 0.6]
# TODO instead use some root or solve popt = curve_fit(ffunc, time, corr, p0=guess, maxfev=10000)[0]
xspace = np.linspace(*time[[0,-1]], 10000) tau_est = 10**popt[-2]
y = ppoly(xspace, *popt) nG_max = popt[0] + popt[1]
tau_est = xspace[np.argmax(y)]
nG_max = np.amax(y)
except: except:
pass pass
if np.isnan(tau_est): if np.isnan(tau_est):