<jupyter_start><jupyter_text>Peak Decoupling through Fitting This dataset contains a particle size distribution curve: columns of the volume % and the particle size. In this exercise, I will resolve the overlapping peaks by fitting a sum of two lognormal distributions. By resolving, the mean peak positions of each underlying moiety can be identified.  Import modules & define formatting<jupyter_code>import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import numpy as np
import math
from matplotlib import rcParams
import pandas as pd
rcParams['font.family']='sans-serif'
rcParams['font.sans-serif']=['Arial']
size=10<jupyter_output><empty_output><jupyter_text>Define a class to read & get data; store the models for fitting<jupyter_code>class AnalyzePSD(object):
    def __init__(self,filename):
        '''input:  filename is of type string, it indicates the name of the file for reading
        Read the file and store it in self.datafile'''
        self.filename=filename
        self.datafile=pd.read_csv(filename,engine='python',sep='\t',header=2)

    def getDataFile(self):
        '''returns the pd.DataFrame of the file that was read in'''
        return self.datafile

    def getVolSize(self):
        '''return the size and volume density data as pd.Series'''
        size=pd.to_numeric(self.datafile['Size Classes (um)'],downcast='float',errors='coerce')
        vol=pd.to_numeric(self.datafile['Volume Density (%)'],downcast='float',errors='coerce')
        return size,vol

    def Pmod(self,x,sigma,mu,c):
        '''define a function to contain the log normal model'''
        denom1=x*sigma*(2*math.pi)**(1/2)
        expterm=((np.log(x)-mu)/sigma)**2
        Ptot=c*(np.exp(-.5*expterm)/denom1)
        return Ptot

    def PmodBimodal(self,x,sigma,mu,c,sigma2,mu2,d):
        '''define a function to contain a bimodal log normal model'''
        denom1=x*sigma*(2*math.pi)**(1/2)
        denom2=x*sigma2*(2*math.pi)**(1/2)
        expterm=((np.log(x)-mu)/sigma)**2
        expterm2=((np.log(x)-mu2)/sigma2)**2
        Ptot=c*(np.exp(-.5*expterm)/denom1)+d*(np.exp(-.5*expterm2)/denom2)
        return Ptot<jupyter_output><empty_output><jupyter_text>Define a Graph object for graphing and formatting<jupyter_code>class Graph(object):
    def __init__(self,fx,fy):
        '''input:  fx and fy are type integers--they are the dimensions of the figure'''
        self.fx=fx
        self.fy=fy
        fig=plt.figure(figsize=(fx,fy),dpi=300)
        
    def graphPSD(self,x,y,sym,clr,samp):
        '''x and y are of the same length (pd.Series type)
        sym is type string identifying the marker symbol
        clr is type string, denoting the color for the marker
        samp is type string, naming the name of the data series for the legend'''
        plt.semilogx(x,y,sym,color=clr,markersize=3,mec=clr,label=samp)
        plt.xlabel('Size ($\mu$m)',fontsize=size)
        plt.ylabel('Volume Density (%)', fontsize=size)
        return
    
    def formatGraph(self):
        '''format the Graph to have ticks on all sides
        incorporate a legend
        and ensure that everything fits inside the dimensions previously noted'''
        plt.tick_params(labelsize=size)
        plt.tick_params(which='major',right='on',direction='in',top='on',length=6)
        plt.tick_params(which='minor',right='on',direction='in',top='on',length=3)
        plt.legend(loc='best',fontsize=size-5)
        plt.tight_layout()
        return<jupyter_output><empty_output><jupyter_text>Call the file and read the data<jupyter_code>file='PSD_Sample1.txt'
dataSet=AnalyzePSD(file)    #instantiate an object of AnalyzePSD
psize,vol=dataSet.getVolSize()
psizefit=psize[34:59]           #select the region for fitting
volfit=vol[34:59]                   #select the region for fitting<jupyter_output><empty_output><jupyter_text>Run the fit using curve_fit from SciPy<jupyter_code>popt,pcov=curve_fit(dataSet.PmodBimodal,psizefit,volfit)<jupyter_output><empty_output><jupyter_text>Print the results to the fit<jupyter_code>print('The results from the fit to the log normal')
print('\nsigma=', popt[0],', mu=',popt[1], ', c=', popt[2],'sigma2=',popt[3],', mu2=',popt[4],', d=',popt[5])<jupyter_output>The results from the fit to the log normal

sigma= 0.3021398865896843 , mu= 0.7411879200328817 , c= 15.088861115070719 sigma2= 0.5245567029884016 , mu2= 1.7576367920222489 , d= 40.479904410888096<jupyter_text>Plot the data and the results from the model<jupyter_code>fig1=Graph(3,3)
fig1.graphPSD(psize,vol,'s','b',file[:-4])   # this is the whole curve
fig1.graphPSD(psizefit,volfit,'.','orchid','Pts for Fit')
#
fig1.graphPSD(psize,dataSet.Pmod(psize,popt[0],popt[1],popt[2]),'--','slategrey','Peak 1 Fit')
fig1.graphPSD(psize,dataSet.Pmod(psize,popt[3],popt[4],popt[5]),'--','orange','Peak 2 Fit')
fig1.graphPSD(psize,dataSet.PmodBimodal(psize,*popt),'-','r','Bimodal Log-Norm Fit')
plt.axis([.1,100,-1,12])
fig1.formatGraph()
plt.show()<jupyter_output><empty_output><jupyter_text>Determine the mean sizes of each PeakFrom the graph above, the peak was resolved into two peaks.  What is the mean particle size associated with each of those peaks?The following calculations will calculate the means for peak 1 and peak 2.  First, find the area under the curve for each resolved peak:<jupyter_code>area1=np.trapz(dataSet.Pmod(psize,popt[0],popt[1],popt[2]),psize)
area2=np.trapz(dataSet.Pmod(psize,popt[3],popt[4],popt[5]),psize)<jupyter_output><empty_output><jupyter_text>Next, determine the cumulative sum of the y values in the distribution.  That is, the cumulative sum of the volume density.<jupyter_code>sum1=dataSet.Pmod(psize,popt[0],popt[1],popt[2]).sum()
sum2=dataSet.Pmod(psize,popt[3],popt[4],popt[5]).sum()<jupyter_output><empty_output><jupyter_text>Normalize each function to the area such that the sum of the y-values equals 1<jupyter_code>norm1=(dataSet.Pmod(psize,popt[0],popt[1],popt[2])/sum1)
norm2=(dataSet.Pmod(psize,popt[3],popt[4],popt[5])/sum2)

normsum1=norm1.sum()
normsum2=norm2.sum()<jupyter_output><empty_output><jupyter_text>Now the mean of each peak can be determined:<jupyter_code>weighted1=(psize*norm1).sum()
weighted2=(psize*norm2).sum()

print('The mean sizes of each peak are the following: peak 1 =',round(weighted1,2), 'um and peak 2 = ',round(weighted2, 2), 'um.')<jupyter_output>The mean sizes of each peak are the following: peak 1 = 2.0 um and peak 2 =  5.05 um.<jupyter_text>Correlation Coefficient of Fit:  R-squaredFinally, how do we quantify how good the fit was to the data? This can be represented by the correlation coefficient.  The closer it is to a value, the better the fit. First, determine the residuals:  the difference between the ydata and the yvalues from the fit.<jupyter_code>residuals=volfit-dataSet.PmodBimodal(psize,*popt)<jupyter_output><empty_output><jupyter_text>Next, square the residuals, and compute its sum.<jupyter_code>residual_sum_squares=np.sum(residuals**2)<jupyter_output><empty_output><jupyter_text>Next, determine the sum of the squares of the data.<jupyter_code>total_sum_squares=np.sum((volfit-np.mean(volfit))**2)<jupyter_output><empty_output><jupyter_text>Finally, the correlation coefficient can be determined through the following.<jupyter_code>r_squared=1-(residual_sum_squares/total_sum_squares)
print('\nR-squared=',round(r_squared,3))<jupyter_output>R-squared= 0.996
{'fraction_non_alphanumeric': 0.13167163370360335, 'fraction_numerical': 0.032782443782172856, 'mean_word_length': 5.167919799498747, 'pattern_counts': {'":': 0, '<': 58, '<?xml version=': 0, '>': 58, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '627994', 'license': 'MIT', 'n_tokens_mistral': 2719, 'n_tokens_neox': 2560, 'n_words': 604, 'path': 'LogNormal_Fitting.ipynb', 'repo_name': 'aarojas20/PSD-Peak-Decoupling', 'size': '7382'}