function coaipca
% COAIPCA, Does Individual in Principal Component Analysis
%      Corneal Topography Analysis System
%   Creates a laege size window with one page PCA
%      Intended to be called by coaipm.m
%      Reads individual file from FEATMPATH and INFILESTR
%      list of files from LISTFILESTR
% Outputs:
%     Only graphics, in current Figure

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




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


%  Set up global variables
%
%  First some general globals
global R0 NMAX MATLABVER 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
%  MATLABVER - version of MATLAB, 4 or 5
%  COMAP - Color Map
global INFILESTR FEATMPATH LISTFILESTR LISTR ;
%  INFILESTR - string with input filename
%  FEATMPATH - path for storing Zernike feature matrices
%  LISTFILESTR - string with input name of file with list of filenames
%  LISTR - string for inputting list file
global IPCAHANDPOS IPCAHAND ;
        %  window positioning vector,
        %  handle for Static PCA Window





%  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 




%  Read in matrix of Zernike feature coefficients for Individual
%
fvfilestr = [FEATMPATH INFILESTR '.mat'] ;
eval(['load ' fvfilestr]) ;
          %  coefficient are now in vector betahat



%  Read in matrix of Zernike feature coefficients for Pop'n
%
mfv = coaflr ;
d = size(mfv,1) ;
n = size(mfv,2) ;



%  Do Principal Components
%
sst = sum(sum(mfv.^2)) ;
%  get mean feature vector, and resids from that
vmean = mean(mfv')' ;
vsd = std(mfv')' ;
ssm = n * sum(vmean.^2) ;
ssmr = sst - ssm ;
mresid = mfv - vec2mat(vmean,n) ;

%  Get cov's and eigen values
mcov = cov(mresid') ;
          %  Get covariance matrix, transpose, since want "coordinates
          %                      as variables"

[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 max(abs())
mpc1 = ev1 * pc1rep ;
          %  "average vector" for display
ipc1 = (betahat - vmean)' * ev1 ;
vipc1 = ev1 * ipc1 ;
          %  projection of individ. onto ev1
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 max(abs())
mpc2 = ev2 * pc2rep ;
          %  "average vector" for display
ipc2 = (betahat - vmean)' * ev2 ;
vipc2 = ev2 * ipc2 ;
          %  projection of individ. onto ev2
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 max(abs())
mpc3 = ev3 * pc3rep ;
          %  "average vector" for display
ipc3 = (betahat - vmean)' * ev3 ;
vipc3 = ev3 * ipc3 ;
          %  projection of individ. onto ev3
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 max(abs())
mpc4 = ev4 * pc4rep ;
          %  "average vector" for display
ipc4 = (betahat - vmean)' * ev4 ;
vipc4 = ev4 * ipc4 ;
          %  projection of individ. onto ev4
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 ;

%  also add mean to individuals
vipc1 = vmean + vipc1 ;
vipc2 = vmean + vipc2 ;
vipc3 = vmean + vipc3 ;
vipc4 = vmean + vipc4 ;



%  summarize power percentages
%
vpopp = 100 * [ss1, ss2, ss3, ss4] / ssmr ;
          %  vector of powers (about the mean), for the population,
          %       as percentages
ssimr = sum((betahat - vmean).^2) ;
          %  Total sum of squares, for individual, about the mean
vindp = 100 * [ipc1^2, ipc2^2, ipc3^2, ipc4^2] / ssimr ;



%  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) ;
ivi = coazr(betahat) ;
ivipc1 = coazr(vipc1) ;
ivipc2 = coazr(vipc2) ;
ivipc3 = coazr(vipc3) ;
ivipc4 = coazr(vipc4) ;




%  Set up window page
%
IPCAHAND = figure ;
          %  create a new window
set(IPCAHAND, 'NumberTitle', 'off');
set(IPCAHAND,'Position', IPCAHANDPOS) ;
set(IPCAHAND,'Name','Individual in PCA') ;



%  Put down plots in 4 x 5 array
%
%  Mean in upper left
%
subplot(4,5,1) ;
pcolor(mxgrid,mygrid,imm) ;
  shading flat ;
  axis('square') ;
  title(['mean, ' LISTFILESTR]) ;
    caxis([0 1]) ;
    colormap(COMAP) ;


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

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


%  Now Full Recon. for Individual
%
subplot(4,5,6) ;
pcolor(mxgrid,mygrid,ivi) ;
  shading flat ;
  axis('square') ;
  title(['Full, ' INFILESTR]) ;
    caxis([0 1]) ;
  colormap(COMAP) ;


ivi = coazr(betahat) ;
ivipc1 = coazr(vipc1) ;


%  Now Projections for Individual
%
for i = 1:4 ;
  subplot(4,5,6+i) ;
  eval(['pcolor(mxgrid,mygrid,ivipc' num2str(i) ') ;']) ;
    shading flat ;
    axis('square') ;
    title(['PC' num2str(i) ', ' INFILESTR]) ;
      caxis([0 1]) ;
    colormap(COMAP) ;
end ;


%  Now do kde's
%
for i = 1:4 ;
  subplot(4,5,16+i) ;
  eval(['[kde, xgrid] = gpkde(pc' num2str(i)  ') ;']) ;
  plot(xgrid,kde,'w') ;
    title('Projections') ;
    vlim = axis ;
          %  Get axis values for drawing vertical lines
    eval(['right = pc' num2str(i) 'rep ;']) ;
    eval(['right = max(right,abs(ipc' num2str(i) ')) ;']) ;
          %  extend boundary if individual is outside
    left = - right ;
    range = right - left ;
    left = left - .05 * range ;
    right = right + .05 * range ;
    bottom = vlim(3) ;
    top = vlim(4) ;
    vcent = (bottom + top) / 2 ;
    axis([left right bottom top]) ;
    hold on ;
      evstr = ['plot([pc' num2str(i) ' pc' num2str(i) ...
                    '], [bottom vcent], ''c'') ;'] ;
      eval(evstr) ;
          %  Show locations of projections

      plot(xgrid,kde,'w') ;
          %  Do a replot to ensure seeing kde

      evstr = ['plot([pc' num2str(i) 'rep pc' num2str(i) ...
                    'rep], [bottom top], ''c--'') ;'] ;
      eval(evstr) ;
      evstr = ['plot([-pc' num2str(i) 'rep -pc' num2str(i) ...
                    'rep], [bottom top], ''c--'') ;'] ;
      eval(evstr) ;
          %  Show locations of representative guys

      evstr = ['plot([ipc' num2str(i) ' ipc' num2str(i) ...
                    '], [bottom top], ''r'') ;'] ;
      eval(evstr) ;
          %  Show location of individual

    hold off ;
end ;



%  Add a bar graph showing powers
%
subplot(4,5,16) ;
  if MATLABVER == 4 ;
    coabar(gca,vpopp,'c',100) ;
          %  put down a shaded bar graph for the population
    hold on ;
%    coabar(gca,vindp,'r',100) ;
      bar(vindp,'r') ;
          %  put down a hollow bar graph for the individual
          %  do this second, so everything shows up well
    hold off ;
  elseif MATLABVER == 5 ;
    bar(vpopp,'c') ;
         %  this use default bar width of 0.8
    hold on ;
      bar(vindp,.5,'r') ;
         %  this use smaller bar width of 0.5,
         %  so these can be contrasted with above
         %  (since they get overlaid)
    hold off ;
  end ;

  title('Relative SS''s') ;
  xlabel('PC') ;
  ylabel('%SS of Mean Res.') ;
    vax = axis ;
    tx = vax(1) + 0.1 * (vax(2) - vax(1)) ;
    ty1 = vax(3) + 0.9 * (vax(4) - vax(3)) ;
    ty2 = vax(3) + 0.8 * (vax(4) - vax(3)) ;
  text(tx,ty1,'Pop''n','Color','c') ;
  text(tx,ty2,'Indiv.','Color','r') ;



%  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,'XTickLabels',[]) ;
set(gca,'XTickMode','manual') ;
set(gca,'XTick',[]) ;
set(gca,'XColor',[0,0,0]) ;
set(gca,'YTickLabelMode','manual') ;
set(gca,'YTickLabels',[]) ;
set(gca,'YTickMode','manual') ;
set(gca,'YTick',[]) ;
set(gca,'YColor',[0,0,0]) ;
          %  Switch off axis stuff
titstr = ['Individual ' INFILESTR ' in PCA for file list: ' LISTFILESTR] ;
titstr = [titstr ',  ' date] ;
titstr = [titstr '  NMAX = ' num2str(NMAX) '  r0 = ' num2str(R0)] ;
text(0,0.5,titstr) ;
          %  put down actual title



