function [veigval,meigvec] = pcaSM(mdat,paramstruct) 
% PCASM, Principal Component Analysis
%     Computes eigenvalues and eigenvectors for data vectors
%     in mdat.  Eigenvectors are orthogonal directions of
%     greatest variability, and eigenvalues are sums of
%     squares of projections of data in those directions
%   Steve Marron's matlab function
%     
% Inputs:
%
%   mdat        - d x n matrix of multivariate data
%                         (each col is a data vector)
%                         d = dimension of each data vector
%                         n = number of data vectors 
%
%   paramstruct - a Matlab structure of input parameters
%                    Use: "help struct" and "help datatypes" to
%                         learn about these.
%                    Create one, using commands of the form:
%
%       paramstruct = struct('field1',values1,...
%                            'field2',values2,...
%                                             ) ;
%
%                          where any of the following can be used,
%                          these are optional, misspecified values
%                          revert to defaults
%
%    fields            values
%
%    npc              number of principal directions to compute
%                     0  (default)  compute "full set" = rank(mdat)
%                            this cannot be bigger than rank of mdat
%                            when doing dual eigen decomposition
%
%    idual            1  force use of dual eigen decomposition
%                     0  (default)  choose smallest
%                     -1  force use of direct eigen decomposition
%
%    iprestd          0  (default)  treat as raw data and standardize,
%                             i.e. subtract mean vector, 
%                                     and mult. by 1/sqrt(n-1)
%                     1  treat as pre-standardized data, i.e.
%                             leave mean and scale alone
%
%    iscreenwrite     0  (default)  no screen writes
%                     1  write to screen to show progress
%
% Outputs:
%
%    veigval  -  vector of 1st npc eigenvalues, sorted in decreasing order
%
%    meigvec  -  vector of corresponding npc eigenvectors
%
%          The covariance matrix of the data admits the decomposition:
%               covariance = meigvec * diag(veigval) * meigvec'
%               i.e. meigvec' * covariance * meigvec = diag(veigval)
%
% Assumes path can find personal function:
%    vec2matSM

%    Copyright (c) J. S. Marron 2001



%  First set all parameter to defaults
%
npc = 0 ;
inpcset = 0 ;
idual = 0 ;
iprestd = 0 ;
iscreenwrite = 0 ;



%  Now update parameters as specified,
%  by parameter structure (if it is used)
%
if nargin > 1 ;   %  then paramstruct is an argument

  if isfield(paramstruct,'npc') ;    %  then change to input value
    npc = getfield(paramstruct,'npc') ; 
    inpcset = 1 ;
          %  record fact that npc was set 
  end ;

  if isfield(paramstruct,'idual') ;    %  then change to input value
    idual = getfield(paramstruct,'idual') ; 
  end ;

  if isfield(paramstruct,'iprestd') ;    %  then change to input value
    iprestd = getfield(paramstruct,'iprestd') ; 
  end ;

  if isfield(paramstruct,'iscreenwrite') ;    %  then change to input value
    iscreenwrite = getfield(paramstruct,'iscreenwrite') ; 
  end ;


end ;  %  of resetting of input parameters



%  get dimensions from data matrix
%
d = size(mdat,1) ;
          %  dimension of each data vector (col)
n = size(mdat,2) ;
          %  number of data vectors (cols)


if iprestd ~= 1 ;    %  then need to standardize input data matrix

  vmean = mean(mdat,2) ;
          %  mean across rows

  mdat = mdat - vec2matSM(vmean,n) ;
          %  recenter data point cloud

  mdat = (1 / sqrt(n - 1)) * mdat ;
          %  puts on "covariance scale", so that covariance matrix
          %  is just "outer product":  mdat * mdat'

end ;

rankmdat = rank(mdat) ;



if idual == 1 ;    %  go with dual computation
  iidual = 1 ;
elseif idual == -1 ;    %  go with direct computation
  iidual = -1 ;
else ;
  if d > n ;
    iidual = 1 ;
          %  dual computation is simplest
  else ;
    iidual = -1 ;
          %  direct computation is simplest
  end ;
end ;



%  set additional parameters
%
eigsopts = struct('disp', 0) ;
          %  Option structure for eigs, this supresses output

if npc == 0 ;    %  Then compute non-zero eigenvalues
  npc = rankmdat ;
end ;



if iscreenwrite == 1 ;
  disp('    pcaSM starting eigendecomposition') ;
end ;



if iidual == -1 ;    %  do direct eigen decomposition

  if npc > d ;
    npc = d ;
    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
    disp('!!!   Warning from pcaSM.m:      !!!') ;
    disp('!!!   can''t do dual eignvalue    !!!') ;
    disp('!!!   computation, since rank    !!!') ;
    disp('!!!   of data matrix too small   !!!') ;
    disp(['!!!   reducing npc to d =  ' num2str(npc) '     !!!']) ;
    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;

  end ;

  sigmahat = mdat * mdat' ;
          %  direct covariance is outer product of normalized data matrix

  if npc < d ;    %  do reduced eigen decomposition
  
    [B,D] = eigs(sigmahat,npc,'LA',eigsopts) ;
          %  partial eigenvalue decomp., using only npc largest
          %  'LR' for Largest Real part

    veigval = diag(D) ;
    meigvec = B ;

  else ;    %  do full eigendecomposition

    [B,D] = eig(sigmahat) ;

    veigval = diag(D) ;
    meigvec = B ;

    %  careful need to do a sort
    %
    [veigval,sortind] = sort(veigval) ;
        %  sorted in increasing order
    veigval = flipud(veigval) ;
    sortind = flipud(sortind) ;
        %  sorted in decreasing order
    meigvec = meigvec(:,sortind) ;
        %  put in same order

  end ;


elseif iidual == 1 ;    %  do dual eigen decomposition

  if npc > rankmdat ;
    npc = rankmdat ;
    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
    disp('!!!   Warning from pcaSM.m:      !!!') ;
    disp('!!!   can''t do dual eignvalue    !!!') ;
    disp('!!!   computation, since rank    !!!') ;
    disp('!!!   of data matrix too small   !!!') ;
    disp(['!!!   reducing npc to: ' num2str(npc) '        !!!']) ;
    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;

  end ;

  sigmastar = mdat' * mdat ;
          %  dual covariance is inner product of normalized data matrix

  [Bstar,Dstar] = eigs(sigmastar,npc,'LA',eigsopts) ;
          %  partial eigenvalue decomp., using only npc largest
          %  'LR' for Largest Real part


  veigval = diag(Dstar) ;
  irDstar = diag(1 ./ sqrt(veigval)) ;
  meigvec = mdat * Bstar * irDstar ;


end ;



if iscreenwrite == 1 ;
  disp('    pcaSM finished with pca calculation') ;
  disp('  ') ;
end ;



