
%--------------------------------------------------------------------%
% script to compute eigenvalues and eigenvectors of time-series data %
%--------------------------------------------------------------------%
% Inputs
%      - time-series data W, W is N by M where the model size N is 
%        assumed to be much larger than the number M of time samples
%      - either var_limit or k to determine how much variance is kept,
%        e.g., var_limit=0.75 keeps enough eigenmodes to explain 75% 
%        of the total variance or k=1:10 keeps ten modes.
% Outputs
%      - z eigenvectors
%      - d eigenvalues
%      - total variance
%--------------------------------------------------------------------%
% Method to compute eigenvalues and eigenvectors of c0 = W*W'
% c0 = W*W' = z*d*z';
% C  = W'*W = z_s*d*z_s' <- same eigenvalues as c0
% W'*c0*W = C^2
% W*W'*c0*W*W' = c0^3 = W*C^2*W' = (W*z)*d^2*(W*z)'
%              = (W*z*sqrt(d))*d^3*(W*z*sqrt(d))'
% z = (W*z*inv(sqrt(d)))
%------------------------------------------------------------------
% returns z and d
%------------------------------------------------------------------
iok = ~isnan(W(:,1));
if nnz(~iok)==0
  C = W'*W; %does this use less memory?
else
  C = W(iok,:)'*W(iok,:);
end

%
[z_s,d]=sorted_eig(C); %sorted eigenvalues and eigenvectors of C
total=trace(d);
if ~exist('k')
  k=find(cumsum(diag(d))/total <= var_limit); %IMPORTANT!!
end
var_expl = sum(diag(d(k,k)))/total;
z_s=z_s(:,k);d=d(k,k); 
%time_series = z_s*diag(1./sqrt(diag(d)));
%z = W*time_series; 
time_series = (z_s*diag(sqrt(diag(d))))';
z=repmat(NaN,size(W,1),length(k));
z(iok,:) = W(iok,:)*time_series'*diag(1./diag(d));

