from math import inf, ceil import re import numpy as np import pandas as pd import matplotlib.pyplot as plt from matplotlib import style from lmfit import models from lmfit.models import GaussianModel import scipy.stats from scipy.signal import savgol_filter from statsmodels.nonparametric.kernel_regression import KernelReg from scipy.signal import find_peaks from matplotlib import gridspec import random import kneed from scipy.signal import peak_widths def generate_model(spec): composite_model = None params = None x = spec['x'] y = spec['y'] x_min = np.min(x) x_max = np.max(x) x_range = x_max - x_min y_max = np.max(y) for i, basis_func in enumerate(spec['model']): prefix = f'm{i}_' model = getattr(models, basis_func['type'])(prefix=prefix) if basis_func['type'] in ['GaussianModel', 'LorentzianModel', 'VoigtModel']: # for now VoigtModel has gamma constrained to sigma model.set_param_hint('sigma', min=1e-6, max=x_range) model.set_param_hint('center', min=x_min, max=x_max) model.set_param_hint('height', min=1e-6, max=1.1*y_max) model.set_param_hint('amplitude', min=1e-6) # default guess is horrible!! do not use guess() default_params = { prefix+'center': x_min + x_range * random.random(), prefix+'height': y_max * random.random(), prefix+'sigma': x_range * random.random() } else: raise NotImplemented(f'model {basis_func["type"]} not implemented yet') if 'help' in basis_func: # allow override of settings in parameter for param, options in basis_func['help'].items(): model.set_param_hint(param, **options) model_params = model.make_params(**default_params, **basis_func.get('params', {})) if params is None: params = model_params else: params.update(model_params) if composite_model is None: composite_model = model else: composite_model = composite_model + model return composite_model, params def _1gaussian(x, amp1,cen1,sigma1): return amp1*(1/(sigma1*(np.sqrt(2*np.pi))))*(np.exp((-1.0/2.0)*(((x-cen1)/sigma1)**2))) def _2gaussian(x, am1,mu1,sigma1, am2,mu2,sigma2): y = am1/(np.sqrt(2*np.pi*(sigma1**2)))*np.exp(-0.5*(x-mu1)**2/(sigma1**2))+ am2/(np.sqrt(2*np.pi*(sigma2**2)))*np.exp(-0.5*(x-mu2)**2/(sigma2**2)) #y = am/(np.sqrt(2*np.pi*(sigma**2)))*np.exp(-0.5*(x-mu)**2/(sigma**2)) return y def gaussian(x, am, sigma, mu): y = am/(np.sqrt(2*np.pi*(sigma**2)))*np.exp(-0.5*(x-mu)**2/(sigma**2)) return y def modelpre(dframe): mf_rframe = dframe[dframe['x']> -inf].reset_index() modelr = GaussianModel() params_r = modelr.guess(mf_rframe['y'], x=mf_rframe['x']) result_r = modelr.fit(mf_rframe['y'], params_r, x=mf_rframe['x']) sigma_r = result_r.best_values['sigma'] center_r = result_r.best_values['center'] amp_r = result_r.best_values['amplitude'] #residen = abs(dframe[(dframe['x']<=center_r) & (dframe['x'])]['y']-gaussian(dframe[(dframe['x']<=center_r) & (dframe['x'])]['x'],amp_r,sigma_r,center_r)-gaussian(dframe[(dframe['x']<=center_r) & (dframe['x']>=0)]['x'],amp_r,sigma_r,center_r)) #resi_cross = dframe['x'].residen.idxmin()] #crosspoint = abs(dframe['y']-gaussian(dframe['x'],amp_r,sigma_r,center_r)) #cros2x = dframe['x'][crosspoint.idxmin()] #fpr_t = sum(dframe[(dframe['x']>=resi_cross) & (dframe['x']<=cros2x)]['y']-gaussian(dframe[(dframe['x']>=resi_cross) & (dframe['x']<=cros2x)]['x'],amp_r,sigma_r,center_r)) #nor_t = sum(dframe[dframe['x']>=resi_cross]['y']) #fnr = scipy.stats.norm(center_r,sigma_r).cdf(resi_cross) #return(resi_cross, cros2x, fnr, fpr_t, nor_t, amp_r, sigma_r, center_r) return(amp_r, sigma_r, center_r) posx=[-10.0, -9.8, -9.6, -9.4, -9.2, -9.0, -8.8, -8.6, -8.4, -8.2, -8.0, -7.8, -7.6, -7.4, -7.199999999999999, -7.0, -6.8, -6.6, -6.4, -6.199999999999999, -6.0, -5.8, -5.6, -5.3999999999999995, -5.199999999999999, -5.0, -4.8, -4.6, -4.3999999999999995, -4.199999999999999, -4.0, -3.8, -3.5999999999999996, -3.3999999999999995, -3.1999999999999993, -3.0, -2.8, -2.5999999999999996, -2.3999999999999995, -2.1999999999999993, -2.0, -1.799999999999999, -1.5999999999999996, -1.4000000000000004, -1.1999999999999993, -1.0, -0.7999999999999989, -0.5999999999999996, -0.3999999999999986, -0.1999999999999993, 0.0, 0.20000000000000107, 0.40000000000000036, 0.6000000000000014, 0.8000000000000007, 1.0, 1.200000000000001, 1.4000000000000004, 1.6000000000000014, 1.8000000000000007, 2.0, 2.200000000000001, 2.4000000000000004, 2.6000000000000014, 2.8000000000000007, 3.0, 3.200000000000001, 3.4000000000000004, 3.6000000000000014, 3.8000000000000007, 4.0, 4.200000000000001, 4.4, 4.600000000000001, 4.800000000000001, 5.0, 5.200000000000001, 5.4, 5.600000000000001, 5.800000000000001, 6.0, 6.199999999999999, 6.400000000000002, 6.600000000000001, 6.800000000000001, 7.0, 7.199999999999999, 7.400000000000002, 7.600000000000001, 7.800000000000001, 8.0, 8.2, 8.400000000000002, 8.600000000000001, 8.8, 9.0, 9.200000000000003, 9.400000000000002, 9.600000000000001, 9.8] posy=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.00010048231511254056, 0.00010048231511254056, 0.00010048231511254056, 0.00025120578778135003, 0.0004521704180064311, 0.000803858520900323, 0.001105305466237942, 0.0018589228295819962, 0.004421221864951757, 0.009093649517684894, 0.018890675241157596, 0.0405446141479099, 0.08490755627009655, 0.16001808681671995, 0.25859123794212224, 0.35776728295819976, 0.432475884244373, 0.47352290996784574, 0.48829381028938906, 0.48889670418006437, 0.4803054662379421, 0.45910369774919624, 0.4148914790996784, 0.33968046623794246, 0.24095659163987143, 0.14142885852090004, 0.06641881028938913, 0.024618167202572278, 0.007284967845659184, 0.002009646302250807, 0.0007033762057877793, 0.00025120578778135095, 0.00010048231511254011, 5.024115755627028e-05, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] #AGAGCATTGA [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.789272030651358e-05, 4.789272030651358e-05, 4.789272030651358e-05, 0.0002873563218390815, 0.0003831417624521078, 0.0005268199233716486, 0.000670498084291188, 0.0008620689655172423, 0.0009578544061302687, 0.0012452107279693502, 0.0018199233716475133, 0.002969348659003829, 0.004693486590038318, 0.007567049808429105, 0.016139846743295036, 0.031752873563218466, 0.06551724137931024, 0.12437739463601544, 0.21149425287356285, 0.3153256704980843, 0.40086206896551757, 0.45560344827586213, 0.4773946360153258, 0.4836206896551725, 0.47959770114942535, 0.4655651340996168, 0.43242337164750977, 0.3738984674329502, 0.2869252873563222, 0.18323754789272037, 0.09722222222222197, 0.04137931034482764, 0.01786398467432947, 0.008524904214559408, 0.003879310344827591, 0.0021551724137931017, 0.0013888888888888913, 0.0008620689655172419, 0.0006226053639846761, 0.0001915708812260539, 4.7892720306513154e-05, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] posdata = {'x':posx, 'y':posy} posdata = pd.DataFrame(posdata) #print(posdata['x']) #peaks, _ = find_peaks(posy, height=0.1) #peaks, _ = find_peaks(-y) #print(peaks) #peakv = [posx[x] for x in peaks] #print(peaks, peakv) spec = { 'x': posx, 'y': posy, 'model':[ {'type': 'GaussianModel'}, {'type': 'GaussianModel'} ] } def plot_to_blog(fig, figure_name): filename = os.path.expanduser(f'{image_dir}/{figure_name}') fig.savefig(filename) return filename model, params = generate_model(spec) output = model.fit(spec['y'], params, x=spec['x']) print(output.best_values) m0_amp = output.best_values['m0_amplitude'] m0_cen = output.best_values['m0_center'] m0_sigma = output.best_values['m0_sigma'] m1_amp = output.best_values['m1_amplitude'] m1_cen = output.best_values['m1_center'] m1_sigma = output.best_values['m1_sigma'] #output.plot(data_kws={'markersize': 1}) plt.scatter(posdata['x'],gaussian(posdata['x'],m0_amp,m0_sigma,m0_cen),s=2 ,label='best fit m0') plt.scatter(posdata['x'],gaussian(posdata['x'],m1_amp,m1_sigma,m1_cen),s=2 ,label='best fit m1') plt.scatter(posdata['x'],_2gaussian(posdata['x'], m0_amp,m0_cen,m0_sigma,m1_amp,m1_cen,m1_sigma),s=2 ,label='best fit m0+m1') posym = _2gaussian(posdata['x'], m0_amp,m0_cen,m0_sigma,m1_amp,m1_cen,m1_sigma) peaks, _ =find_peaks(posym, height=max(posym)*0.1) _, _,_, posyw = peak_widths(posym, peaks, rel_height=1) print(posyw[0]) kneedle = kneed.KneeLocator(posx[peaks[0]:ceil(posyw[0])], posy[peaks[0]:ceil(posyw[0])], S=1.0, curve="convex", direction="decreasing") kneepos = kneedle.elbow print(posyw) print(peaks) print(kneepos) print(np.trapz(posy[0:70], posx[0:70], dx=1)) plt.plot(posx, posy) plt.vlines(kneepos, 0, 0.1) #fig.show() #plot_to_blog(fig, 'xrd-fitting-two-gaussian-noise-lmfit-spec.png') #popt_2gauss, pcov_2gauss = scipy.optimize.curve_fit(_2gaussian, posx, posy, p0=[1,0,1,2,0,1], maxfev=5000) #perr_2gauss = np.sqrt(np.diag(pcov_2gauss)) #pars_1 = popt_2gauss[-inf:0] #pars_2 = popt_2gauss[0:inf] #gauss_peak_1 = _1gaussian(posx, *pars_1) #gauss_peak_2 = _1gaussian(posx, *pars_2) #fig = plt.figure(figsize=(4,3)) #gs = gridspec.GridSpec(1,1) #ax1 = fig.add_subplot(gs[0]) #ax1.plot(posx, gauss_peak_1, "g") #ax1.fill_between(posx, gauss_peak_1.min(), gauss_peak_1, facecolor="green", alpha=0.5) #ax1.plot(posx, gauss_peak_2, "y") #ax1.fill_between(posx, gauss_peak_2.min(), gauss_peak_2, facecolor="yellow", alpha=0.5) #resi_cross, cros2x, fnr, fpr_t, nor_t, amp_r, sigma_r, center_r = modelpre(posdata) #amp_r, sigma_r, center_r = modelpre(posdata) #print(cros2x) #plt.scatter(posdata['x'],posdata['y']-gaussian(posdata['x'],amp_r,sigma_r,center_r),s=1 ,label='residual') #plt.scatter(posdata['x'],gaussian(posdata['x'],amp_r,sigma_r,center_r),s=2 ,label='best fit') #plt.vlines([resi_cross], 0, max(posy), #linestyles="dashed",linewidth=0.5, colors='black') #plt.plot(posx, posy, label = 'pos_AN') plt.legend() plt.show()