function outstruct = 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
% (otherwise will reset to rank(mdat))
%
% 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
% (note assumes 1/sqrt(n-1)
% normalization has been done)
%
% iscreenwrite 0 (default) no screen writes
% 1 write to screen to show progress
%
% viout vector of indices 0 or 1 for various types of output
% 1 in a given entry will include that
% in output structure
% Entry 1: veigval - npc x 1 vector of 1st npc
% eigenvalues, sorted in decreasing order
% Need to multiply by (n-1) for these to
% work like sums of squares
% Entry 2: meigvec - d x npc matrix of corresponding
% eigenvectors
% These are orthonormal
% Entries are typically called "loadings"
% Entry 3: vmean - d x 1 mean vector
% (returns assumed value of zeros(1,n)
% for iprestd = 1)
% Entry 4: mmeanresid - d x n matrix of
% residuals from mean
% (returns assumed value of
% mdat * sqrt(n-1), for iprestd = 1)
% Entry 5: mpc - npc x n matrix of
% "principal components"
% i.e. of coefficients of projections
% of data onto eigenvectors,
% i.e. of "scores".
% Note can "recover data" from
% vmean + meigvec * mpc
% Common "PCi" projection plot comes
% from meigvec(:,i) * mpc(i,:)
% Entry 6: a3proj - d x n x npc 3d array of
% projections of data onto individual
% eigenvectors.
% Use this for PCi projection plots
% Entry 7: sstot - 1 x 1 scalar, with
% overall Sum of Squares of mdat
% Entry 8: ssmresid - 1 x 1 scalar, with
% mean residual Sum of Squares
% (divide by sstot, and multiply
% by 100, to get percent of sstot)
% Entry 9: vpropSSmr - npc x 1 vector of
% proportions of SS,
% relative to meanresid SS
% (multiply by 100 to get percent)
% (use cumsum, for "cumulative power")
% Entry 10: vpropSSpr - npc x 1 vector of
% proportions of SS,
% relative to Previous Residual SS
% (relative to mean resid, for 1st)
% (multiply by 100 to get percent)
% When viout has length < 10, will pad with 0's
% Default is [1, 1]
%
% Output:
%
% outstruct - structure with outputs, as requested in viout above
%
% Note: unpack these with commands of the form:
% veigval = getfield(outstruct,'veigval') ;
%
% 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 ;
viout = [1, 1] ;

% 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') ;
% 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 ;

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

end ; % of resetting of input parameters

% Readjust viout as needed
%
maxviout = 10 ;
% largest useful size for viout
if size(viout,1) > 1 ; % if have more than one row
viout = viout' ;
end ;
if size(viout,1) == 1 ; % then have row vector

if length(viout) < maxviout ; % then pad with 0s
viout = [viout zeros(1,maxviout - length(viout))] ;
end ;

else ; % invalid viout, so indicate and quit

disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
disp('!!! Error from pcaSM: !!!') ;
disp('!!! Invalid viout !!!') ;
disp('!!! Terminating exceution !!!') ;
disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
return ;

end ;

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

% Calculate total sum of squares
%
if viout(7) == 1 ;

sstot = sum(sum(mdat.^2)) ;

end ;

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

if viout(4) == 1 | ...
viout(8) == 1 | ...
viout(9) == 1 | ...
viout(10) == 1 ; % need to output mean residuals
mmeanresid = mdat ;
end ;

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

else ;

if viout(3) == 1 ; % need to output actual mean (despite assumption)
vmean = zeros(d,1) ;
% assumed mean vector
end ;

if viout(4) == 1 | ...
viout(8) == 1 | ...
viout(9) == 1 | ...
viout(10) == 1 ; % need to output actual mmeanresid
% (despite assumption)
mmeanresid = mdat * sqrt(n - 1) ;
% assumed mean residuals
end ;

if iscreenwrite == 1 ;
disp('!!! pcaSM assuming that data are standardized !!!') ;
end ;

end ;

rankmdat = rank(mdat) ;

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

if npc > rankmdat ;
disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
disp('!!! Warning from pcaSM: npc too large !!!') ;
disp('!!! resetting to rank of mdat !!!') ;
disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
npc = rankmdat ;
end ;

if (idual ~= 1) & (idual ~= -1) ; % then need to use inputs
% to choose idual
if d > n ;
idual = 1 ;
% dual computation is simplest
else ;
idual = -1 ;
% direct computation is simplest
end ;

end ;

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

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

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

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 idual == 1 ; % do dual eigen decomposition

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 viout(5) == 1 ;
mpc = sqrt(n - 1) * meigvec' * mdat ;
% npc x n matrix of projections of data on eigenvectors
end ;

if viout(6) == 1 ;

% create 3 dim array versions
%
a3meigvec = [] ;
for i = 1:n ;
a3meigvec = cat(3,a3meigvec,meigvec) ;
end ;
a3meigvec = permute(a3meigvec,[1 3 2]) ;
% essentially "transpose" in 3d
% need to swap 2nd and 3rd dimensions
% to have give d x n x npc result

a3mpc = [] ;
for i = 1:d ;
a3mpc = cat(3,a3mpc,mpc) ;
end ;
a3mpc = permute(a3mpc,[3 2 1]) ;
% essentially "transpose" in 3d
% need to swap 1st and 3rd dimensions
% to have give d x n x npc result

a3proj = times(a3meigvec,a3mpc) ;
% high d version of usual .* operation, gives:
% d x n x npc 3d array of projections of data
% onto individual eigenvectors.

end ;

% Calculate needed sum of squares
%
if viout(8) == 1 | ...
viout(9) == 1 | ...
viout(10) == 1 ;

ssmresid = sum(sum(mmeanresid.^2)) ;

end ;

if viout(9) == 1 | ...
viout(10) == 1 ;

vpcSS = (n - 1) * veigval ;
% sum of squares of pc's in each eigen-direction
vpropSSmr = vpcSS / ssmresid ;

end ;

if viout(10) == 1 ;

prevresidSS = ssmresid ;
% previous residual Sum of Squares
residSS = ssmresid - vpcSS(1) ;
% current residual Sum of Squares
vpropSSpr = residSS / prevresidSS ;

for iev = 2:npc ;
prevresidSS = residSS ;
residSS = prevresidSS - vpcSS(iev) ;
if residSS < (10^(-12) * ssmresid) ; % then any difference with 0
% is regarded as round off errors
% so set to 0, for output
residSS = 0 ;
end ;
vpropSSpr = [vpropSSpr; (residSS / prevresidSS)] ;

end ;

end ;

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

% Construct output structure
%
outstruct = struct('dummy',[]) ;
% put in dummy field, to avoid empty structure

if viout(1) == 1 ;
outstruct = setfield(outstruct,'veigval',veigval) ;
end ;

if viout(2) == 1 ;
outstruct = setfield(outstruct,'meigvec',meigvec) ;
end ;

if viout(3) == 1 ;
outstruct = setfield(outstruct,'vmean',vmean) ;
end ;

if viout(4) == 1 ;
outstruct = setfield(outstruct,'mmeanresid',mmeanresid) ;
end ;

if viout(5) == 1 ;
outstruct = setfield(outstruct,'mpc',mpc) ;
end ;

if viout(6) == 1 ;
outstruct = setfield(outstruct,'a3proj',a3proj) ;
end ;

if viout(7) == 1 ;
outstruct = setfield(outstruct,'sstot',sstot) ;
end ;

if viout(8) == 1 ;
outstruct = setfield(outstruct,'ssmresid',ssmresid) ;
end ;

if viout(9) == 1 ;
outstruct = setfield(outstruct,'vpropSSmr',vpropSSmr) ;
end ;

if viout(10) == 1 ;
outstruct = setfield(outstruct,'vpropSSpr',vpropSSpr) ;
end ;

outstruct = rmfield(outstruct,'dummy') ;