als_ka.m

GUI for Multivariate Image Analysis of 4-dimensional data (Matlab Code)

function[copt,sopt]=als_ka(d, Ncomp,nexp,datamod, N);
% [copt,sopt]=als_ka(d, Ncomp,nexp,datamod, N);
%
%  function for 3-dimensional MCR with following contraints:
% - nonnegativity in both absorption and concentration modes
% - fast NNLS is used in both modes
% - no unimodality contraints
% - no closure contraints
% - no shape contraints
% - no normalization
% - no triliniar structure
% - random initialization is used
%
% convergence criteria:
% - 50 iterations
% - tolerance of 0.1
%
%input data:
%  d(nxm) - original data set with augmented separate experiments:
%   - data_col=[a1; a2; a3; a4; a5]; % column-wise
%   - data_row=[a1 a2 a3 a4 a5]; % row-wise
% x0 - initial estimation (ncomp x m)
% nexp - number of combined matrices/experiments
% datamod=1 for column-wise
% datamod=2 for row-wise
% N-number of rows/columns in each separate data set
%
% the code written by Roma Tauler and Anna de Juan (e-mail [e-mail]) 
% is modified by fixing constraints
%
% modified by K.Artyushkova
% [e-mail]
% February 2005

%% check data
% check dimensions of initial estimates
[nrow,ncol]=size(d);
x0=rand(Ncomp, ncol); % generate random initial estimates 
[nrow2,ncol2]=size(x0);
if nrow2==nrow,	nsign=ncol2; ils=1;end
if ncol2==nrow, nsign=nrow2; x0=x0'; ils=1; end 
if ncol2==ncol, nsign=nrow2; ils=2;end
if nrow2==ncol, nsign=ncol2; x0=x0'; ils=2; end
if ils==1,
	conc=x0;
	[nrow,nsign]=size(conc);
	abss=conc\d;
end
if ils==2,
	abss=x0;
	[nsign,ncol]=size(abss);
	conc=d/abss;
end


%% INITIALIZATIONS
nit=50;
tolsigma=0.1;
isp=ones(nexp,nsign);
csel=[];
ssel=[];
vclos1=0;
vclos2=0;
scons='y'; % all the spectra matrices the same constraints
ccons='y'; % all the concentration matrices the same constraints
niter=0;% iterations counter
idev=0;% divergence counter
idevmax=10;% maximum number of diverging iterations
answer='n'; % default answer
ineg=0;% used for non-negativity constraints
imod=0;% used for unimodality constraints
iclos0=0;% used for closure constraints
iassim=0;% used for shape constraints
matr=1;% one matrix
matc=1;% one matrix
vclos1n=0;% used for closure constraints
vclos2n=0;% used for closure constraints
inorm=0;% no normalization (when closurfe is applied)
type_csel=[]; %no equality/lower than constraints in concentrations 
type_ssel=[]; %no equality/lower than constraints in spectra
if datamod == 1 % columnwise
   matr = 1;
   ncinic(1)=1;
   ncfin(1)=ncol;
   matc = nexp;		
   nrinic(1)=1;
   for i=1:matc,
       nrfin(1)=N;
      if i>1,
         nrinic(i)=nrfin(i-1)+1;
         nrfin(i)=nrinic(i)+N-1;
      end
      ncinic(i)=1;
		 ncfin(i)=ncol;
	end
end
if datamod == 2 %row-wise
   matc = 1;
   nrinic(1)=1;
   nrfin(1)=nrow;
   matr = nexp;
   ncinic(1) = 1;
	for i=1:matr,
		ncsol(i)=N;
		ncfin(1)=ncsol(1);
		if i>1, ncinic(i)=ncfin(i-1)+1;
			ncfin(i)=ncinic(i)+ncsol(i)-1; end

	end
end	

%% contsraints set-up:
wcons=1; % nonnegativity;
if matc>1
ccons = 'y';  %the same constraints to all matrices
end
if matr>1
scons = 'y';  %the same constraints to all matrices
end
c1 = find(wcons == 1);
ineg=2; % both concnetration and spectra non-negative
ialgs=2; % fast nnls algorithm;
spneg = ones(nsign,matr);
ialg=2; % fast nnls algorithm;
cneg = ones(matc,nsign);
inorm=0; % no normalization
ishape=0; % not trilinear shape
trildir=99;
spetric=zeros(1,nsign);
spetris=zeros(1,nsign);
gr='n'; % no graphic output
c3 = find(wcons == 3);
cons4 = find(wcons == 4);
cons5 = find(wcons == 5);

%% REPRODUCTION OF THE ORIGINAL DATA MATRIX BY PCA
% dn is the experimental matrix and d is the PCA reproduced matrix
dn=d;
[u,s,v,d,sd]=pcarep(dn,nsign);
sstn=sum(sum(dn.*dn));
sst=sum(sum(d.*d));
sigma2=sqrt(sstn);


%% D) STARTING ALTERNATING CONSTRAINED LEAST SQUARES OPTIMIZATION
while niter < nit
niter=niter+1;

% E) ESTIMATE CONCENTRATIONS (ALS solutions)
conc=d/abss;

% CONSTRAIN APPROPRIATELY THE CONCENTRATIONS
% non-negativity
if ~isempty(c1)
if ineg==1|ineg==2
for i =1:matc
kinic=nrinic(i);
kfin=nrfin(i);
conc2=conc(kinic:kfin,:);
if ialg==0
	for k=1:nsign,
	if cneg(i,k) ==1
		for j=1:kfin+1-kinic,
			if conc2(j,k)<0.0,
				conc2(j,k)=0.0;
			end
		end
	end
	end
end
if ialg==1
	for j=kinic:kfin
		if cneg(i,:) == ones(1,size(isp,2))
		x=lsqnonneg(abss',d(j,:)');
		conc2(j-kinic+1,:)=x';
	end
	end
end
if ialg==2  
	for j=kinic:kfin
		if cneg(i,:) == ones(1,size(isp,2)) 
		x=fnnls(abss*abss',abss*d(j,:)');
		conc2(j-kinic+1,:)=x';
	end
	end
end
conc(kinic:kfin,:) = conc2;
end
end
end

% trilinearity
if ishape>=1
if trildir==1|trildir==3
   	for j=1:nsign,
   		if spetric(j)==1,
   			[conc(:,j),t]=trilin(conc(:,j),matc,ishape);
            totalconc(j,1:matc)=t;
            if totalconc(j,1)>0,
	   			rt(j,1:matc)=totalconc(j,1:matc)./totalconc(j,1);
            else
               rt(j,1:matc)=totalconc(j,1:matc);
            end
         end
   	end
end
end

% zero concentration species
if matc>1
for i=1:matc,
	for j=1:nsign,
		if isp(i,j)==0,
				conc(nrinic(i):nrfin(i),j)=zeros(nrsol(i),1);
		end
	end
end
end

% unimodality
for i = 1:matc
kinic=nrinic(i);
kfin=nrfin(i);
conc2=conc(kinic:kfin,:);
	if imod==1|imod==3,
		for ii=1:nsign,
			if spmod(i,ii)==1,
				conc2(:,ii)=unimod(conc2(:,ii),rmod,cmod);
			end
        end
    end
conc(kinic:kfin,:)=conc2;
end

% EQUALITY CONSTRAINTS IN CONC
if ~isempty(cons4)
   if type_csel==0,conc(iisel)=csel(iisel);end
   if type_csel==1 
      for ii=1:size(iisel),
         if conc(iisel(ii))>csel(iisel(ii)),conc(iisel(ii))=csel(iisel(ii));end
      end
   end
end

% closure
if ~isempty(c3)
if dc == 1
for i = 1:matc
kinic=nrinic(i);
kfin=nrfin(i);
conc2=conc(kinic:kfin,:);
if iclos(i)==1 | iclos(i)==2,
if tclos1(i) == 0 
	vclos1n=vclos1(kinic:kfin,1);
end
if iclos(i) ==2 & tclos2(i)==0
	vclos2n=vclos2(kinic:kfin,1);
end
[conc2]=closure(conc2,iclos(i),sclos1(i,:),iclos1(i),tclos1(i),tclos2(i),sclos2(i,:),iclos2(i),vclos1n,vclos2n);
end
conc(kinic:kfin,:)=conc2;
end
end
end


% QUANTITATIVE INFORMATION FOR THREE-WAY DATA SETS
% recalculation of total and ratio concentrations if ishape=0 or niter=1
if ishape==0 | niter==1,
   	for j=1:nsign,
   		for inexp=1:matc,
   			totalconc(j,inexp)=sum(conc(nrinic(inexp):nrfin(inexp),j));
         end
         if totalconc(j,1)>0,
	   			rt(j,1:matc)=totalconc(j,1:matc)./totalconc(j,1);
         else
               rt(j,1:matc)=totalconc(j,1:matc);
         end
   	end
end

% areas under concentration profiles
area=totalconc;

%% ESTIMATE SPECTRA (ALS solution)
abss=conc\d;

% non-negative spectra
if ineg ==2 |ineg==3,
for i = 1:matr
kinic = ncinic(i);
kfin = ncfin(i);
abss2 = abss(:,kinic:kfin);
	if ialgs==0,
		for k=1:nsign,
		if spneg(k,i)==1
			for j=1:kfin+1-kinic,
    				if abss2(k,j)<0.0,
    					abss2(k,j)=0.0;
    				end
			end
		end
		end
    end
	if ialgs==1,
		for j=kinic:kfin,
			if spneg(:,i)== ones(size(isp,2),1)
			abss2(:,j-kinic+1)=lsqnonneg(conc,d(:,j));
			end
		end
    end
	if ialgs==2, 
		for j=kinic:kfin,
			if spneg(:,i)== ones(size(isp,2),1)			
			abss2(:,j-kinic+1)=fnnls(conc'*conc,conc'*d(:,j));
			end			
		end
	end
abss(:,kinic:kfin)=abss2;
end
end


% trilinearity
if ishape>=1,
if trildir==2|trildir==3
   	for j=1:nsign,
   		if spetris(j)==1,
   			[absst,t]=trilin(abss(j,:)',matr,ishape);
			abss(j,:)=absst';
   		end
   	end
end
end

% constrain the unimodality of spectra
for i = 1:matr
kinic = ncinic(i);
kfin = ncfin(i);
abss2 = abss(:,kinic:kfin);
if imod==2|imod==3,
	for j=1:nsign,
		if spsmod(i,j)==1
		dummy=unimod(abss2(j,:)',smod,cmod);
		abss2(j,:)=dummy';
		end
	end
end
abss(:,kinic:kfin)=abss2;
end

% EQUALITY CONSTRAINTS FOR SPECTRA
if ~isempty(cons5)
   if type_ssel==0,abss(jjsel)=ssel(jjsel);end
   if type_ssel==1
       for jj=1:size(jjsel),
         if abss(jjsel(jj))>ssel(jjsel(jj)),abss(jjsel(jj))=ssel(jjsel(jj));end
      end
   end
end

% closure in spectra (in case of inverted analysis D'=SC)
if ~isempty(c3)
if dc==2,
for i = 1:matr
kinic = ncinic(i);
kfin = ncfin(i);
abss2 = abss(:,kinic:kfin);
if iclos(i)==1 | iclos(i)==2,
if tclos1(i) == 0 
	vclos1n=vclos1(kinic:kfin,1);
end
if iclos(i) ==2 & tclos2(i)==0
	vclos2n=vclos2(kinic:kfin,1);
end
abst = closure(abss2',iclos(i),sclos1(i,:),iclos1(i),tclos1(i),tclos2(i),sclos2(i,:),iclos2(i),vclos1n,vclos2n);
abss2=abst';
end
abss(:,kinic:kfin) = abss2;
end
end
end

% NORMALIZATION OF SPECTRA
% equal heighth
if inorm==1,
	maxabss=max(abss');
		for i=1:nsign,
			abss(i,:)=abss(i,:)./maxabss(i);
		end
end
% equal length
if inorm==2, abss=normv2(abss); end

%% CALCULATE RESIDUALS
res=d-conc*abss;
resn=dn-conc*abss;

%% OPTIMIZATION RESULTS
u=sum(sum(res.*res));
un=sum(sum(resn.*resn));
sigma=sqrt(u/(nrow*ncol));
sigman=sqrt(un/(nrow*ncol));
change=((sigma2-sigma)/sigma);
if change < 0.0,
	idev=idev+1;
else,
	idev=0;
end
change=change*100;
sstd(1)=sqrt(u/sst)*100;
sstd(2)=sqrt(un/sstn)*100;
r2=(sstn-un)/sstn;

% If change is positive, the optimization is working correctly
if change>0 | niter==1,
	sigma2=sigma;
	copt=conc;
	sopt=abss;
	sdopt=sstd;
	ropt=res;
	rtopt=rt';
	itopt=niter;
        areaopt=area;
	r2opt=r2;   
end

% test for convergence within maximum number of iterations allowed
if abs(change) < tolsigma,
   
%  finish the iterative optimization because convergence is achieved
   figure(1)
   subplot(2,1,1);plot(copt);title('conc profile in optimal iteration');
   subplot(2,1,2);plot(sopt');title('pure spectra in optimal iteration');
   pause(3)
close(1)
   return         % 1st return (end of the optimization, convergence)
end
%  finish the iterative optimization if divergence occurs 20 times consecutively
if idev > 20,
    figure(1)
   subplot(2,1,1);plot(copt);title('conc profile in optimal iteration');
   subplot(2,1,2);plot(sopt');title('pure spectra in optimal iteration');
   pause(3)
close(1)
   return          % 2nd return (end of optimization, divergence)      
   
end
% this end refers to number of iterations initially proposed exceeded
end

% finish the iterative optimization if maximum number of allowed iterations is exceeded
figure(1)
subplot(2,1,1);plot(copt);title('conc profile in optimal iteration');
subplot(2,1,2);plot(sopt');title('pure spectra in optimal iteration');
pause(3)
close(1)
return          % 3rd return (end of optimization, number of iterations exceeded)