function [evals,evecs,efuns,edata]=pca(D,meantrend,varnorm,k_rank,method,reportflag)
%PCA - Principal Component Analysis of data matrix
% PCA performs principal component (EOF) analysis on input 
% time series, using either the SVD method or the COVARIANCE
% MATRIX method, depending on the size of the problem. See Emery 
% and Thompson or Preisendorfer for details.
% 
% Input: D - a matrix of observation time series, one series 
% per column. In the case that the columns represent 
% different physical quantities, the varnorm flag below
% can be used to normlize the column variance, by 
% dividing each series by it's demeaned std.
% meantrend - binary flag (0|1) specifying whether (or not)
% to remove mean and linear trends from the column
% data; the default is to remove both mean and
% linear trend.
% varnorm - binary flag (0|1) specifying whether (or not)
% the code should perform variance normaliation
% of the observation matrix D columns. The default
% is NOT to normalize the variance.
% k_rank - number of eigenvalues/vectors/functions to return.
% If k_rank>rank(D), then this option does nothing.
% If k_rank==0, then all eig's are returned; this is
% the default. Inputting k_rank < rank(D) amounts
% to returning the rank-k approximation to the data 
% matrix. k_rank can also be the string 'SIG', which 
% tells the routine to determine the number of 
% statistically signifigant modes and return that
% k_rank approximation. Specifying k_rank to anything 
% other than 0 forces the method to be 'SVD'.
% method - override the method selection; either 'COV','SVD'
% or 'SEL' to let the code select the method
% report - binary flag (0|1) to write a "report" to stdout 
% of the "success" of the analysis. Default is no report. 
%
% Output:
% evals - vector of eigenvalues of the covariance matrix (D'*D)
% evecs - matrix of eigenvectors of the covariance matrix
% efuns - matrix of eigenfunctions (amplitude functions, e.g.)
% edata - k_rank approximation to the data matrix; this matrix
% will only be returned if 0 < k_rank < rank(D).
%
% NOTES: 1) Use COMP_VAR_EL_PARAMS to rotate vector series onto principal
% axes if data resprsents vector series.
%
% Call as: [evals,evecs,efuns]=pca(D);
% OR: [evals,evecs,efuns]=pca(D,meantrend);
% OR: [evals,evecs,efuns]=pca(D,meantrend,varnorm);
% OR: [evals,evecs,efuns,edata]=pca(D,meantrend,varnorm,k_rank);
% OR: [evals,evecs,efuns,edata]=pca(D,meantrend,varnorm,k_rank,method);
% OR: [evals,evecs,efuns,edata]=pca(D,meantrend,varnorm,k_rank,method,report);
%

if nargin==0
disp(' Call as: [evals,evecs,efuns]=pca(D);');
disp(' [evals,evecs,efuns]=pca(D,meantrend);');
disp(' [evals,evecs,efuns]=pca(D,meantrend,varnorm);');
disp(' [evals,evecs,efuns,edata]=pca(D,meantrend,varnorm,k_rank);');
disp(' [evals,evecs,efuns,edata]=pca(D,meantrend,varnorm,k_rank,method);');
disp(' [evals,evecs,efuns,edata]=pca(D,meantrend,varnorm,k_rank,method,report);');
evals=[];evecs=[];efuns=[];
return
end

%% PROCESS INPUT ARGUMENTS
if isstr(D)&strcmp(lower(D),'refs')
disp('Preisendorfer, PCA in Meteorology and Oceanography, Elsevier, 1988.');
disp('Overland and Preisendorfer, Signifigance test for PCA, Monthly Weather Review, 1982.');
disp('Kelly, JPO, 1988');
disp('Emery and Thompson, Data Analysis Methods in Physical Oceanography, 1998.')
evals=[];evecs=[];efuns=[];edata=[];
return
end
if exist('meantrend')
if meantrend~=0 & meantrend~=1
   error('MEANTREND flag to PCA must be 0|1.')
end
else
meantrend=1;
end

if exist('varnorm')
if varnorm~=0 & varnorm~=1
   error('VARNORM flag to PCA must be 0|1.')
end
else
varnorm=0;
end
sig_test='no';
if exist('k_rank')
if isstr(k_rank)
   if ~strcmp(lower(k_rank),'sig')
   error('if K_RANK to PCA is a string, it must be the string ''SIG''.')
   end
sig_test='yes';
end
if k_rank<0
   error('K_RANK to PCA >=0 or the string ''SIG''.')
end
else
k_rank=0;
end

if exist('method')
if ~isstr(method)|((~strcmp(method,'COV'))&...
       (~strcmp(method,'SVD'))&...
       (~strcmp(method,'SEL')))
   error('METHOD to PCA must be ''COV'' or ''SVD'' or ''SEL''')
end
else
method='SEL';
end
if isstr(k_rank),method='SVD';,end
% Check the number of outgoing arguments.
if nargout<4
error('Rank-k data matrix requested, but <4 output arguments provided.')
end

if exist('reportflag')
if reportflag~=0 & reportflag~=1
   error('REPORT flag to PCA must be 0|1.')
end
else
reportflag=0;
end

% Get number of series (spatial points, variables, etc.) (M, 
% column count) and length of observations (N, row count).
% Recall that the observation matrix is input with time
% proceeding down the rows, and location (variables) across the 
% columns
[M,N]=size(D'); % NOTE THE TRANSPOSE !!! 

% The relative sizes of M and N determine whether we compute
% the EOFs from the covariance matrix or whether we use the
% SVD method, unless overridden by user on input. 
if strcmp(method,'SEL')
if M>N
method='SVD';
else
method='COV';
end
end

% Regardless of the method, we form a "data" matrix DT, which is
% de-meaned and linearly de-trended, if the user has specified
% meantrend==1.

DM=D;
if meantrend==1
% Demean data matrix D
if reportflag
disp(' ')
disp(['Before: Mean = ',num2str(mean(D))])
disp(['De-meaning, De-trending data'])
end
[DM,MEAN]=demean(DM);

% Some de-trending needs to happen here; atleast linear removal
% Compute linear trend for each variable
[DM,TRENDCOEFFS]=detrend2(DM);
end

% So far, observations have been de-meaned and linearly
% de-trended. If the observation series
% represent different physical quantities (temp and sal, for 
% example), then the user should have specified varmorm==1, so ...
if varnorm
STD=std(DM);
if reportflag
disp(' ')
disp(['Normalizing variance.'])
disp(['Before: Min STD = ',num2str(min(STD))])
disp([' Max STD = ',num2str(max(STD))])
end
[DM,STDD]=normvar(DM);
if reportflag
STD=std(DM);
disp(['After : Min STD = ',num2str(min(STD))])
disp([' Max STD = ',num2str(max(STD))])
end
end

% At this point, we transpose the data matrix to put the
% time along the columns;
DM=DM';

switch method
case 'COV'
% Form DM times DM-transpose, DDT
DDT=(1/(N-1))*DM*DM';% This is the covariance matrix
flops(0);tic;
[U,d]=eig(DDT);
FLOPS=flops;
TOC=toc;
[evals,iperm]=sort(diag(d));
iperm=flipud(iperm); % Flip to get descending order
evals=flipud(evals); 
evecs=U(:,iperm); % Permute the columns of U
efuns=DM'*evecs;
edata=[]; % Rank-k data cannot be returned if method=='COV'
case 'SVD'
flops(0);tic;
if reportflag,disp('Performing SVD...'),end 
[U,S,V]=svd(DM);
FLOPS=flops;
TOC=toc;
% now, we compute the eigenvalues from the singular values;
% they are the squares of the SVs, scaled by the length of 
% the time series -1, as in the covariance, matrix.
evals=(1/(N-1))*diag(S).^2;
% evecs are in U
evecs=U;
% Eigenfunctions
efuns=V*S';

DMrank=length(diag(S)>eps);

% Rank-k Approximation to the input data matrix
if isstr(k_rank) % then it must be 'SIG', and determine the
% statistically signifigant singular values.
     
if reportflag,disp('Determining signifigance...'),end
% set k_rank according to the following signifigance test.
% The sig. test is based on a Monte Carlo simulation.
% See Overland and Preisendorfer, Monthly Weather Review,
% Jan. '82, Vol. 110, No. 1 for details, and notation.
% Form the normalized eigenvalue statistic:
T=evals/sum(evals); % eqn 1

% Form 100 replicates of a length DMrank by 
% M gaussian random noise matrix
for r=1:100
F=randn(DMrank,M);
delta(r,:)=(eig(F*F'))';
end
for r=1:100
UU(r,:)=delta(r,:)/sum(delta(r,:));
end
for j=1:M
UU(:,j)=sort(UU(:,j));
end
ik=find(T>UU(95,:)');
k_rank=ik(length(ik));
end
if k_rank>0 & (k_rank<=min(M,N))
%return the k_rank approximation to the input data matrix.
edata=(U(:,1:k_rank)*S(1:k_rank,1:k_rank)*V(:,1:k_rank)')';
if varnorm % add back in the std
edata=edata.*(ones(size(edata(:,1)))*STDD);   
end
if meantrend % add back in the means and trends
   tt=(1:length(edata(:,1)))';tt=tt-mean(tt);
   edata=edata+tt*TRENDCOEFFS;
edata=edata+ones(size(edata(:,1)))*MEAN;
end
else % k_rank=0
edata=[];
end
end

if nargout==1 | nargout==0
efuns=[];evecs=[];
elseif nargout==2
efuns=[];
elseif nargout==3
edata=[];
end

if reportflag
disp(' ')
disp(['Method = ' upper(method)]);
disp(['Flops = ' int2str(FLOPS)]);
disp(['O(MNN) = ' int2str(M*N*N)]);
disp(['Time = ' num2str(TOC)]);
disp(['Matrix Rank = ' num2str(DMrank)]);
disp(' ')
if strcmp(sig_test,'yes')
disp(['Signifigant Rank (95%-conf level) = ' int2str(k_rank)])
end
if k_rank>0
disp(['Rank-' int2str(k_rank) ' data matrix approximation returned in 4th output argument.'])
disp(' ')
end
maxerr=max(max(abs(evecs*evecs')-eye(size(evecs))));
disp(['Max error in orthogonality of eigenvectors: ' num2str(maxerr)])
disp(' ')
disp(' Variance Table')
if strcmp(sig_test,'yes')
disp(' (** indicates signifigance)')
end
disp(' Mode E-val %Var ')
nr=min(M,N);
pvar=cumsum(evals)/sum(evals)*100;
for i=1:nr
if i<=k_rank & strcmp(sig_test,'yes')
fprintf(' ** %d %.3f %6.2f\n',[ i evals(i) pvar(i)])
else
fprintf(' %d %.3f %6.2f\n',[ i evals(i) pvar(i)])
end
end
fprintf('Total %.3f\n',sum(evals(1:nr)))
end

%
% Brian O. Blanton
% Department of Marine Sciences
% Ocean Processes Numerical Modeling Laboratory
% 12-7 Venable Hall
% CB# 3300
% University of North Carolina
% Chapel Hill, NC
% 27599-3300
%
% 919-962-4466
% blanton@marine.unc.edu
%
% Summer 1999
%