# -*- coding: utf-8 -*- """ Created on Thu Aug 31 12:26:03 2023 @author: Willger """ import numpy as np from scipy.interpolate import interp1d from scipy.optimize import least_squares,minimize import pandas as pd from matplotlib import pyplot as plt from scipy.signal import find_peaks from scipy import interpolate def find_marker(mixturename): colors= plt.get_cmap('jet',16) if mixturename=='Acetone-Cyclohexane': marker='o'#$CyHex$' color='tab:grey' elif mixturename=='Acetone-Water': marker='D'#$H2O$' color='tab:blue' elif mixturename=='Acetone-Acetic acid': marker='s'#'Ac.ac.' color='tab:green' elif mixturename=='Acetone-ACN': marker='p'#'tab:orange' color=colors(11) return marker,color #%% preprocessing def SlicedRamanshift(ramanshift,idx_start, idx_end): """ For extracting the required region of a ramanshift. Parameters ---------- ramanshift : List The ramanshift for each spectrum. idx_start : Integer Starting index of the spectrum and ramanshift. idx_end : Integer Ending index of the spectrum and ramanshift. Returns ------- sliced_ramanshift : List The sliced ramanshift according to the starting and ending index choosen. """ sliced_ramanshift = ramanshift[idx_start:idx_end] return(sliced_ramanshift) def baseline_correction(spectra, ramanshift, basepoints): """ baseline correction Parameters ---------- spectra : array raw spectrum ramanshift : array ramashift basepoints : array basepoints Returns ------- list baseline corrected spectrum """ if isinstance(spectra, list): spectra = np.array(spectra) idx_sample = np.zeros(len(basepoints)) val_sample = np.zeros(len(basepoints)) sampling_with = 2 for i in range(len(basepoints)): idx_sample[i] = int((np.abs(ramanshift - basepoints[i])).argmin()) val_sample[i] = min(spectra[int(idx_sample[i] - sampling_with):int(idx_sample[i] + sampling_with)]) baseline_func = interp1d(basepoints, val_sample, fill_value='extrapolate',kind='linear') baseline = baseline_func(ramanshift) spectra_bc = spectra - baseline # plt.figure() # plt.plot(ramanshift,baseline) # plt.plot(ramanshift,spectra) # plt.show() return ([spectra_bc]) def norm_spectra(spectra, ramanshift, idx_start,idx_end): if isinstance(spectra, list): spectra = np.array(spectra).flatten() sum = np.sum(spectra[idx_start:idx_end]) spectra_norm = spectra[idx_start:idx_end] / sum return list(spectra_norm) def norm_spectra_c1(spectra,ramanshift, idx_start,idx_end): spectra=np.vstack(np.array(spectra)).astype(np.float64).reshape(len(spectra),1043) spectra_c1=spectra[-1] sum_c1 = np.sum(spectra_c1[idx_start:idx_end]) spectra_norm = spectra[:,idx_start:idx_end] / sum_c1 sum_norm_c1=np.sum(spectra_norm,axis=1) frame={ 'spectra_norm_c1':pd.Series(list(spectra_norm)), 'sum_norm_c1':pd.Series(list(sum_norm_c1)) } return pd.DataFrame(frame) def calculate_int_mol(spectra): spectra=np.vstack(np.array(spectra)).astype(np.float64) int_mol=np.sum(spectra,axis=1) return int_mol def calculate_mol_spectra(x1,k,spectra): spectra=np.vstack(np.array(spectra)).astype(np.float64) k=np.array(k) x1=np.array(x1) spectra_mol=(x1+(1-x1)*k)*spectra.T return pd.DataFrame({'spectra_mol':list(spectra_mol.T)}) def calc_k(ramanshift,data_norm,x1,idx_fit_start,idx_fit_end): spectra=np.vstack(np.array(data_norm)).astype(np.float64) a0=3 fit='Spectral_difference' if fit =='Spectral_difference': k=np.zeros_like(x1) for i in range (len(x1)): funtomin=lambda a_0:spectra[i,:]-(x1[i]/(x1[i]+(1-x1[i])*a_0)*spectra[-1,:]+(1-x1[i]/(x1[i]+(1-x1[i])*a_0))*spectra[0,:]) res=least_squares(funtomin,a0) k[i]=res.x[0] index_mix=np.array(x1.where((x1 !=0)&(x1 !=1)).dropna().index) k_mix=k[index_mix] k_mean=np.array([np.mean(k_mix)]*len(x1)) x1=np.array(x1) spectra_id=((1-x1)*k_mean)*np.array([spectra[0,:],]*len(x1)).T+x1*np.array([spectra[-1,:],]*len(x1)).T else: k=np.zeros_like(x1) for i in range (len(x1)): funtomin_area=lambda a_0:np.abs(np.sum(spectra[i,:])-np.sum((x1[i]/(x1[i]+(1-x1[i])*a_0)*spectra[-1,:]+(1-x1[i]/(x1[i]+(1-x1[i])*a_0))*spectra[0,:]))) res=minimize(funtomin_area,a0) k[i]=res.x[0] k_mean=np.array([np.mean(k)]*len(x1)) x1=np.array(x1) spectra_id=((1-x1)*k_mean)*np.array([spectra[0,:],]*len(x1)).T+x1*np.array([spectra[-1,:],]*len(x1)).T frame={'spectra_id':pd.Series(list(spectra_id.T)), 'k':pd.Series(list(k_mean))} res=pd.DataFrame(frame) return(res) def id_spectra(x1,data_norm): x1=np.array(x1) spectra=np.vstack(np.array(data_norm)).astype(np.float64) spectra_id=np.zeros_like(spectra) for i in range (0,len(spectra)): spectra_id[i,:]=x1[i]*spectra[-1,:]+(1-x1[i])*spectra[0,:] return list(spectra_id) def part_specs(spectra_in,x1): k=2 spectra=np.vstack(np.array(spectra_in)).astype(np.float64) spectra=np.squeeze(spectra) x1=np.array(x1) part_specs_c1=np.zeros_like(spectra) part_specs_c1[:]=np.nan part_specs_c2=np.zeros_like(spectra) part_specs_c2[:]=np.nan n_pix=len(spectra[0,:]) for i in range (len(x1)): if x1[i]==0: pass elif x1[i]==1: part_specs_c1[i,:]=spectra[i,:] else: for j in range(n_pix): if i=len(x1): part_diff=(spectra[i+1,j]-spectra[i-k,j])/(x1[i+1]-x1[i-k]) else: part_diff=(spectra[i+k,j]-spectra[i-k,j])/(x1[i+k]-x1[i-k]) part_specs_c1[i,j]=spectra[i,j]+(1-x1[i])*part_diff for i in range (len(x1)): if (1-x1[i])==0: pass elif (1-x1[i])==1: part_specs_c2[i,:]=spectra[i,:] else: for j in range(n_pix): if i=len(x1): part_diff=(spectra[i+1,j]-spectra[i-k,j])/(x1[i+1]-x1[i-k]) else: part_diff=(spectra[i+k,j]-spectra[i-k,j])/(x1[i+k]-x1[i-k]) part_specs_c2[i,j]=spectra[i,j]-x1[i]*part_diff frame={'part_specs_c1':list(part_specs_c1), 'part_specs_c2':list(part_specs_c2) } return pd.DataFrame(frame) #%% voigt fit OH_end=3800 OH_start=3100 def pseudo_voigt(rs, paras,nr_peaks): """ Calculates pseudo voigt profile of a spectrum or single peak in the intensity form Parameters ---------- rs : array float array: nr. of peaks * ramanshfitvector paras : array of floats parametres of pseudo voigt profile: paras[:,0]: Lorentz percentages paras[:,1]: peakcenters paras[:,2]: peakwidth paras[:,3]: intensity nr_peaks : int nr. of peaks Returns ------- v_profile_bc : array float baeline corrected voigt fit """ 'single profile' if nr_peaks==1: x=rs-paras[1] v_profile=(paras[0]*paras[3]*2/paras[2]*1/np.pi*1/((x/(paras[2]*0.5))**2+1)+(1-paras[0])*paras[3]*2/paras[2]*np.sqrt(np.log(2)/np.pi)*np.exp(-np.log(2)*(x/(0.5*paras[2]))**2)).flatten() else: """multiple profiles""" paras=np.reshape(paras,(nr_peaks,4)) nr_peaks=len(paras[:,1]) x=rs-paras[:,1] v_profile=np.sum(paras[:,0]*paras[:,3]*2/paras[:,2]*1/np.pi*1/((x/(paras[:,2]*0.5))**2+1)+(1-paras[:,0])*paras[:,3]*2/paras[:,2]*np.sqrt(np.log(2)/np.pi)*np.exp(-np.log(2)*(x/(0.5*paras[:,2]))**2),axis=1) return (v_profile) """calculate derivative of spectrum""" def deriv(spectra,rs): """ calculates derivative of spectra according ramanshift Parameters ---------- spectra : float array spectrum rs : float array ramanshift Returns ------- derivative of spectrum """ derivate=np.zeros_like(rs) for i in range(2,len(rs)-1): derivate[i]=(spectra[i+1]-spectra[i-1])/(rs[i+1]-rs[i-1]) return derivate def start_paras_pure(peaks,spectra,rs,max_peakwidth): """ get initial parameters for pure compound voigt profile fit Parameters ---------- peaks : float array peakcenters detected by find_peaks spectra : float array meaured spectrum rs : float array raman shift/cm^-1 Returns ------- bounds_ reshped: array with shape (nr. of peaks, 8); first column containing lower und upper bounds for lorentz-percentage, peakcenter, peakwidth, intensity pars_0: contains initial values for each peak """ nr_peaks=len(peaks) #%% initial values "peakcenter, intensity, lorentz percentage" peakcenter_0=rs[peaks] intensity_0=spectra[peaks]*40 lor_0=np.full((nr_peaks,1),0.5).flatten() "peak width" gamma_0=np.empty_like(intensity_0) for i in range (nr_peaks): # fingerprint region: if rs[peaks[i]]<2000: gamma_0[i]=15 # CH-bond: elif (2000left_pos and spec[left] > left_min+left_fwhm_threshold: left -= 1 left_fwhm_idx = interpolate_spec(left, spec[left], left + 1, spec[left + 1],left_min+ left_fwhm_threshold) #if left > 0 else left left_fwhm=interp_rs(rs,left_fwhm_idx) right=peak while right right_min+right_fwhm_threshold: right += 1 right_fwhm_idx = interpolate_spec(right, spec[right], right - 1, spec[right - 1],right_min+ right_fwhm_threshold) #if left > 0 else left right_fwhm=interp_rs(rs,right_fwhm_idx) right_fwhms.append(right_fwhm) left_fwhms.append(left_fwhm) real_cen = left_fwhm+(right_fwhm-left_fwhm)/2 real_cens.append(real_cen) return real_cens, left_fwhms, right_fwhms def find_interpolated_points(pos_pure, h_pure, rs, spec): # Sicherstellen, dass die Arrays sortiert sind sorted_indices = np.argsort(rs) rs = np.array(rs)[sorted_indices] spec = np.array(spec)[sorted_indices] # Finden der benachbarten Punkte um pos_pure idx_right = np.searchsorted(rs, pos_pure, side='right') idx_left = idx_right - 1 if idx_left < 0 or idx_right >= len(rs): raise ValueError("pos_pure liegt außerhalb des Bereichs von rs") # Suche nach Punkten links von pos_pure, die h_pure umschließen for i in range(idx_left, -1, -1): if spec[i]< h_pure: # if (spec[i] - h_pure) * (spec[i + 1] - h_pure) < 0: x1, y1 = rs[i], spec[i] x2, y2 = rs[i + 1], spec[i + 1] x_left = x1 + (h_pure - y1) * (x2 - x1) / (y2 - y1) break else: x_left = None # Suche nach Punkten rechts von pos_pure, die h_pure umschließen for i in range(idx_right, len(rs) - 1): if spec[i]< h_pure: #if (spec[i] - h_pure) * (spec[i + 1] - h_pure) < 0: x1, y1 = rs[i], spec[i] x2, y2 = rs[i - 1], spec[i - 1] x_right = x1 + (h_pure - y1) * (x2 - x1) / (y2 - y1) break else: x_right = None if x_right==None or x_left==None or np.abs(pos_pure-x_right)>100 or np.abs(pos_pure-x_left)>100: pos=np.nan else: pos=x_left+(x_right-x_left)/2 return x_left, x_right,pos def find_zeros_and_sums(arr,spectra, err_lim=0,spec_lim=0,dist=0): zero_indices = [] sums_between_zeros = [] max_indices = [] # Suche nach Vorzeichenwechseln (Nullstellen) for i in range(len(arr) - 1): if arr[i] * arr[i + 1] < 0: # Vorzeichenwechsel zero_indices.append(i) # Berechnung der Summen zwischen den Nullstellen for i in range(len(zero_indices) - 1): start, end = zero_indices[i], zero_indices[i + 1] segment = arr[start+1:end+1] # Werte zwischen den Nullstellen sum_segment = np.sum(np.abs(segment)) # Betrag der Summe # Nur speichern, wenn Summe den Grenzwert überschreitet if sum_segment > err_lim: local_max_index = start + 1 + np.argmax(segment) # Index des lokalen Maximums if spectra[local_max_index]>spec_lim: max_indices.append(local_max_index) sums_between_zeros.append(sum_segment) if max_indices: if len(max_indices)>1: idx_keep=[i for i in range(len(max_indices)-1) if (max_indices[i+1]-max_indices[i])>dist] idx_keep.append(len(max_indices)-1) filtered_sums = [sums_between_zeros[i] for i in idx_keep] filtered_max_indices = [max_indices[i] for i in idx_keep] else: filtered_sums = sums_between_zeros filtered_max_indices = max_indices else: filtered_sums=[] filtered_max_indices=[] return filtered_sums,filtered_max_indices def find_voigtprofile_pure(spectra,rs): """ finds pure compound voigt profiles Parameters ---------- spectra : array pure compound spectrum rs : array raman shift/cm^-1 Returns ------- paras_sorted : numpy array shape n_peakx4 pseudo voigt paras: lorentz racton , central peak position, peak with, peak area TYPE DESCRIPTION. """ deriv_1=deriv(spectra,rs) deriv_2=deriv(deriv_1,rs) deriv_2_rev=-deriv_2 all_peaks,_=find_peaks(deriv_2_rev,prominence=6e-07,height=1e-8)#all_peaks,_=find_peaks(deriv_2_rev,prominence=6e-07,height=1e-8) 'drop peaks with very small intensitites (also noise)' peaks_uncorr=np.array(all_peaks[spectra[all_peaks]>0.0004])#0.00035 ' OH bond' idx_oh_start=int((np.abs(rs - OH_start)).argmin()) if spectra[np.argmin(np.abs(rs-3302))]>0.0005: non_oh_peaks=peaks_uncorr[np.where(peaks_uncorr0.0005: oh_peaks.append(np.argmin(np.abs(rs-3115))) if spectra[np.argmin(np.abs(rs-3400))]>0.0005: oh_peaks.append(np.argmin(np.abs(rs-3400))) if spectra[np.argmin(np.abs(rs-3600))]>0.0002: oh_peaks.append(np.argmin(np.abs(rs-3600))) peaks=np.concatenate((non_oh_peaks,oh_peaks),axis=0) else: peaks=peaks_uncorr rs_int,left_fwhms, right_fwhms=find_real_cen(rs,deriv_2_rev,peaks) 'if OH bond is not found, find at least one peak from non-derivaed spectrum' nr_peaks=len(peaks) if nr_peaks==0: peaks,properties=find_peaks(spectra,prominence=5e-08,height=1e-4) nr_peaks=len(peaks) bounds_reshaped,pars_0=start_paras_pure(peaks,spectra,rs,200) pars_0_lq=pars_0.flatten() rs_arr=np.array([rs,]*nr_peaks).T res=least_squares(funtomin,pars_0_lq,bounds=bounds_reshaped.T,args=(spectra,rs_arr,nr_peaks),max_nfev=5000)#jac=an_jac, paras_new=res.x add_profs=True j=0 while add_profs==True and j<3: "add profiles" voigt_prof=pseudo_voigt(rs_arr,paras_new,nr_peaks) error=spectra-voigt_prof error_sums, index=find_zeros_and_sums(error,spectra,err_lim=0.003,spec_lim=0.0003) if error_sums: for i in range(len(error_sums)): if np.min(np.abs(peaks-index[i]))>5: peaks=np.append(peaks,index[i]) rs_int.append(rs[index[i]]) nr_peaks=len(peaks) bounds_reshaped,pars_0=start_paras_pure(peaks,spectra,rs,200) pars_0_lq=pars_0.flatten() rs_arr=np.array([rs,]*nr_peaks).T res=least_squares(funtomin,pars_0_lq,bounds=bounds_reshaped.T,args=(spectra,rs_arr,nr_peaks),max_nfev=5000)#jac=an_jac, paras_new=res.x j+=1 else: add_profs=False voigt_prof=pseudo_voigt(rs_arr,paras_new,nr_peaks) error=spectra-voigt_prof paras=np.reshape(paras_new,(nr_peaks,4)) paras_sorted=paras[paras[:,1].argsort()] plt.figure() for i in range(nr_peaks): singel_prof=pseudo_voigt(rs,paras[i,:],1) plt.plot(rs,singel_prof,color='b') rs_arr=np.array([rs,]*nr_peaks).T plt.plot(rs,pseudo_voigt(rs_arr,paras,nr_peaks)) plt.plot(rs,spectra) plt.plot(rs,error) plt.scatter(rs[index],error[index],color='r') plt.show() return paras_sorted,res.cost def mix_bounds_start_pc(paras_0,nr_peaks_mix,peakshift,peakshift_OH): """ get bounds of peakcenters for mixture fit Parameters ---------- paras_0 : float array shape(nr.peaks,2) initial parameters; voigt paras from last fit resp. pure compound paras for first mixture fit nr_peaks_mix : int total number of voigt profiles in compound 1 and 2 peakshift : float maximal peakshift in cm^-1 peakshift_OH : float maximal peakshift in cm^-1 for OH bond Returns ------- bounds_stacked : float array shape(nr.peaks,2) DESCRIPTION. """ 'bounds peak center' OH_exist=np.any(paras_0[:,1]>OH_start) if OH_exist == True: idx_OH_start=np.argmin(np.abs(paras_0[:,1]-OH_start)) lb_peak_center=np.concatenate([paras_0[:idx_OH_start,1]-peakshift,paras_0[idx_OH_start:,1]-peakshift_OH]) ub_peak_center=np.concatenate([paras_0[:idx_OH_start,1]+peakshift,paras_0[idx_OH_start:,1]+peakshift_OH]) else: lb_peak_center=paras_0[:,1]-peakshift ub_peak_center=paras_0[:,1]+peakshift bounds_stacked=np.column_stack([lb_peak_center,ub_peak_center]) 'avoid that lower and upper bound are equal in case peakshift is 0' bounds_stacked[bounds_stacked[:,0]>=bounds_stacked[:,1],1]=bounds_stacked[bounds_stacked[:,0]>=bounds_stacked[:,1],1]+1 return bounds_stacked def mix_bounds_start_wi(paras_0,nr_peaks_mix,widthchange,widthchange_OH): """ get bounds of peakwidth for mixture fit Parameters ---------- paras_0 : float array shape(nr.peaks,2) initial parameters; voigt paras from last fit resp. pure compound paras for first mixture fit nr_peaks_mix : int total number of voigt profiles in compound 1 and 2 widthchange : float maximal percentage of allowed widthchange widthchange_OH : float maximal percentage of allowed widthchange of OH bond Returns ------- bounds_stacked : TYPE DESCRIPTION. """ 'bounds width' max_peak_width=150 max_peak_width_OH=450 OH_exist=np.any(paras_0[:,1]>OH_start) if OH_exist == True: idx_OH_start=np.where(paras_0[:,1]>OH_start)[0][0] lb_peak_width=np.concatenate([paras_0[:idx_OH_start,2]-widthchange*paras_0[:idx_OH_start,2],paras_0[idx_OH_start:,2]-widthchange_OH*paras_0[idx_OH_start:,2]]) ub_not_OH=paras_0[:idx_OH_start,2]+widthchange*paras_0[:idx_OH_start,2] ub_not_OH[ub_not_OH>max_peak_width]=max_peak_width ub_OH=paras_0[idx_OH_start:,2]+widthchange_OH*paras_0[idx_OH_start:,2] ub_OH[ub_OH>max_peak_width]=max_peak_width_OH ub_peak_width=np.concatenate([ub_not_OH,ub_OH]) else: lb_peak_width=paras_0[:,2]-widthchange*paras_0[:,2] ub_peak_width=paras_0[:,2]+widthchange*paras_0[:,2] ub_peak_width[ub_peak_width>max_peak_width]=max_peak_width bounds_stacked=np.column_stack([lb_peak_width,ub_peak_width]) 'avoid that lower and upper bound are equal in case peakshift is 0' bounds_stacked[bounds_stacked[:,0]>=bounds_stacked[:,1],0]=bounds_stacked[bounds_stacked[:,0]>=bounds_stacked[:,1],1]-15 bounds_stacked[bounds_stacked<0]=0 return bounds_stacked def funtomin_pc_wt(paras,paras_tot,spec,rs_arr,idx_fit,idx_non_fit,error_min): """ objective function to optimize the paramters of the profiles which are fitted modeled spectrum is sum of fitted profiles + fixed profiles Parameters ---------- paras : float array (len=2*nr. profiles fitted) array containing the center and width of the profiles whcih shall be fitted paras_tot : float array (4* len(profiles)) all parameters of all voigt profiles spec : float array mixture spectrum rs_arr : float array (len(rs)*nr.profiles) raman shift array idx_fit : int list (len: nr. of fitted profiles) list containing the indices of the profiles which shall be fitted idx_non_fit : int list (len: nr. of fixed profiles) list containing the indices of the profiles which are fixed Returns ------- error: float array len(rs) difference between measured and fitted spectrum """ nr_peaks_fit=len(idx_fit) paras=np.reshape(paras,(nr_peaks_fit,2)) x=rs_arr[:,0:nr_peaks_fit]-paras[:,0] x_fix=rs_arr[:,nr_peaks_fit:]-paras_tot[idx_non_fit,1] v_profile_fit=np.sum(paras_tot[idx_fit,0]*paras_tot[idx_fit,3]*2/paras[:,1]*1/np.pi*1/((x/(paras[:,1]*0.5))**2+1)+(1-paras_tot[idx_fit,0])*paras_tot[idx_fit,3]*2/paras[:,1]*np.sqrt(np.log(2)/np.pi)*np.exp(-np.log(2)*(x/(0.5*paras[:,1]))**2),axis=1) v_profile_fix= np.sum(paras_tot[idx_non_fit,0]*paras_tot[idx_non_fit,3]*2/paras_tot[idx_non_fit,2]*1/np.pi*1/((x_fix/(paras_tot[idx_non_fit,2]*0.5))**2+1)+(1-paras_tot[idx_non_fit,0])*paras_tot[idx_non_fit,3]*2/paras_tot[idx_non_fit,2]*np.sqrt(np.log(2)/np.pi)*np.exp(-np.log(2)*(x_fix/(0.5*paras_tot[idx_non_fit,2]))**2),axis=1) error=spec-v_profile_fit-v_profile_fix return error def funtomin_pc(paras,paras_tot,spec,rs_arr,idx_fit,idx_non_fit): """ objective function to optimize the paramters of the profiles which are fitted modeled spectrum is sum of fitted profiles + fixed profiles Parameters ---------- paras : float array (len=2*nr. profiles fitted) array containing the center and width of the profiles whcih shall be fitted paras_tot : float array (4* len(profiles)) all parameters of all voigt profiles spec : float array mixture spectrum rs_arr : float array (len(rs)*nr.profiles) raman shift array idx_fit : int list (len: nr. of fitted profiles) list containing the indices of the profiles which shall be fitted idx_non_fit : int list (len: nr. of fixed profiles) list containing the indices of the profiles which are fixed Returns ------- error: float array len(rs) difference between measured and fitted spectrum """ nr_peaks_fit=len(idx_fit) x=rs_arr[:,0:nr_peaks_fit]-paras x_fix=rs_arr[:,nr_peaks_fit:]-paras_tot[idx_non_fit,1] v_profile_fit=np.sum(paras_tot[idx_fit,0]*paras_tot[idx_fit,3]*2/paras_tot[idx_fit,2]*1/np.pi*1/((x/(paras_tot[idx_fit,2]*0.5))**2+1)+(1-paras_tot[idx_fit,0])*paras_tot[idx_fit,3]*2/paras_tot[idx_fit,2]*np.sqrt(np.log(2)/np.pi)*np.exp(-np.log(2)*(x/(0.5*paras_tot[idx_fit,2]))**2),axis=1) v_profile_fix= np.sum(paras_tot[idx_non_fit,0]*paras_tot[idx_non_fit,3]*2/paras_tot[idx_non_fit,2]*1/np.pi*1/((x_fix/(paras_tot[idx_non_fit,2]*0.5))**2+1)+(1-paras_tot[idx_non_fit,0])*paras_tot[idx_non_fit,3]*2/paras_tot[idx_non_fit,2]*np.sqrt(np.log(2)/np.pi)*np.exp(-np.log(2)*(x_fix/(0.5*paras_tot[idx_non_fit,2]))**2),axis=1) error=spec-v_profile_fit-v_profile_fix return error def funtomin_wt(paras,paras_tot,spec,rs_arr,idx_fit,idx_non_fit): """ objective function to optimize the paramters of the profiles which are fitted modeled spectrum is sum of fitted profiles + fixed profiles Parameters ---------- paras : float array (len=2*nr. profiles fitted) array containing the center and width of the profiles whcih shall be fitted paras_tot : float array (4* len(profiles)) all parameters of all voigt profiles spec : float array mixture spectrum rs_arr : float array (len(rs)*nr.profiles) raman shift array idx_fit : int list (len: nr. of fitted profiles) list containing the indices of the profiles which shall be fitted idx_non_fit : int list (len: nr. of fixed profiles) list containing the indices of the profiles which are fixed Returns ------- error: float array len(rs) difference between measured and fitted spectrum """ nr_peaks_fit=len(idx_fit) x=rs_arr[:,0:nr_peaks_fit]-paras_tot[idx_fit,1] x_fix=rs_arr[:,nr_peaks_fit:]-paras_tot[idx_non_fit,1] v_profile_fit=np.sum(paras_tot[idx_fit,0]*paras_tot[idx_fit,3]*2/paras_tot[idx_fit,1]*1/np.pi*1/((x/(paras_tot[idx_fit,1]*0.5))**2+1)+(1-paras_tot[idx_fit,0])*paras_tot[idx_fit,3]*2/paras_tot[idx_fit,1]*np.sqrt(np.log(2)/np.pi)*np.exp(-np.log(2)*(x/(0.5*paras_tot[idx_fit,1]))**2),axis=1) v_profile_fix= np.sum(paras_tot[idx_non_fit,0]*paras_tot[idx_non_fit,3]*2/paras_tot[idx_non_fit,2]*1/np.pi*1/((x_fix/(paras_tot[idx_non_fit,2]*0.5))**2+1)+(1-paras_tot[idx_non_fit,0])*paras_tot[idx_non_fit,3]*2/paras_tot[idx_non_fit,2]*np.sqrt(np.log(2)/np.pi)*np.exp(-np.log(2)*(x_fix/(0.5*paras_tot[idx_non_fit,2]))**2),axis=1) error=spec-v_profile_fit-v_profile_fix return error def an_jac_mix(paras, paras_tot,intensity,rs_arr,idx_fit,idx_non_fit,error_min): """ calculation of analytical jacobian matrix in order to improve the speed of the calculation includes the partial derivatives of the partial voigt profiles according to center and with for all profiles whcih are fitted; fixed parameters are not relevant here Parameters ---------- paras : float array (len=2*nr. profiles fitted) array containing the center and width of the profiles whcih shall be fitted paras_tot : float array (4* len(profiles)) (not relevant) intensity : float array (not relevant) rs_arr : float array (len(rs)*nr.profiles) raman shift array idx_fit : int list (len: nr. of fitted profiles) list containing the indices of the profiles which shall be fitted idx_non_fit : int list (len: nr. of fixed profiles) (not relevant) Returns ------- jacobian: float array shape (len(rs)2*nr. of fitted peaks*) """ nr_peaks_fit=len(idx_fit) paras=np.reshape(paras,(nr_peaks_fit,2)) rs=rs_arr[:,0] jac=np.zeros((nr_peaks_fit*2,len(rs))) a=np.sqrt(np.log(2)/np.pi) for i in range (nr_peaks_fit): x=np.array(rs-paras[i,0]) exp_term=np.exp(-np.log(2)*(x/(0.5*paras[i,1]))**2) jac[i*2,:]=-(paras_tot[i,0]*paras_tot[i,3]*-16*paras[i,1]*x/(np.pi*(4*x**2+paras[i,1]**2)**2)+(1-paras_tot[i,0])*paras_tot[i,3]*-16*np.log(2)**1.5*x/(np.sqrt(np.pi)*paras[i,1]**3)*exp_term)# del peakcenter jac[i*2+1,:]=-2*paras_tot[i,3]*paras_tot[i,0]*(paras[i,1]**2-4*x**2)/(np.pi*(paras[i,1]**2+4*x**2)**2)+2*a*paras_tot[i,3]*(paras_tot[i,0]-1)*(paras[i,1]**2-8*np.log(2)*x**2)*exp_term/paras[i,1]**4 # del peakwidth jacobian=-jac.T return jacobian def find_idx_fit(voigt_paras_pure,peaks_non_fit): """ function to get indices of profiles which are fitted within the mixture Parameters ---------- voigt_paras_pure : numpy array shape n_peak x 4 pseudo voigt parameters peaks_non_fit : numpy array approximate peak centers of peaks which are fixed Returns ------- idx_fit : numpy array index of profiles that are fitted within voigt_paras_pure """ idx_non_fit=np.zeros_like(peaks_non_fit) for i in range(len(peaks_non_fit)): idx_non_fit[i]=np.argmin(np.abs(voigt_paras_pure-peaks_non_fit[i])) idx_fit = np.where(~np.isin(np.arange(len(voigt_paras_pure)), idx_non_fit))[0] return idx_fit def voigt_paras(spectra,rs,x1,peaks_non_fit_c1,peaks_non_fit_c2): """ Voigtfit of mixture molar spectra: lorentz percentage: is set peakcenter: is fit in first step peakwidth: is fit in second step intensity: is set by known mole fraction Parameters ---------- spectra : DataFrame mixture Raman spectra at different compositions rs : float array ramanshift x1 : Series mol fraction compound 1 pure_voigt_data_path : string data path to folder where pure compound parameters are saved compound_1 : string name of compound 1 compound_2 : string name of compound 1 Returns ------- DataFrame containing: 'paras_c1': fitted voigt parameters compound 1 (array len(x1)x 4) 'paras_c2':fitted voigt parameters compound 2 (array len(x1)x 4) 'voigt_profil_c1':fitted voigt profiles compound 1(array len(x1)x len(rs)) 'voigt_profil_c2':fitted voigt profiles compound 1(array len(x1)x len(rs)) 'voigt_fit':fitted voigt profiles compound 1 an 2 (=> mixture voigt fit) (array len(x1)x len(rs)) 'voigt_ex_c1':excess voigt profiles compound 1(array len(x1)x len(rs)) 'voigt_ex_c2':excess voigt profiles compound 2(array len(x1)x len(rs)) """ spectra=np.vstack(np.array(spectra)).astype(np.float64) x1=(np.vstack(np.array(x1)).astype(np.float64)).flatten() """ Check if pure profile parameters exist; if yes load if not fit them""" voigt_paras_c2_pure,error_pure_c2=find_voigtprofile_pure(spectra[0,:],rs) voigt_paras_c1_pure,error_pure_c1=find_voigtprofile_pure(spectra[-1,:],rs) idx_fit_c1=find_idx_fit(voigt_paras_c1_pure[:,1],peaks_non_fit_c1) idx_fit_c2=find_idx_fit(voigt_paras_c2_pure[:,1],peaks_non_fit_c2) """ prepare the fit of the mixture spectra; """ nr_peaks_c1=len(voigt_paras_c1_pure[:,1]) nr_peaks_c2=len(voigt_paras_c2_pure[:,1]) nr_peaks_mix=nr_peaks_c1+nr_peaks_c2 nr_points=len(x1) rs_arr=np.array([rs,]*nr_peaks_mix).T rs_arr_c1=np.array([rs,]*nr_peaks_c1).T rs_arr_c2=np.array([rs,]*nr_peaks_c2).T all_voigt_paras_c1=np.empty((nr_points,nr_peaks_c1,4)) all_voigt_paras_c2=np.empty((nr_points,nr_peaks_c2,4)) error=np.zeros_like(x1) idx_pure_c1=np.where(x1==1) idx_pure_c2=np.where(x1==0) all_voigt_paras_c1[idx_pure_c2,:,:]=np.zeros_like(voigt_paras_c1_pure) all_voigt_paras_c2[idx_pure_c1,:,:]=np.zeros_like(voigt_paras_c2_pure) all_voigt_paras_c1[:]=np.nan all_voigt_paras_c2[:]=np.nan all_voigt_paras_c1[idx_pure_c1,:,:]=voigt_paras_c1_pure all_voigt_paras_c2[idx_pure_c2,:,:]=voigt_paras_c2_pure voigt_paras_c1=voigt_paras_c1_pure voigt_paras_c2=voigt_paras_c2_pure int_c1=voigt_paras_c1_pure[:,3] int_c2=voigt_paras_c2_pure[:,3] """ The profiles with a width of max_peakwidth shall be excluded from the fit. Therfore, the corresponding indices have to be found. """ all_pars=np.concatenate((voigt_paras_c1_pure,voigt_paras_c2_pure),axis=0) idx_fit = np.concatenate((idx_fit_c1,nr_peaks_c1+ idx_fit_c2)) idx_non_fit=np.where(~np.isin(np.arange(nr_peaks_mix), idx_fit))[0] ''' Fit of parameters between x1=0.5-1 parameters of c1 are adjusted on the partial spectra of c1 prameters of c2 are adjusted on the mixture spectrum-fitted spectrum of c1 ''' "start at equimolar composition" idx_equimol=np.abs(x1-0.5).argmin() first=True for i in range(idx_equimol,nr_points,1): if x1[i] == x1[i-1]: all_voigt_paras_c1[i,:,:]=all_voigt_paras_c1[i-1,:,:] all_voigt_paras_c2[i,:,:]=all_voigt_paras_c2[i-1,:,:] elif x1[i]==1 or x1[i]==0: pass else: peakshift_c1=np.abs(x1[i]-x1[i-1])*150 #15 peakshift_c2=np.abs(x1[i]-x1[i-1])*150 #50 peakshift_OH=40 widthchange_c1=np.abs(x1[i]-x1[i-1])*10 widthchange_c2=np.abs(x1[i]-x1[i-1])*10 widthchange_OH=np.abs(x1[i]-x1[i-1])*100 if i==idx_equimol: peakshift_c1=20#4 peakshift_c2=20#4 peakshift_OH=125 widthchange_OH=0.8 #widthchange_c1=0.1 #widthchange_c2=0.1 widthchange_c1=0.2 widthchange_c2=0.2 voigt_fit=[] "Get initial parameters and bounds of centers and widths for the fit" bounds_reshaped_pc_c1=mix_bounds_start_pc(voigt_paras_c1,nr_peaks_c1,peakshift_c1,peakshift_OH) bounds_reshaped_wi_c1=mix_bounds_start_wi(voigt_paras_c1,nr_peaks_c1,widthchange_c1,widthchange_OH) bounds_reshaped_pc_c2=mix_bounds_start_pc(voigt_paras_c2,nr_peaks_c2,peakshift_c2,peakshift_OH) bounds_reshaped_wi_c2=mix_bounds_start_wi(voigt_paras_c2,nr_peaks_c2,widthchange_c2,widthchange_OH) bounds_reshaped_pc=np.concatenate([bounds_reshaped_pc_c1, bounds_reshaped_pc_c2]) bounds_reshaped_width=np.concatenate([ bounds_reshaped_wi_c1, bounds_reshaped_wi_c2]) bounds_stacked=np.column_stack([bounds_reshaped_pc,bounds_reshaped_width]) pars_0=np.concatenate([voigt_paras_c1[:,1:3],voigt_paras_c2[:,1:3]]) "exclude peaks that are fixed" bounds_stacked_fit=bounds_stacked[idx_fit,:] bounds_reshaped=np.reshape(bounds_stacked_fit,(2*len(idx_fit),2)) pars_0_fit=pars_0[idx_fit,:].flatten() mask = pars_0_fit > bounds_reshaped[:, 1] pars_0_fit[mask] = bounds_reshaped[:, 1][mask] * 0.9 "Prepare fixed parameters" paras_tot_c1= np.column_stack([voigt_paras_c1[:,0:3],x1[i]*voigt_paras_c1_pure[:,3]]) paras_tot_c2= np.column_stack([voigt_paras_c2[:,0:3],(1-x1[i])*voigt_paras_c2_pure[:,3]]) paras_tot=np.append(paras_tot_c1,paras_tot_c2,axis=0) "fit" error_min=x1[i]*error_pure_c1+(1-x1[i])*error_pure_c2 res=least_squares(funtomin_pc_wt,pars_0_fit,jac=an_jac_mix,bounds=bounds_reshaped.T,args=(paras_tot,spectra[i,:],rs_arr,idx_fit,idx_non_fit,error_min),ftol=1e-12,gtol=1e-12,xtol=1e-12)# error[i] =res.cost "reshape fitted paras and store them " all_paras_fitted=np.reshape(res.x,(len(idx_fit),2)) all_pars[idx_fit,1:3]=all_paras_fitted voigt_paras_c1=np.column_stack([all_pars[:nr_peaks_c1,0:3],x1[i]*int_c1]) voigt_fit_c1=pseudo_voigt(rs_arr_c1,voigt_paras_c1,nr_peaks_c1) voigt_paras_c2=np.column_stack([all_pars[nr_peaks_c1:,0:3],(1-x1[i])*int_c2]) voigt_fit_c2=pseudo_voigt(rs_arr_c2,voigt_paras_c2,nr_peaks_c2) if i==idx_equimol: pars1_equ=voigt_paras_c1 pars2_equ=voigt_paras_c2 voigt_fit=voigt_fit_c1+voigt_fit_c2 all_voigt_paras_c1[i,:,:]= voigt_paras_c1 all_voigt_paras_c2[i,:,:]= voigt_paras_c2 if first==True: #i=i-1 first=False ''' Fit of parameters between x1=0.5-0 parameters of c2 are adjusted on the partial spectra of c2 prameters of c1 are adjusted on the mixture spectrum-fitted spectrum of c2 ''' voigt_paras_c1=pars1_equ voigt_paras_c2=pars2_equ for i in range(idx_equimol-1,0,-1): if x1[i] == x1[i+1]: all_voigt_paras_c1[i,:,:]=all_voigt_paras_c1[i+1,:,:] all_voigt_paras_c2[i,:,:]=all_voigt_paras_c2[i+1,:,:] elif x1[i]==1 or x1[i]==0: pass else: peakshift_c1=np.abs(x1[i]-x1[i-1])*150 #15 peakshift_c2=np.abs(x1[i]-x1[i-1])*150 #50 peakshift_OH=40 widthchange_c1=np.abs(x1[i]-x1[i-1])*10 widthchange_c2=np.abs(x1[i]-x1[i-1])*10 widthchange_OH=np.abs(x1[i]-x1[i-1])*100 "Get initial parameters and bounds of centers and widths for the fit" bounds_reshaped_pc_c1=mix_bounds_start_pc(voigt_paras_c1,nr_peaks_c1,peakshift_c1,peakshift_OH) bounds_reshaped_wi_c1=mix_bounds_start_wi(voigt_paras_c1,nr_peaks_c1,widthchange_c1,widthchange_OH) bounds_reshaped_pc_c2=mix_bounds_start_pc(voigt_paras_c2,nr_peaks_c2,peakshift_c2,peakshift_OH) bounds_reshaped_wi_c2=mix_bounds_start_wi(voigt_paras_c2,nr_peaks_c2,widthchange_c2,widthchange_OH) bounds_reshaped_pc=np.concatenate([bounds_reshaped_pc_c1, bounds_reshaped_pc_c2]) bounds_reshaped_width=np.concatenate([ bounds_reshaped_wi_c1, bounds_reshaped_wi_c2]) bounds_stacked=np.column_stack([bounds_reshaped_pc,bounds_reshaped_width]) pars_0=np.concatenate([voigt_paras_c1[:,1:3],voigt_paras_c2[:,1:3]]) "exclude peaks that are fixed" bounds_stacked_fit=bounds_stacked[idx_fit,:] bounds_reshaped=np.reshape(bounds_stacked_fit,(2*len(idx_fit),2)) pars_0_fit=pars_0[idx_fit,:].flatten() mask = pars_0_fit > bounds_reshaped[:, 1] pars_0_fit[mask] = bounds_reshaped[:, 1][mask] * 0.9 "Prepare fixed parameters" paras_tot_c1= np.column_stack([voigt_paras_c1[:,0:3],x1[i]*voigt_paras_c1_pure[:,3]]) paras_tot_c2= np.column_stack([voigt_paras_c2[:,0:3],(1-x1[i])*voigt_paras_c2_pure[:,3]]) paras_tot=np.append(paras_tot_c1,paras_tot_c2,axis=0) "fit" error_min=x1[i]*error_pure_c1+(1-x1[i])*error_pure_c2 res=least_squares(funtomin_pc_wt,pars_0_fit,jac=an_jac_mix,bounds=bounds_reshaped.T,args=(paras_tot,spectra[i,:],rs_arr,idx_fit,idx_non_fit,error_min),ftol=1e-12,gtol=1e-12)# error[i] =res.cost "reshape fitted paras and store them " all_paras_fitted=np.reshape(res.x,(len(idx_fit),2)) all_pars[idx_fit,1:3]=all_paras_fitted voigt_paras_c1=np.column_stack([all_pars[:nr_peaks_c1,0:3],x1[i]*int_c1]) voigt_fit_c1=pseudo_voigt(rs_arr_c1,voigt_paras_c1,nr_peaks_c1) voigt_paras_c2=np.column_stack([all_pars[nr_peaks_c1:,0:3],(1-x1[i])*int_c2]) voigt_fit_c2=pseudo_voigt(rs_arr_c2,voigt_paras_c2,nr_peaks_c2) all_voigt_paras_c1[i,:,:]= voigt_paras_c1 all_voigt_paras_c2[i,:,:]= voigt_paras_c2 voigt_fit_c1=np.empty((nr_points,len(rs))) voigt_fit_c2=np.empty_like(voigt_fit_c1) voigt_fit=np.empty_like(voigt_fit_c2) for i in range (nr_points): if x1[i]==0: voigt_fit_c1[i,:]=np.zeros((len(rs),1)).flatten() voigt_fit_c2[i,:]=pseudo_voigt(rs_arr_c2,(all_voigt_paras_c2[i,:,:]).flatten(),nr_peaks_c2) elif x1[i]==1: voigt_fit_c2[i,:]=np.zeros((len(rs),1)).flatten() voigt_fit_c1[i,:]=pseudo_voigt(rs_arr_c1,all_voigt_paras_c1[i,:,:],nr_peaks_c1) else: voigt_fit_c1[i,:]=pseudo_voigt(rs_arr_c1,all_voigt_paras_c1[i,:,:],nr_peaks_c1) voigt_fit_c2[i,:]=pseudo_voigt(rs_arr_c2,all_voigt_paras_c2[i,:,:],nr_peaks_c2) voigt_fit[i,:]=voigt_fit_c1[i,:]+voigt_fit_c2[i,:] list_c1=[] for i in range(nr_points): list_c1.append([]) list_c1[i].append(all_voigt_paras_c1[i,::]) list_c2=[] for i in range(nr_points): list_c2.append([]) list_c2[i].append(all_voigt_paras_c2[i,::]) fitted_True_c1 = ['True' if i in idx_fit_c1 else 'False' for i in range(len(voigt_paras_c1_pure[:, 0]))] fitted_True_c2 = ['True' if i in idx_fit_c2 else 'False' for i in range(len(voigt_paras_c2_pure[:, 0]))] frame={'paras_c1':pd.Series(list_c1), 'paras_c2':pd.Series(list_c2), 'fitted_c1':pd.Series(fitted_True_c1), 'fitted_c2':pd.Series(fitted_True_c2), 'voigt_profil_c1':pd.Series(list(voigt_fit_c1)), 'voigt_profil_c2':pd.Series(list(voigt_fit_c2)), 'voigt_fit':pd.Series(list(voigt_fit)), 'error':pd.Series(error) } res=pd.DataFrame(frame) return(res)