% erpPCA - Unrestricted covariance-based PCA with Varimax rotation
%          (cf. Kayser J, Tenke CE, Clin Neurophysiol, 2006;117(2):369-80)
%
% Usage: Y = erpPCA_Varimax4m_progress (X, varargin)
%
% Generic PCA-Varimax implementation emulating the PCA agorithms used by 
% BMDP-4M (Dixon, 1992) and SPSS 10.0 FACTOR (http://www.spss.com/tech/
% stat/Algorithms/11.5/factor.pdf). It expects a data matrix (X) of ERP 
% waveforms, with ERPs (cases) as rows and sample points (variables) as 
% columns. The routine returns the unrotated (LU) and Varimax-rotated (LR) 
% factor loadings as a variables-by-factors matrix, the rotated factor 
% scores (FSr) as a cases-by-factors matrix, and Eigenvalues and explained 
% variance as a variables-by-variance matrix (VT), with four columns 
% consisting of Eigenvalues and percentage of explained variance before 
% and after rotation.
%
% erpPCA employs Varimax4M (max. 100 iterations, 0.0001 convergence criterion,
% Kaiser's normalization; MatLab code by $jk available on request), which 
% emulates algorithms described by Harman (1967, pp. 304-308) as implemented 
% in BMDP-4M (Dixon, 1992, pp. 602-603).
%
% Output arguments:
% -----------------
%  Y     = structure consisting of:
%
%    LU  = unrotated factor loadings data matrix <vars> x <factors>
%    LR  = rotated factor loadings data matrix   <vars> x <factors>
%    FSr = rotated factor scores data matrix    <cases> x <factors>
%    VT  = variance table data matrix         <factors> x <columns>
%          column 1: unrotated factors: Eigenvalue
%          column 2: unrotated factors: explained variance
%          column 3: rotated factors: Eigenvalue
%          column 4: rotated factors: explained variance
%    G   = Varimax Simplicity G criterium matrix:
%          column 1: <#> iteration
%          column 2: Simplicity G criterion
%          column 3: difference to iteration (#-1>
%          column 4: (progressively adjusted) convergance <nTol> criterium
%    mu  = grand mean (across <vars>)
%  sigma = SD for each variable
%          (Programmer's Note: "std(X)" differs from "sqrt(diag(cov(X)))")
%  FSCFr = rotated and scaled factor score coefficients
%  rank  = matrix rank
%
% Input arguments:
% ---------------
% X      = 2-D data matrix <cases> x <vars>
%
% varargin (argument pairs consisting of <string> and <var>  
% --------
% 'Kaiser' = apply Kaiser's normalization to Varimax rotation (1 = yes [default]; 0 = no)
% 'iMax'   = specify max. number of Varimax iterations (default: 100)
% 'nTol'   = specify iteration tolerance threshold (default: -1e-06)
% 'fMax'   = max. number of factors (maxF) with factor scores (FS) returned
% 'totPer' = fMax with FS returned as percentage of total variance explained (default: 99%)
% 'facPer' = fMax with FS returned as percentage of variance explained per factor
% 'mu'     = provide grand mean (across <vars>) [default: calculate (requires more memory)]
% 'sigma'  = provide SD (across <vars>) [default: calculate (requires more memory)]
% 'Matrix' = PCA association matrix (specify one):
%                 'cov' = covariance [default]
%                  'cp' = cross-products
%            'corrcoef' = correlation
% 'rank'   = force matrix rank
%
% Copyright (C) 2003 by Jürgen Kayser (Email: kayserj@nyspi.columbia.edu)
% GNU General Public License (http://www.gnu.org/licenses/gpl.txt)
% Updated: $Date: 2003/07/08 14:00:00 $ $Author: jk $
%
% Updated: $Date: 2013/12/26 15:56:00 $ $Author: jk $
%   - limit computation of rotated factor scores
% Updated: $Date: 2017/05/10 18:22:00 $ $Author: jk $
%   - added progress bar for computing factor scores
% Updated: $Date: 2017/05/12 23:14:00 $ $Author: jk $
%   - improved code to compute factor scores
%   - revised input syntax (all optional arguments via <varargin>)
%   - added options 'totPer' and 'facPer' to limit factor score computation by variance
%   - added option 'Kaiser' to control Kaiser's normalization in Varimax rotation
%   - revised output syntax (all variables embedded in structure)
%   - added <mu>, <sigma>, <FSCFr> to output (can be subsequentially used
%     to compute factor scores via DeriveFactor Scores.m)
% Updated: $Date: 2017/05/19 12:29:00 $ $Author: jk $
%   - added options to specify <mu> and <sigma> to save memory demand
%   - clear variables from memory when no longer used
% Updated: $Date: 2018/08/30 10:35:00 $ $Author: jk $
%   - added options to specify association matrix type
% Updated: $Date: 2019/04/22 15:45:00 $ $Author: jk $
%   - fixed documentation
% Updated: $Date: 2021/01/19 20:21:00 $ $Author: jk $
%   - added option 'rank'
%   - added matrix rank <rk> to output

function Y = erpPCA_Varimax4m_progress (X, varargin)
tstart = tic;
try
 [cases, vars] = size(X);             % get dimensions of input data matrix     1
catch
 cases = 0;
 vars = 0;
end;
if (cases < 1) | (vars < 1)
 Y = [];   
 help erpPCA_Varimax4m_progress   
 return   
end;
if ~isempty(varargin)
 try, g = struct(varargin{:}); 
 catch, error('Argument error in the {''param'', value} sequence'); end; 
end;
try g.iMax;   catch g.iMax = []; end;
try g.nTol;   catch g.nTol = []; end;
try g.fMax;   catch g.fMax = []; end;
try g.totPer; catch g.totPer = []; end;
try g.facPer; catch g.facPer = []; end;
try g.Kaiser; catch g.Kaiser = []; end;
try g.mu; catch g.mu = []; end;
try g.sigma; catch g.sigma = []; end;
try g.Matrix; catch g.Matrix = []; end;
try g.rank; catch g.rank = []; end;
if isempty(g.nTol); nTol = -1e-06; else nTol = g.nTol; end;
if isempty(g.iMax); iMax = 100; else iMax = g.iMax; end;
if isempty(g.fMax); fMax = 0; else fMax = g.fMax; end;
if isempty(g.totPer); fPer = 99; else fPer = g.totPer; end;
if ~isempty(g.facPer); fMax = -1; fPer = g.facPer; end;
if isempty(g.Kaiser); Kaiser = 1; else Kaiser = g.Kaiser; end;
if isempty(g.Matrix); Matrix = 'cov'; else Matrix = g.Matrix; end;
if ~(strcmpi(Matrix,'cov')|strcmpi(Matrix,'cp')|strcmpi(Matrix,'corrcoef'))
  help erpPCA_Varimax4m_progress   
  error(sprintf('Argument error in the {''Matrix'', ''%s''} sequence',Matrix));
  return   
end    
Y.LU = [];
Y.LR = [];
Y.FSr = [];
Y.VT = [];
Y.G = [];
Y.mu = [];
Y.sigma = [];
Y.FSCFr = [];
if strcmpi(Matrix,'cov')
 disp('Covariance Matrix');   
 D = cov(X);                         % compute covariance matrix               2
elseif strcmpi(Matrix,'cp')
 disp('Cross-Products Matrix');   
 D = cp(X);                          % compute cross-products matrix           2
elseif strcmpi(Matrix,'corrcoef')
 disp('Correlation Matrix');   
 D = corrcoef(X);                    % compute correlation matrix              2
end
[EM, EV] = eig(D);                   % determine Eigenvectors and Eigenvalues  3
sigma = sqrt(diag(D));
clear D;
UL = EM * sqrt(EV);                  % determine unrotated factor loadings     4
[u, ux] = sort(diag(EV)');           % sort initial Eigenvalues, keep indices  5
u = fliplr(u); ux = fliplr(ux);      % set descending order                    6
LU = UL(:,ux);                       % sort unrotated factor loadings          7
rk = rank(corrcoef(X),1e-4);         % estimate the number of singular values  8
if ~isempty(g.rank) 
  rk = min(rk,g.rank);               % assign specified rank if less than observed rank 
end;
disp(sprintf('Matrix Rank: %d',rk));   
LU(:,rk+1:end) = [];                 % ... and remove all linearly dependent   9
u = u(1:rk); ux = ux(1:rk);          % ... components and their indices       10
s = ones(1,rk);                      % current sign of loading vectors        11
s(abs(max(LU)) < abs(min(LU))) = -1; % determine direction of loading vectors 12
LU = LU .* repmat(s,size(LU,1),1);   % redirect loading vectors if necessary  13
%RL = Varimax4M_progress(LU,200,1e-4,1);       % Varimax-rotate factor loadings         14
[RL,Y.G] = Varimax4M_progress(LU,iMax,nTol,Kaiser);       % Varimax-rotate factor loadings         14
EVr = sum(RL .* RL);                 % compute rotated Eigenvalues            15
[r, rx] = sort(EVr);                 % sort rotated Eigenvalues, keep indices 16
r = fliplr(r); rx = fliplr(rx);      % set descending order                   17
LR = RL(:,rx);                       % sort rotated factor loadings           18
s = ones(1,size(LR,2));              % current sign of loading vectors        19
s(abs(max(LR)) < abs(min(LR))) = -1; % determine direction of loading vectors 20
LR = LR .* repmat(s,vars,1);         % redirect loading vectors if necessary  21
tv = trace(EV);                      % compute total variance                 22
Y.VT = [u' 100*u'/tv ...             % table explained variance for unrotated 23
        r' 100*r'/tv ];              % ... and Varimax-rotated components     24
FSCFr = LR * inv(LR' * LR);          % compute rotated FS coefficients        25
FSCFr = FSCFr .* ...                 % rescale rotated FS coefficients by     26
   repmat(sigma,1,rk);               % ... the corresponding SDs              27
%   repmat(sqrt(diag(D)),1,rk);       % ... the corresponding SDs              27
if ~isempty(g.mu); 
 mu = g.mu; 
else
 mu = mean(X);                       % compute grand mean
end;
if ~isempty(g.sigma); 
 sigma = g.sigma;
else
 sigma = std(X);                     % compute SD
end;
%mu = mean(X); sigma = std(X);        % compute Mean and SD for each variable  28
%mu = mean(X); sigma = std(X,1);        % compute Mean and SD for each variable  28
%mu = mean(X); sigma = sigma';        % compute Mean and SD for each variable  28
%Xc = X - repmat(mu,cases,1);         % remove grand mean                      29
X = (X - repmat(mu,cases,1)) ./ repmat(sigma,cases,1);  % remove grand mean and normalize by SD
if fMax == 0
  for i = 1:size(Y.VT,1); 
    rk2 = i;                         % compute FSr for factors explaining <fPer>% total variance
    if sum(Y.VT(1:i,4)) >= fPer; break; end;
  end;
elseif fMax < 1
  rk2 = 0;  
  for i = 1:size(Y.VT,1); 
    if Y.VT(i,4) >= fPer; 
      rk2 = i;                       % compute FSr for factors explaining at least <fPer>% variance
    else break; end;
  end;
else
  rk2 = min(rk,fMax);                % compute FSr for requested number of factors
end
Y.LU = LU;
Y.LR = LR;
Y.FSr = zeros(cases,rk2);            % claim memory to speed computations     30
Y.mu = mu;
Y.sigma = sigma;
Y.FSCFr = FSCFr;
Y.rank = rk;
fprintf('Factor scores [');          % intialize progress bar
%Xc = Xc ./ repmat(sigma,cases,1);
%for m = 1:vars; Xc(:,m) = Xc(:,m) / sigma(m); end;
FSCFr = FSCFr';
for n = 1:cases;
 if mod(n-1,fix(cases/50)+1) == 0; fprintf('*'); end; % show progress
 for m = 1:rk2;                      % compute rotated factor scores from     31
%  Y.FSr(n,m) = sum( Xc(n,:) .* FSCFr(m,:) );   
  Y.FSr(n,m) = sum( X(n,:) .* FSCFr(m,:) );   
%  FSr(n,m) = sum( (Xc(n,:) ./ ...    % ... the normalized raw data and        32
%       sigma) .* FSCFr(:,m)' );      % ... the corresponding rescaled         33
 end; 
end;                                 % ... factor score coefficients          34
fprintf('] %.3g s\n',toc(tstart));     % finalize progress bar
