disp('Running MATLAB script file vire4.m') ;
%

%    Reads already stored Zernike Decomps, 
%    and does PCA on set of all individuals
%



ipart = 2 ;      %  1 - read through data, construct matrix
                 %  2 - Do PCA on individuals



iR0 = 2 ;        %  1 - R0 = 4.2
                 %  2 - R0 = 3.0


ips = 1 ;
ipe = 4 ;
          %  loop through all patients


if iR0 == 1 ;
  outsuf = 'a' ;
  R0 = 4.2 ;
elseif iR0 == 2 ;
  outsuf = 'b' ;  
  R0 = 3.0 ;
end ;


ngrid = 64 ;     %  Construct surfaces as ngrid x ngrid array
xgrid = linspace(-R0, R0, ngrid) ;
ygrid = xgrid ;	        %  Vectors of image coordinate values 

global R0 NMAX IBASIS IRECENT INDATPATH HICURVE ;
NMAX = 10 ;
IBASIS = 1 ;
IRECENT = 0 ;


INDATPATH = '\matlab\cornean\indata\' ;
          %  path for input Keratron files
[COMAP, HICURVE] = coacmap ;




if ipart == 1 ;    %  then do PCA on set of individuals


  clf ;

  %  Do appropriate file reads
  %
  mout = [] ;
  for ip = ips:ipe ;

    if ip == 1 ;    
      patname = 'Hickson' ;
      featmpath = '\matlab\cornean\zfm\virehick\' ;
      listr = [featmpath 'hickall.lst'] ;
    elseif ip == 2 ;    
      patname = 'Meadows' ;
      featmpath = '\matlab\cornean\zfm\viremead\' ;
      listr = [featmpath 'meadall.lst'] ;
    elseif ip == 3 ;    
      patname = 'Tesoro' ;
      featmpath = '\matlab\cornean\zfm\vireteso\' ;
      listr = [featmpath 'tesoall.lst'] ;
    elseif ip == 4 ;    
      patname = 'Yi' ;
      featmpath = '\matlab\cornean\zfm\vireyi\' ;
      listr = [featmpath 'yiall.lst'] ;
    end ;



    %  read files, using lines from   coalwzr.m
    %
    fid = fopen(listr) ;
            %  assigns a file handle, (read only assumed)

    if (fid == -1) ;      %  then file not found, so report
      disp('   ') ;
      disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
      disp('!!!   Warning from vire1.m:   can''t find file:') ;
      disp(['!!!       ' listr]) ;
      disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
    end ;


    %  Read in list of filenames
    %
    minfilestr = fscanf(fid,'%s',[8,inf]) ;
              %  formatted read, assumes single column
    fclose(fid) ;

    minfilestr = minfilestr' ;
            %  each row corresponds to one filename


    %  loop through files (one for each image)
    %
    ie = size(minfilestr,1) ;
    for i = 1:ie ;

      infilestr = minfilestr(i,:) ;

      disp(['  Loading ' num2str(i) ' of ' num2str(ie) ...
                               ' for ' patname ', file: ' infilestr ]) ;


      zmistr = [featmpath infilestr '.mat'] ;
            %  full path and filename 
      eval(['load ' zmistr]) ;
            %  loads this Zernike feature matrix and parameter values

      vout = betahat' ;
            %  turn column vector of coefficients into a row vector

%      filenum = str2num(infilestr(1:6)) ;
%      filenum = filenum + (str2num(infilestr(8))) / 10 ;
            %  get first 6 numbers as an integer, 
            %  then add last digit after decimal

%      vout = [filenum ip vout] ;
            %  insert ASCII numerical file number, and patient number


      mout = [mout; vout] ;
            %  tack on new row to output matrix (this patient)


    end ;


  end ;


    sfilestr = ['\matlab\cornean\vr4a' outsuf] ;
  eval(['save ' sfilestr ' mout']) ;
          %  Save as .mat file  



elseif ipart == 2 ;    % then do PCA 

  clf ;

    sfilestr = ['\matlab\cornean\vr4a' outsuf] ;
  eval(['load ' sfilestr]) ;
          %  Load mout Matrix from above


  %  Do PCA, using lines from coaspca.m
  %


  %  Set local parameters
  %
  ngrid = 64 ;     %  Construct surfaces as ngrid x ngrid array
  irmp = 0 ;    %  0 - center at mean    
                %  1 - spherical mean
                %  2 - elliptical mean
  ircp = 0 ;    %  0 - usual covariance
                %  1 - spherical covariance
                %  2 - elliptical covariance



%  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





  mfv = mout' ;
          %  mfv wants data vectors as columns
  d = size(mfv,1) ;
  n = size(mfv,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 



  %  Main Principal Component Block
  %
  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 = 'PC' ;
  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 = 'SPC' ;
  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 = 'EPC' ;
  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) ;






  %  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 ', All']) ;
        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 ' + ' covstr 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 ' - ' covstr num2str(i)]) ;
        caxis([0 1]) ;
      colormap(COMAP) ;
  end ;


  %  Now do kde's
  %
  for i = 1:4 ;
    subplot(4,4,12+i) ;
    eval(['[kde, xgrid] = gpkde(pc' num2str(i)  ') ;']) ;
    plot(xgrid,kde,'w') ;
      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 ;
        evstr = ['plot([pc' num2str(i) ' pc' num2str(i) ...
                      '], [bottom top], ''c'') ;'] ;
        eval(evstr) ;
          %  Show locations of projections
        plot(xgrid,kde,'w') ;
          %  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 = ['Static PCA, using ' meanstr ' and ' covstr ...
                         'A, for All'] ;
  titstr = [titstr ',  ' date] ;
  titstr = [titstr '  NMAX = ' num2str(NMAX) '  r0 = ' num2str(R0)] ;
  text(0,0.5,titstr) ;
          %  put down actual title



  %  Print resulting PCA Plot
  %
  orient landscape ;
  eval(['print -dpsc \matlab\cornean\vr4b' outsuf '.ps']) ;


  %  Write eigenvectors into a file
  %
  mout = meigvec(:,(d-3):d) ;
        %  4 eigenvectors with largest eigenvalues
  mout = [(1:size(mout,1))', mout ] ;
        %  Add column of Zernike indices

  fid = fopen(['\matlab\cornean\vr4b' outsuf '.out'],'wt') ;

    titlstr1 = 'Output from the MATLAB Script vire4.m,   ' ;
    titlstr1 = [titlstr1,date] ;
  cntbytes = fprintf(fid,'%1s\n\n\n',titlstr1) ;
          %  '%1s/n' says "print string, followed by a line feed"
          %       (Note, this only goes to next line)

    titlstr2 = ['Zernike #     Weights 1     Weights 2' ...
                     '     Weights 3     Weights 3'] ;
  cntbytes = fprintf(fid,'%1s\n\n',titlstr2) ;

    forstr = '      %3.0f   %11.8f   %11.8f   %11.8f   %11.8f\n' ;
  cntbytes = fprintf(fid,forstr,mout') ;
          %  Careful here, because this reads down columns first,
          %  So feed it the transpose of the usual form

  fclose(fid) ;




end ;


