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

%    Reads already stored Zernike Decomps, does PCA
%    and saves as ASCII files, for e.g. SAS analysis
%



ipart = 1 ;      %  1 - Create Ascii files
                 %  2 - Some Sum of Square Analyses

ipat = 0 ;       %  0 - loop through all
                 %  1 - Hickson
                 %  2 - Meadows
                 %  3 - Tesoro 
                 %  4 - Yi     


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



npc = 20 ;       %  Number of Principal Components to Keep


irmp = 0 ;       %  0 - ordinary sample mean
                 %  1 - Spherical L1 Mean
                 %  2 - Elliptical L1 Mean

ircp = 0 ;       %  0 - ordinary covraince PCA
                 %  1 - Spherical PCA
                 %  2 - Elliptical PCA


%  Initialize parameters
%
if ipat == 0 ;     %  then loop through all patients
  ips = 1 ;
  ipe = 4 ;
else ;             %  then just do input patients
  ips = ipat ;
  ipe = ipat ;
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 ;

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

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





if ipart == 1 ;    %  then create big ASCII file with results


  %  Do appropriate file reads
  %
  mbetahat = [] ;
  vfilenum = [] ;
  vip = [] ;
  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 vire2.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

      mbetahat = [mbetahat, betahat] ;



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

      vip = [vip; ip] ;



    end ;

  end ;




  %  Main Principal Component Block
  %
  d = size(mbetahat,1) ;
  n = size(mbetahat,2) ;

  sst = sum(sum(mbetahat.^2)) ;

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

  %  Get residuals from mean
  %
  mresid = mbetahat - 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
  



  mev = meigvec(:,(d - npc + 1):d) ;
          %  matrix of eigenvectors with largest eigenvalues
  mpc = mresid' * mev ;
          %  matrix of projections of data onto 
          %  eigenvectors with largest eigenvalues
  mpc = mpc(:,npc:-1:1) ;
          %  reverse the columns, to put the projections
          %  corresponding to the largest eigenvalues on the left


  %  Next read in patient information using lines from vire1r.m
  %
  fid = fopen(['\matlab\cornean\vire1' outsuf '.in'],'rt') ;
          %  'rt' is for "read only" and "text"

  forstr = '%f %f %f %f %f' ;
          %  format for reading:
          %      file no.
          %      pat. no.
          %      day no.
          %      session no.
          %      eye no.
  ninpar = 5 ;
  ncoeff = 66 ;
  for iformat = 1:ncoeff ;
    forstr = [forstr ' %f'] ;
  end ;

  [mdata,cnt] = fscanf(fid,forstr) ;
  fclose(fid) ;

  ncol = ninpar + ncoeff ;
  nrow = cnt / ncol ;

  mdata = reshape(mdata,ncol,nrow)' ;
          %   since reshape works column wise

  vid = mdata(:,3) ;
          %  vector of indices for days
  vis = mdata(:,4) ;
          %  vector of indices for sessions
  vie = mdata(:,5) ;
          %  vector of indices for eyes




  mout = [vfilenum vip vid vis vie mpc] ;
          %  glue together final output matrix




  %  Save as ASCII file, using lines from gptable.m
  %
  fid = fopen(['\matlab\cornean\vire2' outsuf '.in'],'wt') ;
          %  'wt' is for "delete contents of this file and open 
          %               for writing" (with 't' for "text").
          %  Use 'ay' to append to what is already there, 'rt' for reading.



  forstr = '%6.1f %1.0f %1.0f %1.0f %1.0f' ;
  for iformat = 1:npc ;
    forstr = [forstr ' %14.11f'] ;
  end ;
  forstr = [forstr '\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) ;


  disp(' ') ;
  disp(' ') ;
  disp(['Finished constructing ASCII file: vire2' outsuf '.in']) ;




end ;


