function [Xeofs,Xpcs,Var,Xrecon]=eof(X,neof)
% function [EOFs,PCs,Var,Xrecon]=eof_analysis(X,neof)
% Wrapper function to perform PCA of a field X
% with TWO spatial dimensions. (This code will also
% work with a SINGLE spatial dimension. But it might
% be easier to directly call 'principal_component_analysis'.)
% This function basically transforms the data matrix into
% a standard 2-d data matrix, taking into account NaNs.
% Input:
%   X: (x,y,t) or (x,t).
%   neof: number of EOF/PC to return
% Output:
%   EOFs: 3-d (x,y,e) or 2-d (x,e) matrix with spatial
%         patterns.
%   PCs:  2-d matrix (t,e) with principal components (scores)
%         in the columns
%   Var:  variance of each principal component
%   Xrecon: 3-d (x,y,t) or 2-d (x,t) matrix with reconstructed X
%           WITHOUT adding back the mean
% Note: Xpc is the projection of X onto the corresponding EOF.
% That is, Xpc(it,ie) = nsum(nsum(X(:,:,it).*Xeofs(:,:,ie)))
% Use this to project a physical field onto EOF space.

% If X only has 2 dimensions, assume second dimension is time.
% Insert 'y' dimension to conform with rest of function which
% assumes X represents multiple realizations of a 2-D field.


if ndims(X)==2
 [nx,nt]=size(X);
 X=reshape(X,[nx 1 nt]);
end

% Flatten 2-d fields into single vector
[nx,ny,nt]=size(X);
Xd=reshape(X,[nx*ny nt]); % (space,time)
% remove NaNs
inn=find(~isnan(Xd(:,1)));
Xd=Xd(inn,:); % (space, time)  

n=size(Xd,1);
%normalize
for i=1:n
% Xd(i,:)=Xd(i,:)-nanmean(Xd(i,:));
% Xd(i,:)=detrend(squeeze(Xd(i,:)),'constant');
% Xd(i,:)=detrend(squeeze(Xd(i,:)));
Xd(i,:)=zscore(detrend(squeeze(Xd(i,:)),'constant'));
end

Xd=Xd';
Xd=double(Xd);
% PCA
if nargout>3
 [EOFs,Xpcs,Var,XR]=principal_component_analysis(Xd,neof);
else
 [EOFs,Xpcs,Var]=principal_component_analysis(Xd,neof);
end  
% Xpcs (t,e)

% Reshape EOFs to physical space
Xeofs=repmat(NaN,[nx*ny neof]);
Xeofs(inn,:)=EOFs;
Xeofs=squeeze(reshape(Xeofs,[nx ny neof]));

Xrecon=repmat(NaN,[nx ny nt]);
for ie=1:neof
 for it=1:nt
   Xrecon(:,:,it)=Xrecon(:,:,it)+Xpcs(it,ie)*Xeofs(:,:,ie);
 end
end
 
if nargout>3
 XR=XR'; % (x,t)
 Xrecon=repmat(NaN,[nx*ny nt]);
 Xrecon(inn,:)=XR;
 Xrecon=squeeze(reshape(Xrecon,[nx ny nt]));
end
function Xpcs=eof_project(Xanom,Xeofs,numd)
% function Xpcs=eof_analysis(Xanom,Xeofs,numd)
% Function to project physical field (1 or 2 spatial
% dimensions) onto EOF.
% Input:
%   Xanom: (x,y,t) or (x,t) - physical ANOMALY field
%   Xeofs: 3-d (x,y,e) or 2-d (x,e) matrix with spatial
%          patterns.
%   numd: if Xanom is (x,t), MUST set ndims=1
% Output:
%   Xpcs: 2-d matrix (t,e) with projection coefficients
%         (principal components) in the columns. Number
%         of projection coefficients returned is equal
%         to the number of EOFs passed in the input argument.

% Samar Khatiwala (spk@ldeo.columbia.edu)

if nargin<3 % default is 2-d
 numd=2;
end

if ndims(Xanom)==2 & numd==1  % data is (x,t)
 Xpcs=Xeofs'*Xanom;
 Xpcs=Xpcs';
else
%   nt=size(Xanom,3)
%   neofs=size(Xeofs,3)
%   Xpcs=repmat(NaN,[nt neofs]);
%   for it=1:nt
%     Xpcs(it,:)=squeeze(nansum(nansum(Xeofs.*repmat(Xanom(:,:,it),[1 1 neofs]))))';
%   end
%   for ie=1:neofs
%     Xpcs(:,ie)=squeeze(nansum(nansum(Xanom.*repmat(Xeofs(:,:,ie),[1 1 nt]))));
%   end  
% Flatten 2-d fields into single vector
 [nx,ny,nt]=size(Xanom);
 neofs=size(Xeofs,3);  
 Xanom=reshape(Xanom,[nx*ny nt]); % (x,t)
 Xeofs=reshape(Xeofs,[nx*ny neofs]); % (x,e)  
% remove NaNs
 inn=find(~isnan(Xanom(:,1)));  
 Xanom=Xanom(inn,:); % (x,t)
 Xeofs=Xeofs(inn,:); % (x,e);
 Xpcs=Xeofs'*Xanom;
 Xpcs=Xpcs';
end


function [EOFs,PCs,Var,Xrecon]=principal_component_analysis(X,neof)
% function [EOFs,PCs,Var]=principal_component_analysis(X,neof)
% Function to do a principal component analysis of
% data matrix X.
% Input:
%   X: (t,x) each row corrsponds to a sample, each column
%      is a variable. (Each column is a time series of a
%      variable.)
%   neof: number of EOF/PC to return
% Output:
%   EOFs: (x,e) matrix with EOFs (loadings) in the columns
%   PCs:  (t,e) matrix with principal components (scores) in the columns
%   Var:  variance of each principal component
%   Xrecon: (t,x) reconstructed X (WITHOUT adding back the mean)
% To reconstruct: Xrecon = PCs*EOFs'
% Notes: (1) This routine will subtract off the mean of each
%            variable (column) before performing PCA.
%        (2) sum(var(X)) = sum(Var) = sum(diag(S)^2/(m-1))

% Samar Khatiwala (spk@ldeo.columbia.edu)

if strcmp(class(X),'single')
 disp('WARNING: Converting input matrix X to class DOUBLE')
end  

% Center X by subtracting off column means
[m,n] = size(X);
%X = X - repmat(mean(X,1),m,1);
r = min(m-1,n); % max possible rank of X

% SVD
if nargin < 2
  [U,S,EOFs]=svds(X,r);
else
  [U,S,EOFs]=svds(X,min(r,neof));
end

% EOFs: (x,e)
% U: (t,e)

% Determine the EOF coefficients
PCs=U*S; % PCs=X*EOFs (t,e)

% compute variance of each PC % Modified
EV=diag(S).^2/sum( diag(S).^2 );
Var = EV*100;

% Note: X = U*S*EOFs'
%       EOFs are eigenvectors of X'*X = (m-1)*cov(X)
%       sig^2 (=diag(S)^2) are eigenvalues of X'*X
%       So tr(X'*X) = sum(sig_i^2) = (m-1)*(total variance of X)

if nargout>3
 Xrecon = PCs*EOFs'; % (t,x)
end