% Copyright 1996 Thomas Mc Avoy and Dong Dong
% This write up describes how to use MATLAB routines that implement the
% nonlinear principal component analysis (NLPCA) method [1] which is
% based on the principal curves algorithm [2]. The package consists of 4
% routines:
% npca1.m: main routine for NLPCA data calculation. 
% local1.m: implement locally weighted regression.
% pjord1.m: project data and reorder the data for locally weighted
%            regression.
% nlpcadem1.m: a program to demonstrate NLPCA 
% After the NLPCA results are calculated, the NNPLS package can be used to
% calculate a neural net, NLPCA model.

% The following references describe the approach used
% [1] Dong, D & T. J. McAvoy, "Nonlinear Principal Component Analysis
%     Based on Principal Curves and Neural Networks",  Computers and
%     Chemical Engineering , 1996, Vol 20, 65-78.
% [2] Hastie, T & Stuetzle, W., "Principal Curves", Journal of the 
%     American Statistical Association , 1989, Vol 84, 502-516.

  function  [NN1,NN2,T,DataR] = npca1(Data,span,ncom)
% This program calculates NLPCA data based on the Principal Curves algorithm.
% Data : A matrix for NLPCA calculation.
% span : Define a local range, usually between 0.4 - 0.6. See below
% ncom : Number of nonlinear principal components.
% T    : Nonlinear principal score matrix
% DataR: Reconstructed data
% NN1  : Data matrix for training the first three-layer neural net
%        NLPCA model
% NN2  : Data matrix for training the second three-layer neural net
%        NLPCA model
% function [NN1,NN2,T,DataR]= npca1(Data,span,ncom)

		[m,n]=size(Data);
		nscorse=zeros(m,ncom);
		if ncom > n
			error(' No. of NLPC must be <= No. of variables')
		end
%	span=round(span*m);% This is approach recommended by Hastie and Stuetzle
% in this code a fixed number of points, 20 with 10 on each side of the point
% of interest are used for locally weighted regression.  If more points are 
% desired span can be changed.  An even number should be used.  It is
% assumed that the number of data points is > 20.

%  Scale data to zero mean and unit variance
        [data,datam,datas]=auto(Data);
%  Calculate the explained variation by linear PC
        [u,ss,v]=svd(data,0);
        s=diag(ss);
        s=s.^2;
		lpercent=s*100/sum(s);
		disp('  ')
		disp('Percent Variance Captured by PCA')
		disp('  ')
		disp('     PC#    This PC    Total ')
		disp([(1:min(n,m))' lpercent cumsum(lpercent)]) 
%  NLPCA calculation begins here
		res=data;
        sgmf=zeros(m,n);
%
        for count=1:ncom
%
%  Initialize with first linear principal component

        [u,s,v]=svd(res,0);
        score=u*s;
        [score(:,1),kkk]=sort(score(:,1));
		res1=res(kkk,:);% res1 give data sorted by score 1
	  	f=score(:,1)*(v(:,1))';
		distsav=sum(sum((f-res).^2))/(m*n);
		distnew=distsav;
		iter=1;

		
% when distnew/distsav<.01 convergence is assumed

	while 1	
		out=sprintf('iteration # = %6.3f',iter);
		disp(out)
	

%  Carry out locally weighted regression
%  Get a better result for f
		lmd=zeros(1,m);
        for i=2:m
           lmd(i)=lmd(i-1)+norm(f(i,:)-f(i-1,:));
        end
		[f1]=local1(lmd,res1,span);
		
		
%  Project data to f and reorder data
	     [lmd,res1,distnew,kkk]=pjord1(f,res1,span,kkk);
		 if abs((distnew-distsav)/distsav)>.01&iter<10
% only allow 10 iterations
		 	distsav=distnew;
			iter=iter+1;
		 else
			disp('convergence')
			disp(' ')
			break
		 end
	end

%  Carry out final locally weighted regression
%  Get final result for f

		lmd=zeros(1,m);
        for i=2:m
           lmd(i)=lmd(i-1)+norm(f(i,:)-f(i-1,:));
        end
		[f1]=local1(lmd,res1,span);

	     
%  Ready for next NPC calculation


	  	nscore(1:m,count)=lmd(1:m)';

		
% the approach above has used sorted data, res1.  to calculate a 
% residual for the next NPC the results have to be unsorted.

		[k1,k2]=sort(kkk);
		f1=f1(k2,:);
		sgmf(1:m,:)=sgmf(1:m,:)+f1;% sigf needs same order as original data
		nscore(1:m,count)=nscore(k2,count);  
      	res=data-sgmf;
	  
%  Calculate explained variance by nonlinear PC

           [ur,sr,vr]=svd(res,0);
           percentcum(count)=(1.-sum(sum(sr.^2))/sum(sum(ss.^2)))*100;
           stex=sprintf('The frist %g NLPCs explain %g percent of variance',count,percentcum(count));
		   disp(stex)
		   disp(' ')
	end
%
		kper=0;
		for i=1:ncom
	    	percent(i)=percentcum(i)-kper;
	    	kper=percent(i)+kper;
		end
        disp('  ')
        disp('Percent Variance Captured by NLPCA')
        disp('  ')
        disp('     NLPC#    This NLPC    Total ')
        disp([(1:ncom)' percent'  percentcum']) 

%  Unscale data
		for i=1:m
	   		sgmff(i,:)=sgmf(i,:).*datas+datam;
		end
%  Save results
		T=nscore;
		DataR=sgmff;
		NN1(:,1:ncom)=nscore;
		NN1(:,ncom+1:ncom+n)=Data;
		NN2(:,1:ncom)=nscore;
		NN2(:,ncom+1:n+ncom)=sgmff;

