% Time Resolved Classical Least Squares (TRCLS). Simulation of the
% separation of three similar components by chromatograpy with an array detector.
% You can change the concentrations of the unknown mixture in line 29, the
% random detector noise in line 18, the spectra in lines 33-48, and the
% chromatographic peaks in lines 62-84. To make it more realistic,
% simulated sample spectra and chromatographic peaks are based on the three
% acetophenone positional isomers: o-methyl (o-MAP), m-methyl (m-MAP), and
% p-methyl (p-MAP), estimated from the graphs in the technical report from
% Shimadazu: https://solutions.shimadzu.co.jp/an/n/en/hplc/jpl217011.pdf
% Note. This script requires the following downloadable functions in the
% Matlab path, all freely available at http://tinyurl.com/cey8rwh
% modelpeaks.m, expgaussian.m, autopeaks.m, halfwidth.m, val2ind.m
clear
clf
format compact
format short g

% -------Initial settings, to change as you wish (lines 18-55) --------
wincrement=1.5; % wavelength increment (= slit width in Shimadzu specs)
w=210:wincrement:320; % wavelength axis 
tau=4; % Exponential time constant (larger = more asymmetrical)
DetNoise=.002; % Detector noise standard deviation (Shimadzu specification)
SmoothWidth=1; % Smoothing of chromatographic peaks to reduce noise
%
% SOLUTION concentrations
% Known concentrations of components in the standard solution [1 3 2]
M=[400 400 400]; 
% Unknown concentrations of components in the simulated unknown sample [1 3 2] 
U=[.1 .2 .4]; 
Concentrations=U
%
% SPECTRAL characteristics of component 1 (based on OrthoMAP)
opeaks=[200 245 290]; % Peak wavelength, nm
omAU=[900 220 40];  % Peak absorbance, mAU
owidth=[20 25 30];  % Peak width, nm (FWHM)
% Spectral characteristics of component 2 (based on MetaMAP)
mpeaks=[200 250 290]; % Peak wavelength, nm
mmAU=[900 260 40]; % Peak absorbance, mAU
mwidth=[20 25 35];% Peak width, nm (FWHM)
% Spectral characteristics of component 3 (based on ParaMAP)
ppeaks=[200 255 290]; % Peak wavelength, nm
pmAU=[500 360 4]; % Peak absorbance, mAU
pwidth=[18 30 60];% Peak width, nm (FWHM)
%
% CHROMATOGRAPHY of Component 1
ort=5.1; % ortho retention time (min) 
ow=.07 ; % ortho peak width (min)
oh=220 ; % ortho peak height (mAU)
% Chromatograpy of Component 2
mrt=5.2; % meta retention time (min) 
mw=.07 ; % meta peak width (min)
mh=220 ; % meta peak height (mAU)
% Chromatograpy of Component 3
prt=5.3; % para retention time (min)
pw=.07 ; % para peak width (min)
ph=260 ; % para peak height (mAU)

% ------------ Calculations ---------------------------
% Each spectrum is modeled as the sum of 3 Gaussians + random detector noise
% Compute spectrum of component 1 (OrthoMAP)
oMAPspect=modelpeaks(w,3,1,omAU,opeaks,owidth,0,0)+DetNoise.*randn(1,length(w));
% Compute spectrum of component 2 (MetaMAP)
mMAPspect=modelpeaks(w,3,1,mmAU,mpeaks,mwidth,0,0)+DetNoise.*randn(1,length(w));
% Compute spectrum of component 3 (ParaMAP)
pMAPspect=modelpeaks(w,3,1,pmAU,ppeaks,pwidth,0,0)+DetNoise.*randn(1,length(w));
MixtureSpectrum=U(1)*oMAPspect+U(2)*mMAPspect+U(3)*pMAPspect+DetNoise.*randn(1,length(w));

figure(1)
clf
% Plot all three spectra in different colors
plot(w,oMAPspect,'b',w,mMAPspect,'g',w,pMAPspect,'r','linewidth',2)
xlabel('Wavelength, nm')
ylabel('Absorbance, mAU')
title ('Simulated spectra of Components 1 (blue), 2 (green), and 3 (red) at 400 ug/mL ')
axis([210 320 0 500])
grid

% ---Chromatography----
% Chromatographic peaks (modeled as Exponentially modified Gaussians)
% Peak of component 1
StartTime=ort-2.*ow;
EndTime=prt+2.*pw;
t=StartTime:.0035:EndTime; % Time axis (~1.5 minutes total duration)
oMAPchrom=oh*expgaussian(t,ort,ow,-tau); % tau is exponential time constant
oMAParea=trapz(t,oMAPchrom); % area by trapezoidal numerical integration
hwo=halfwidth(t,oMAPchrom);
rto=t(val2ind(oMAPchrom,max(oMAPchrom))); % Retention time
% Peak of component 2
mMAPchrom=mh*expgaussian(t,mrt,mw,-tau);% tau is exponential time constant
mMAParea=trapz(t,mMAPchrom); % area by trapezoidal numerical integration
hwm=halfwidth(t,mMAPchrom);
rtm=t(val2ind(mMAPchrom,max(mMAPchrom))); % Retention time
% Peak of component 3
pMAPchrom=ph*expgaussian(t,prt,pw,-tau);% tau is exponential time constant
pMAParea=trapz(t,pMAPchrom); % area by trapezoidal numerical integration
hwp=halfwidth(t,pMAPchrom);
rtp=t(val2ind(pMAPchrom,max(pMAPchrom))); % Retention time

TrueRetentionTimes=[rto rtm rtp];

% Plot all three chromatographic peaks in different colors
figure(2)
clf
chrom=oMAPchrom+mMAPchrom+pMAPchrom+DetNoise.*randn(1,length(t))'; % Sum of the three peaks
plot(t,oMAPchrom,'b',t,mMAPchrom,'g',t,pMAPchrom,'r',t,chrom,'k.','linewidth',2)
axis([min(t) max(t) -20 max(chrom).*1.1])
grid
ylabel('Absorbance, mAU')
xlabel('Time, minutes')
title('Simulated chromatographic peaks of separate components at 400 ug/mL')
maxy=1.05*max(chrom); 
text(ort,maxy,'1');text(mrt,maxy,'2');text(prt,maxy,'3')
A=[oMAPspect;mMAPspect;pMAPspect];
H=[oMAPchrom';mMAPchrom';pMAPchrom'];
E=A'./M; % Analytical sensitivity


figure(3)
clf
% Alternative calculation #1, lines 113-146, to animate the evolution
%  of the chromatographic profiles in the subplots of Figure 3.
C=zeros(3,length(t))';
for x=1:length(t)
    Do=(oMAPchrom(x)*oMAPspect).*U(1)./(M(1)); % oMAP at time x
    Dm=(mMAPchrom(x)*mMAPspect).*U(2)./(M(2)); % mMAP at time x
    Dp=(pMAPchrom(x)*pMAPspect).*U(3)./(M(3)); % pMAP at time x
    % Use unknown concentrations U to weight the
    % contribution of each component in the spectrum
    % of the sample and add DetNoise mAU random noise.
    % Spectral profile at the detector at time t:
    D=U(1)*Do + U(2)*Dm + U(3)*Dp + DetNoise.*randn(1,length(w)); 
    % Use analytical sensitivity E to compute
    % concentrations C of each component at time t(x)
    C(x,:)=D/E'; % C = concentrations of each component at time t(x)
    % Plot the chromatographic profile of each separate component
    CS=C*A; % The Chromatographic-Spectroscopic matrix, CS
    
    figure(3)
    clf
    plot(w,D,'k','linewidth',2)
    grid
    ylabel('Absorbance, mAU')
    xlabel('Wavelength, nm')
    title(['Spectral profile at the detector at time t = ' num2str(t(x),3) ' minutes' ])
    axis([min(w) max(w) -.01.*max(MixtureSpectrum) max(MixtureSpectrum)/3])  
    
    figure(4)
    clf
    CompDet(x)=rank(CS',100);
    plot(t(1:x),CompDet,'linewidth',2)
    % if x>3,contour(CS');end
    grid
    ylabel('Rank of H"*A matrix')
    xlabel('Time, minutes')
    title('Number of components detected in sample (by matrix rank)')

    figure(5)
    clf
    subplot(2,2,1)
    y1=fastsmooth(C(:,1),SmoothWidth,3);
    plot(t,y1,'b','linewidth',1.5)
    grid
    xlabel('Time, minutes')
    ylabel('Absorbance, mAU')
    title(['Component 1 = ' num2str(U(1)) 'ug/mL' ])
    
    subplot(2,2,2)
    y2=fastsmooth(C(:,2),SmoothWidth,3);
    plot(t,y2,'g','linewidth',1.5)
    grid
    ylabel('Absorbance, mAU')
    xlabel('Time, minutes')
    title(['Component 2 = ' num2str(U(2)) 'ug/mL' ])
    
    subplot(2,2,3)
    y3=fastsmooth(C(:,3),SmoothWidth,3);
    plot(t,y3,'r','linewidth',1.5)
    ylabel('Absorbance, mAU')
    xlabel('Time, minutes')
    title(['Component 3 = ' num2str(U(3)) 'ug/mL' ])
    grid
    
    subplot(2,2,4)
    yobsd=[oMAPchrom mMAPchrom pMAPchrom]*U'+DetNoise.*randn(1,length(t))';
    plot(t,y1+y2+y3,'k',t,y1,'b',t,y2,'g',t,y3,'r','linewidth',1.5)
    % plot(t(1:x),yobsd(1:x),'k','linewidth',1.5)
    % axis([min(t) max(t) -.1.*max(yobsd) 1.1.*max(yobsd)])
    grid
    ylabel('Absorbance, mAU')
    xlabel('Time, minutes')
    title('Chromatogram of sample mixture')

end

% Alternative calculation #2, lines 151-171, matrix math method
% rather than a loop: shorter and faster, but without animation
% C=(H'.*U)+DetNoise.*randn(1,length(t))';
% % Plot the chromatographic profile of each separate component
% figure(3)
% subplot(2,2,1)
% y1=fastsmooth(C(:,1),SmoothWidth,3);
% plot(t,y1,'b')
% ylabel('mAU')
% xlabel('Time, minutes')
% title('Component 1')
% subplot(2,2,2)
% y2=fastsmooth(C(:,2),SmoothWidth,3);
% plot(t,y2,'g')
% ylabel('mAU')
% xlabel('Time, minutes')
% title('Component 2')
% subplot(2,2,3)
% y3=fastsmooth(C(:,3),SmoothWidth,3);
% plot(t,y3,'r')
% ylabel('mAU')
% xlabel('Time, minutes')
% title('Component 3')
    

% The CS matrix contains both the chromatographic (C) and Spectrocopic (S)
% data in one matrix. To dispay a 3-D surface or a contour plot, type
% "mesh(CS)" or "contour(CS)".
% To plot the spectrum at time "time" (in minutes), type 
% "plot(w,CS(val2ind(t,time),:))".
% To plot the chromatogram at any single wavelength "nm", type
% "plot(t,CS(:,val2ind(w,nm)))".

% Report and plot results
ComponentsDetected=rank(H'*A)
HalfWidths=[halfwidth(t,y1) halfwidth(t,y2) halfwidth(t,y3)];
to1=t(val2ind(y1,max(y1))); % Retention time 1 measured by CLS
to2=t(val2ind(y2,max(y2))); % Retention time 2 measured by CLS
to3=t(val2ind(y3,max(y3))); % Retention time 3 measured by CLS
disp('     Component 1  Component 2  Component 3 ')
TrueRetentionTimes
CLSRetentionTimes=[to1 to2 to3]
disp(' ')
maxy=.9*max(y1+y2+y3);
text(ort,maxy,'1');text(mrt,maxy,'2');text(prt,maxy,'3')
TruePeakAreas=[oMAParea mMAParea pMAParea].*U
CLSPeakAreas=[trapz(t,y1) trapz(t,y2) trapz(t,y3)]./U
CLSPercentError=100.*[(CLSPeakAreas(1)-TruePeakAreas(1))./CLSPeakAreas(1) (CLSPeakAreas(2)-TruePeakAreas(2))./CLSPeakAreas(2) (CLSPeakAreas(3)-TruePeakAreas(3))./CLSPeakAreas(3)]
disp(' ')
% disp('Attempt to detect and measure peaks directly using autopeaks.m')
try
    P=autopeaks(t',chrom,2*ComponentsDetected);
    sizeP=size(P);
    if sizeP(1)==3
        disp('Perpendicular Drop Method, using PerpDropAreas.m')
        PerpDropRetentionTimes=[P(1,2) P(2,2) P(3,2)]
        PerpDropAreas=[P(1,5) P(2,5) P(3,5)].*U
        PerpDropPercentError=-100.*[(TruePeakAreas(1)-PerpDropAreas(1))./TruePeakAreas(1) (TruePeakAreas(2)-PerpDropAreas(2))./TruePeakAreas(2) (TruePeakAreas(3)-PerpDropAreas(3))./TruePeakAreas(3)] 
%         MeanAbsCLSerror=mean(abs(PercentError));
%         MeanAbsPDerror= mean(abs(PerpDropPercentError));
          disp(['Area measurement by CLS is better than perp drop by factor of ' num2str(mean(abs(PerpDropPercentError./CLSPercentError)),2)])
    end
    if sizeP(1)==2
        disp('Only 2 peaks detected by autopeaks.m. ')
    end
    if sizeP(1)==1
        disp('Only 1 peak detected by autopeaks.m.')
    end   
    if sizeP(1)>3
        disp('More than 3 peaks detected by autopeaks.m.')
    end  
catch
    disp('Perpendicular drop method failed.')
end

figure(6)
contourf(CS)
xlabel('Wavelength') 
ylabel('Time, minutes')
zlabel('Absorbance, mAU')
title('Contour plot of sample (Z axis is absorbance)')
