function coaspc2(iarg,rmp,rcp)
% COASPCA, Does Static Principal Component Analysis for 2 Populations
%      Corneal Topography Analysis System
%   Creates a large size window with one page PCA
%      Intended to be called by coa2pm.m
%      Reads list of files from LISTFILESTR1 and LISTFILESTR2
% Inputs:
%     iarg - Type of Analysis
%             1 - standard PCA 
%                    (eigen analysis of usual covariance matrix)
%             2 - Intra Class Spread normalized PCA
%                    (eigen analysis of cov matrix, modulo intra 
%                        class spread)
%             3 - Principle Discriminant Analysis
%                    (eigen analysis of inverse of within class cov
%                        times between class cov)
%     rmp - Robust Mean Parameter:
%             0 - (or not specified) use the ordinary sample mean
%             1 - Use the Spherical L1 location estimate
%             2 - Use the Elliptical (rescaled by MADs) L1 location estimate
%     rcp - Robust Covariance Parameter:
%             0 - (or not specified) use the ordinary covariance matrix
%             1 - Use the Spherical covariance matrix
%             2 - Use the Elliptical (rescaled by MADs) covariance matrix
% Outputs:
%     Only graphics, in current Figure

%    Copyright (c) J. S. Marron 1998
%    Department of Statistics, University of North Carolina



%  Set parameters and defaults according to number of input arguments
%
if nargin == 1 ;    %  1 argument input, use default
  irmp = 0 ;      %  default mean robustness parameter
else ;              %  was specified, use that
  irmp = rmp ;
end ;

if nargin <= 2 ;    %  at most 2 argument input, use default
  ircp = 0 ;      %  default covarince robustness parameter
else ;              %  was specified, use that
  ircp = rcp ;
end ;




%  Set local parameters
%
ngrid = 64 ;     %  Construct surfaces as ngrid x ngrid array
    


%  Set up global variables
%
%  First some general globals
global R0 NMAX IBASIS IRECENT COMAP ;
%  R0 - radius of analysis
          %    this gets reloaded when a feature vec is reloaded from file
%  NMAX - maximum number of Zernike columns
          %    this gets reloaded when a feature vec is reloaded from file
%  IBASIS - 1 for full (triangle), 2 for radial (3 row)
%  IRECENT - To recenter (or not) decompositions
%  COMAP - Color Map

global SPCAHANDPOS SPC2HAND FEATMPATH SUBDIRSTR LISTR ;
        %  window positioning vector,
        %  handle for Static PCA Window
        %  path for storing Zernike feature matrices
        %  SUBDIRSTR - current subdirectory (of both indata and outputs)
        %  string for inputting list file
global LISTFILESTR1 LISTFILESTR2 SETLIST1PATH SETLIST2PATH ;
        %  list file strings
        %  strings to set paths to lists
global PC2STR PC2LOAD1 PC2LOAD2 ;
        %  PC2STR - String with title info for 2 population analyses
        %  PC2LOAD1 - 1st Population Loadings for 2 population analyses
        %  PC2LOAD2 - 2nd Population Loadings for 2 population analyses



%  Read in matrix of Zernike feature coefficients, Population 1
%
disp('      Reading Files for Population 1') ;
eval(SETLIST1PATH) ;
          %  sets LISTR to Full Path List, for Population 1
mfv1 = coaflr ;
d1 = size(mfv1,1) ;
n1 = size(mfv1,2) ;



%  Read in matrix of Zernike feature coefficients, Population 2
%
disp('      Reading Files for Population 2') ;
eval(SETLIST2PATH) ;
          %  sets LISTR to Full Path List, for Population 2
mfv2 = coaflr ;
d2 = size(mfv2,1) ;
n2 = size(mfv2,2) ;




%  Calculate grid quantities
%
xgrid = linspace(-R0, R0, ngrid) ;
ygrid = xgrid ;	        %  Vectors of image coordinate values 
[mxgrid, mygrid] = meshgrid(xgrid,ygrid) ;
              		%  Matrices of image coordinate values 





%  Combine into one big data set
if d1 == d2 ;
  d = d1 ;
else ;
  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
  disp('!!!   Error from coaspc2.m:          !!!') ;
  disp('!!!   Different Number of Zernikes   !!!') ;
  disp('!!!      for these 2 lists           !!!') ;
  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
end ;
mfv = [mfv1, mfv2] ;
n = n1 + n2 ;




%  Main Principal Component Block (combined data set)
%
sst = sum(sum(mfv.^2)) ;

%  Get mean feature vector
%
if irmp == 0 ;        %  then do original non-robust sample mean
  vmean = mean(mfv')' ;
  meanstr = 'mean' ;
elseif irmp == 1 ;    %  then do spherical robust sample mean
  vmean = rmean(mfv')' ;
  meanstr = 'S. L1' ;
elseif irmp == 2 ;    %  then do elliptical robust sample mean
  vmad = gpmad(mfv',0)' ;
          %  column vector of MADs of each row
  smfv = mfv .* vec2mat(1 ./ vmad, n) ;
          %  make each coordinate have MAD 1
  vmean = rmean(smfv')' ;
  vmean = vmean .* vmad ;
  meanstr = 'E. L1' ;
end ;

%  Get residuals from mean
%
vsd = std(mfv')' ;
ssm = n * sum(vmean.^2) ;
ssmr = sst - ssm ;
mresid = mfv - vec2mat(vmean,n) ;

%  Get covariance matrix
%
if ircp == 0 ;        %  then do original non-robust covariance matrix
  mcov = cov(mresid') ;
          %  Get covariance matrix, transpose, since want "coordinates
          %                      as variables"
  covstr = '' ;
elseif ircp == 1 ;    %  then do spherical robust sample mean
  vrad = sqrt(sum(mresid .^2)) ;
          %  Transpose, since want "coordinates as variables"
          %  vector of radii of each curve
  sphereresid = mresid' ./ vec2mat(vrad',d) ;
          %  make each curve have "length" one
  mcov = cov(sphereresid) ;
          %  Get covariance matrix of "sphered" data
  covstr = 'S' ;
elseif ircp == 2 ;    %  then do elliptical robust sample mean
  vmad = gpmad(mresid',0)' ;
          %  column vector of MADs of each row
  smresid = mresid .* vec2mat(1 ./ vmad, n) ;
          %  make each coordinate have MAD 1
  vrad = sqrt(sum(smresid .^2)) ;
  sphereresid = smresid' ./ vec2mat(vrad',d) ;
          %  Transpose, since want "coordinates as variables"
          %  vector of radii of each curve
          %  make each curve have "length" one
  ellipsresid = sphereresid .* vec2mat(vmad', n) ;
          %  return each coordinate to original scale
  mcov = cov(ellipsresid) ;
          %  Get covariance matrix of "ellipsed" data
  covstr = 'E' ;
end ;



if iarg == 1 ;
  PC2STR = [covstr 'PCA'] ;
elseif iarg == 2 ;    % then normalize by "within class cov"
  PC2STR = [covstr 'ICS'] ;


  %  Get individual sample means
  %
  if irmp == 0 ;        %  then do original non-robust sample mean
    vmean1 = mean(mfv1')' ;
    vmean2 = mean(mfv2')' ;
  elseif irmp == 1 ;    %  then do spherical robust sample mean
    vmean1 = rmean(mfv1')' ;
    vmean2 = rmean(mfv2')' ;
  elseif irmp == 2 ;    %  then do elliptical robust sample mean
    vmad1 = gpmad(mfv1',0)' ;
    vmad2 = gpmad(mfv2',0)' ;
          %  column vector of MADs of each row
    smfv1 = mfv1 .* vec2mat(1 ./ vmad1, n1) ;
    smfv2 = mfv2 .* vec2mat(1 ./ vmad2, n2) ;
          %  make each coordinate have MAD 1
    vmean1 = rmean(smfv1')' ;
    vmean2 = rmean(smfv2')' ;
    vmean1 = vmean1 .* vmad1 ;
    vmean2 = vmean2 .* vmad2 ;
  end ;

elseif iarg == 3 ;    % then get principle discriminant "cov matrix"
  PC2STR = [covstr 'PDA'] ;

  %  Get individual sample means
  %
  if irmp == 0 ;        %  then do original non-robust sample mean
    vmean1 = mean(mfv1')' ;
    vmean2 = mean(mfv2')' ;
  elseif irmp == 1 ;    %  then do spherical robust sample mean
    vmean1 = rmean(mfv1')' ;
    vmean2 = rmean(mfv2')' ;
  elseif irmp == 2 ;    %  then do elliptical robust sample mean
    vmad1 = gpmad(mfv1',0)' ;
    vmad2 = gpmad(mfv2',0)' ;
          %  column vector of MADs of each row
    smfv1 = mfv1 .* vec2mat(1 ./ vmad1, n1) ;
    smfv2 = mfv2 .* vec2mat(1 ./ vmad2, n2) ;
          %  make each coordinate have MAD 1
    vmean1 = rmean(smfv1')' ;
    vmean2 = rmean(smfv2')' ;
    vmean1 = vmean1 .* vmad1 ;
    vmean2 = vmean2 .* vmad2 ;
  end ;

  mresid1 = mfv1 - vec2mat(vmean1,n1) ;
  mresid2 = mfv2 - vec2mat(vmean2,n2) ;

  %  Get individual covariance matrices
  %
  if ircp == 0 ;        %  then do original non-robust covariance matrix
    mcov1 = cov(mresid1') ;
    mcov2 = cov(mresid2') ;
          %  Get covariance matrix, transpose, since want "coordinates
          %                      as variables"
  elseif ircp == 1 ;    %  then do spherical robust sample mean
    vrad1 = sqrt(sum(mresid1 .^2)) ;
    vrad2 = sqrt(sum(mresid2 .^2)) ;
          %  Transpose, since want "coordinates as variables"
          %  vector of radii of each curve
    sphereresid1 = mresid1' ./ vec2mat(vrad1',d) ;
    sphereresid2 = mresid2' ./ vec2mat(vrad2',d) ;
          %  make each curve have "length" one
    mcov1 = cov(sphereresid1) ;
    mcov2 = cov(sphereresid2) ;
          %  Get covariance matrix of "sphered" data
  elseif ircp == 2 ;    %  then do elliptical robust sample mean
    vmad1 = gpmad(mresid1',0)' ;
    vmad2 = gpmad(mresid2',0)' ;
          %  column vector of MADs of each row
    smresid1 = mresid1 .* vec2mat(1 ./ vmad1, n1) ;
    smresid2 = mresid2 .* vec2mat(1 ./ vmad2, n2) ;
          %  make each coordinate have MAD 1
    vrad1 = sqrt(sum(smresid1 .^2)) ;
    vrad2 = sqrt(sum(smresid2 .^2)) ;
    sphereresid1 = smresid1' ./ vec2mat(vrad1',d) ;
    sphereresid2 = smresid2' ./ vec2mat(vrad2',d) ;
          %  Transpose, since want "coordinates as variables"
          %  vector of radii of each curve
          %  make each curve have "length" one
    ellipsresid1 = sphereresid1 .* vec2mat(vmad1', n1) ;
    ellipsresid2 = sphereresid2 .* vec2mat(vmad2', n2) ;
          %  return each coordinate to original scale
    mcov1 = cov(ellipsresid1) ;
    mcov2 = cov(ellipsresid2) ;
          %  Get covariance matrix of "ellipsed" data
  end ;


  mcovb = vmean1 - vmean2 ;
  mcovb = mcovb * mcovb' ;

  mcovw = (n1 / n) * mcov1 + (n2 / n) * mcov2 ;
            %  Average of "within class covariances"

  mcov = pinv(mcovw) ;
  mcov = (mcov + mcov') / 2 ;
            %  symmetrize, since psuedo inv is not quite symmetric
%max(max(abs(mcov - mcov')))
%max(max(abs(mcovb - mcovb')))

  mcov = pinv(mcovw) * mcovb ;
          %  original line, appears to be numerically unstable,
          %  in particular result is not symmetric
  mcov = (mcov + mcov') / 2 ;
            %  symmetrize, since not quite symmetric

%  mcov = sqrtm(mcov) ;
%  mcov = mcov * mcovb * mcov ;
          %  a symmetrized version


%max(max(abs(mcov - mcov')))
%mcov(1:16,1:16)

end ;

% Do eigen analysis
%
[meigvec,meigval] = eig(mcov) ;
veigval = (diag(meigval)) ;
          %  Get eigenvalues and eigenvectors

[veigval,vind] = sort(veigval) ;
meigvec = meigvec(:,vind) ;
          %  Sort to avoid possible ordering problems
  
ev1 = meigvec(:,d) ;
          %  eigenvector with largest eigenvalue
pc1 = mresid' * ev1 ;
          %  vector of projections of data onto ev1
pc1rep = max(abs(pc1 - mean(pc1))) ;
          %  "most representative pc", defined as MAD
mpc1 = ev1 * pc1rep ;
          %  "average vector" for display
mresid1 = (mresid - ev1 * pc1') ;
ss1r = sum(sum(mresid1.^2)) ;
ss1 = ssmr - ss1r ;

ev2 = meigvec(:,d-1) ;
          %  eigenvector with largest eigenvalue
pc2 = mresid1' * ev2 ;
          %  vector of projections of data onto ev2
pc2rep = max(abs(pc2 - mean(pc2))) ;
          %  "most representative pc", defined as MAD
mpc2 = ev2 * pc2rep ;
          %  "average vector" for display
mresid2 = (mresid1 - ev2 * pc2') ;
ss2r = sum(sum(mresid2.^2)) ;
ss2 = ss1r - ss2r ;

ev3 = meigvec(:,d-2) ;
          %  eigenvector with largest eigenvalue
pc3 = mresid2' * ev3 ;
          %  vector of projections of data onto ev3
pc3rep = max(abs(pc3 - mean(pc3))) ;
          %  "most representative pc", defined as MAD

mpc3 = ev3 * pc3rep ;
          %  "average vector" for display
mresid3 = (mresid2 - ev3 * pc3') ;
ss3r = sum(sum(mresid3.^2)) ;
ss3 = ss2r - ss3r ;

ev4 = meigvec(:,d-3) ;
          %  eigenvector with largest eigenvalue
pc4 = mresid3' * ev4 ;
          %  vector of projections of data onto ev4
pc4rep = max(abs(pc4 - mean(pc4))) ;
          %  "most representative pc", defined as MAD
mpc4 = ev4 * pc4rep ;
          %  "average vector" for display
mresid4 = (mresid3 - ev4 * pc4') ;
ss4r = sum(sum(mresid4.^2)) ;
ss4 = ss3r - ss4r ;



%  Form mean +- representative eigenvectors
mpc1p = vmean + mpc1 ;
mpc1m = vmean - mpc1 ;
mpc2p = vmean + mpc2 ;
mpc2m = vmean - mpc2 ;
mpc3p = vmean + mpc3 ;
mpc3m = vmean - mpc3 ;
mpc4p = vmean + mpc4 ;
mpc4m = vmean - mpc4 ;




%  Create images
%
imm = coazr(vmean) ;
impc1p = coazr(mpc1p) ;
impc1m = coazr(mpc1m) ;
impc2p = coazr(mpc2p) ;
impc2m = coazr(mpc2m) ;
impc3p = coazr(mpc3p) ;
impc3m = coazr(mpc3m) ;
impc4p = coazr(mpc4p) ;
impc4m = coazr(mpc4m) ;




%  Set up window page
%
SPC2HAND = figure ;
          %  create a new window
set(SPC2HAND, 'NumberTitle', 'off');
set(SPC2HAND,'Position', SPCAHANDPOS) ;
set(SPC2HAND,'Name',['Static ' PC2STR]) ;


%  Put down plots in 4 x 4 array
%
%  Means first
%
for i = 1:4 ;
  subplot(4,4,4+i) ;
  pcolor(mxgrid,mygrid,imm) ;
    shading flat ;
    axis('square') ;
    title([meanstr ', ' LISTFILESTR1 ', ' LISTFILESTR2]) ;
      caxis([0 1]) ;
    colormap(COMAP) ;
  hold on ;
    if IRECENT == 0 ;  %  then everything is centered at apex normals
      plot(0,0,'+k');
    elseif IRECENT == 1 ;  %  then everything is recentered
      plot(0,0,'ok');
    end ;
    tx = -.95 * R0 ;
    ty = .7 * R0 ;
    eval(['pcss = 100 * ss' num2str(i) ' / ssmr ;']) ;
    ssstr = num2str(pcss) ;
    tstr = [ssstr '% of MR'] ;
    text(tx,ty,tstr) ;
  hold off ;
end ;


%  Now mean +- Pc's
%
for i = 1:4 ;
  subplot(4,4,i) ;
  eval(['pcolor(mxgrid,mygrid,impc' num2str(i) 'p)']) ;
    shading flat ;
    axis('square') ;
    title([meanstr ' + ' PC2STR num2str(i)]) ;
      caxis([0 1]) ;
    colormap(COMAP) ;

  subplot(4,4,8+i) ;
  eval(['pcolor(mxgrid,mygrid,impc' num2str(i) 'm)']) ;
    shading flat ;
    axis('square') ;
    title([meanstr ' - ' PC2STR num2str(i)]) ;
      caxis([0 1]) ;
    colormap(COMAP) ;
end ;


%  Now do kde's
%
for i = 1:4 ;
  subplot(4,4,12+i) ;
  eval(['pccur = pc' num2str(i) ';']) ;
          %  Matrix of current PC loadings
  pccur1 = pccur(1:n1) ;
          %  Current PC Loadings, for Population 1
  pccur2 = pccur(n1+1:n) ;
          %  Current PC Loadings, for Population 2
  [kde1, xgrid1] = gpkde(pccur1) ;
  [kde2, xgrid2] = gpkde(pccur2) ;
  plot(xgrid1,kde1,'r',xgrid2,kde2,'c') ;
    title('Projections') ;
    vlim = axis ;
          %  Get axis values for drawing vertical lines
    eval(['left = -pc' num2str(i) 'rep ;']) ;
    eval(['right = pc' num2str(i) 'rep ;']) ;
    range = right - left ;
    left = left - .05 * range ;
    right = right + .05 * range ;
    bottom = vlim(3) ;
    top = vlim(4) ;
    axis([left right bottom top]) ;
    hold on ;
      plot([pccur2 pccur2],[bottom top], 'c') ;
      plot([pccur1 pccur1],[bottom top], 'r') ;
          %  Show locations of projections 
          %  (pop'n 2 first, since likely more spread)
      plot(xgrid1,kde1,'r',xgrid2,kde2,'c') ;
          %  Do a replot to ensure seeing kde
    hold off ;
end ;




%  Add a title
%
tx = 0.05 ;
ty = 0.94 ;
tw = 0.8 ;
th = 0.05 ;
          %  set position for new axes, containing only title
axes('Position',[tx,ty,tw,th]) ;
          %  create new axes
set(gca,'XTickLabelMode','manual') ;
set(gca,'XTickLabel',[]) ;
set(gca,'XTickMode','manual') ;
set(gca,'XTick',[]) ;
set(gca,'XColor',[0,0,0]) ;
set(gca,'YTickLabelMode','manual') ;
set(gca,'YTickLabel',[]) ;
set(gca,'YTickMode','manual') ;
set(gca,'YTick',[]) ;
set(gca,'YColor',[0,0,0]) ;
          %  Switch off axis stuff
set(gca,'visible','off') ;
          %  Turn off surrounding box
titstr = ['2 Population Static (' meanstr ')' PC2STR ' for file lists: ' ...
               LISTFILESTR1 LISTFILESTR2] ;
titstr = [titstr ',  ' date] ;
titstr = [titstr '  NMAX = ' num2str(NMAX) '  r0 = ' num2str(R0)] ;
text(0,0.5,titstr) ;
          %  put down actual title



%  Save 2 population analysis Loadings
%
mload = [pc1 pc2 pc3 pc4] ;
PC2LOAD1 = mload(1:n1,:) ;
PC2LOAD2 = mload(n1+1:n,:) ;

