Home › Code Snippets

Self Modeling Curve Resolution (SMCR) Procedure

Posted December 12th, 2011 by Igor
in
Average rating
(0 votes)

My name is David Wilcox, and I am a graduate student in physical chemistry at Purdue University under Prof. Dor Ben-Amotz. We use self-modeling curve resolution (SMCR) to study small perturbations of water's vibrational mode frequencies around dissolved solutes (relative to bulk water). SMCR is one of the multivariate curve resolution algorithms which we use to obtain aqueous hydration shell spectra.

// Input: datamat is a 2D matrix whose columns are spectra obtained from two-component 
// solutions of different concentration. For example, the first spectrum can be that of a 
// pure solvent and the remaining spectra can have various (at least one) solute
// concentrations.
// There is no need to normalize the spectral measurements -- this is done
// inside the function.
// Output: component1 and component2, one of which is the minimum area solute correlated
// spectrum, while the other is the minimum area spectrum pertaining to the pure solvent.
// The pure solvent spectrum itself can be re-generated from the SMCR scores pertaining
// to the pure solvent spectrum. For example, if the first spectrum (in column 0) is the pure 
// solvent then the following command will re-generate the pure solvent spectrum (solvent).
// from the SMCR scores and PC loadings.
// solvent=score1[0]*PCLoadings1+score2[0]*PCLoadings2
// Requires: IP 6.22 or newer.
// Reference: Lawton and Sylvestre, Technometrics, volume 13, page 617, 1971
// Procedure written by: David Wilcox (wilcoxds at purdue dot edu)
 
 
Function SMCR(datamat)
          wave datamat
 
  		// the following command normalizes and transposes the input matrix.
  		MatrixOP/O/FREE inMat=scaleCols(datamat,(powr(sumcols(datamat),-1)^t))^t
  		Variable i,numRows
 
           //Calculate scores and loadings from SVD 
          MatrixSVD inMat				               
          wave M_U, W_W, M_VT 
          MatrixTranspose M_VT
          MatrixOP/O W_W=diagonal(W_W)
          MatrixOP /FREE/O PCLoadings1=col(M_VT,0)
          MatrixOP/FREE/O PCLoadings2=col(M_VT,1)
          MatrixOP/FREE/O aug_scores = M_U x W_W
          MatrixOP/FREE/O score1=col(aug_scores ,0)		// PC1 scores for each input specxtrum
          MatrixOP/FREE/O score2=col(aug_scores ,1)		// PC2 scores for each input specxtrum
 
          //Find minimum positive and maximum negative of loadings ratio (slopes)
          MatrixOP/O/FREE Lratio = PCLoadings1/PCLoadings2
          numRows=DimSize(Lratio,0)
          variable m1,m2
          i=0
          m1=-10^100
          m2=10^100
          do
                if (Lratio[i] > 0)
          	    m1 = max(-Lratio[i],m1)
          	elseif (Lratio[i] < 0)
          	    m2 = min(-Lratio[i],m2)
          	endif
          	i +=1
          while (i<numRows)
 
          // Linear regression of scores ratio
          CurveFit /Q line, score2/x=score1
          wave W_coef
          Variable m=W_coef[1]
          Variable b=W_coef[0]
 
          // Intersection of scores and loadings
          make/FREE/O/N=(2,2)  x1x2num={{b,1},{0,1}}
          make/FREE/O/N=(2,2)  x1y1denom={{-m,1},{-m1,1}}
          MatrixOP/FREE/O x1= Det(x1x2num)/Det(x1y1denom) // (x1,y1) is one end point of the PC score line
          make/FREE/O/N=(2,2)  y1num={{-m,b},{-m1,0}}
          MatrixOP/FREE/O y1= Det(y1num)/Det(x1y1denom)
          make/FREE/O/N=(2,2)  x2y2denom={{-m,1},{-m2,1}}
          MatrixOP/FREE/O x2= Det(x1x2num)/Det(x2y2denom) // (x2,y2) is one end point of the PC score line
          make/FREE/O/N=(2,2)  y2num={{-m,b},{-m2,0}}
          MatrixOP/FREE /O y2  = Det(y2num)/Det(x2y2denom) 
 
          // Make pure component spectra
          Variable x1v=x1[0],x2v=x2[0],y1v=y1[0],y2v=y2[0]
          MatrixOP/O	component1 = x1v*PCLoadings1 + y1v*PCLoadings2
          MatrixOP/O	component2 = x2v*PCLoadings1 + y2v*PCLoadings2
          KillWaves/Z M_U,W_W,M_VT,W_coef,W_sigma
end

Back to top