disp('MATLAB Script File:  cosd14.m') ;
%
%   For the Cornea Study, this does a "surfaces as data" approach
%       2nd Generation, using complete data sets
%   This does standard principle components 
%       but only for the 0th Fourier coefficients 
%       for the time differences in the PRK study
%       (i.e. in spherical direction only)
%     Color maps are now replaced by curve overlays
%
format compact ;

autop = 1 ;   %  1 to create a postscript file

idif = 612 ;    %  01    pre - 1 mo.
               %  06    pre - 6 mo.
               %  012   pre - 12 mo.
               %  16    1 mo. - 6 mo.
               %  112    1 mo. - 12 mo.
               %  612    6 mo. - 12 mo.

izero = 1 ;   %  0 - Use full feature matrix
              %  1 - Use only Fourier 0 terms



r0 = 4.2 ;    %  Usual radius of analysis

nrgrid = 43 ;
prgrid = linspace(0,4.2,nrgrid)' ;
          %  grid for plotting profiles          
ptgrid = zeros(nrgrid,1) ;



disp('Working on PRK Population ') ;
%  First load in file info from file prkdkey.dat, 
%  which was modified from fcpairs.txt
%  (everything but the data was edited out).
%    1st col is difference code (same as idif above)
%    2nd col is Patient number
%    3nd col is 1st file id (add last 4 digits to "fm42" to get file suffix)
%    4th col is 2nd file id (add last 4 digits to "fm42" to get file suffix)
fid = fopen('\matlab\steve\cornea\fm\prkdkey.dat','r') ;
  frmtstr = '%g %g %g %g' ;
    mfilein = fscanf(fid,frmtstr,[4,inf]) ;
    mfilein = mfilein' ;
          %  since data rows are read in as columns
  fclose(fid) ;

vdcode = mfilein(:,1) ;
          %  vector of difference codes
vpatid = mfilein(:,2) ;
          %  vector of patient id's
vfileid1 = mfilein(:,3) ;
          %  vector of 1st file id's
vfileid2 = mfilein(:,4) ;
          %  vector of 2nd file id's
mfilein = 0 ;
          %  won't be using this again


if idif == 01 ;
  imagstr = 'PRK Pre-1' ;
  vaxis1 = [0 4.2 -40 40] ;
  vaxis2 = [0 4.2 -20 20] ;
elseif idif == 06 ;
  imagstr = 'PRK Pre-6' ;
  vaxis1 = [0 4.2 -40 40] ;
  vaxis2 = [0 4.2 -20 20] ;
elseif idif == 012 ;
  imagstr = 'PRK Pre-12' ;
  vaxis1 = [0 4.2 -40 40] ;
  vaxis2 = [0 4.2 -20 20] ;
elseif idif == 16 ;
  imagstr = 'PRK 1-6' ;
  vaxis1 = [0 4.2 -40 40] ;
  vaxis2 = [0 4.2 -15 15] ;
elseif idif == 112 ;
  imagstr = 'PRK 1-12' ;
  vaxis1 = [0 4.2 -40 40] ;
  vaxis2 = [0 4.2 -15 15] ;
elseif idif == 612 ;
  imagstr = 'PRK 6-12' ;
  vaxis1 = [0 4.2 -40 40] ;
  vaxis2 = [0 4.2 -15 15] ;
end ;
vflag = (vdcode == idif) ;



disp(['Loading feature matrices for ' imagstr]) ;

vpatid = vpatid(vflag) ;
vfileid1 = vfileid1(vflag) ;
vfileid2 = vfileid2(vflag) ;
nimag = length(vpatid) ;


%  Loop through feature matrix files, and combine into big data matrix
mfv = [] ;
for iimag = 1:nimag ;

  %  Get first part of difference
  fid1 = vfileid1(iimag) ;
  if fid1 < 100 ;
    fid1s = ['00' num2str(fid1)] ;
  elseif fid1 < 1000 ;
    fid1s = ['0' num2str(fid1)] ;
  else ;
    fid1s = num2str(fid1) ;
  end ;
  eval(['load \matlab\steve\cornea\fm\fm42' fid1s '.mat']) ;
          %  This load Jinting's previously saved stuff
          %  which is a number of variables, especially fm,
          %  a 19 x 7 "feature matrix"
  dfm = fm ;

  %  Get 2nd part of difference
  fid2 = vfileid2(iimag) ;
  if fid2 < 100 ;
    fid2s = ['00' num2str(fid2)] ;
  elseif fid2 < 1000 ;
    fid2s = ['0' num2str(fid2)] ;
  else ;
    fid2s = num2str(fid2) ;
  end ;
  eval(['load \matlab\steve\cornea\fm\fm42' fid2s '.mat']) ;
          %  This load Jinting's previously saved stuff
          %  which is a number of variables, especially fm,
          %  a 19 x 7 "feature matrix"
  dfm = dfm - fm ;


  mfv = [mfv dfm(:)] ;
          %  put matrix into a column vector and add to right
          %  Including the feature matrix for this image
end ;





%  Now zero out terms in the feature matrix,
if izero == 0 ;     %  Use Full Feature Matrix
  zstr = ' Full FM, ' ;
elseif izero == 1 ;     %  Use only Fourier 0 terms
  zstr = ' Only F0, ' ;
  mask = [ones(1,7); zeros(18,7)] ; 
          %  1's where want to keep coefficients in 19x7 form
  mask = vec2mat(mask(:),size(mfv,2)) ;
          %  First put into a long column, then extend horizontally
  mfv = mfv .* mask ;
end ;



disp(['feature matrix has size ' num2str(size(mfv,1)) ' rows ']) ;
disp(['and ' num2str(size(mfv,2)) ' cols']) ;




disp('Working on eigenector decomposition') ;
d = size(mfv,1) ;
n = size(mfv,2) ;
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 = median(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 = median(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 = median(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 = median(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 ;


fclose('all') ;





disp('Inverting feature vectors to get profile curves') ;
moutr = [] ;
for ipat = 1:size(mfv,2) ;    %  loop through patients and reconstruct
                              %        raw profiles
  fmraw = reshape(mfv(:,ipat),19,7) ;
  [mhr,mhpr,mhppr] = corhpp(fmraw',prgrid,ptgrid,r0) ;
  voutr = -mhppr ./ (1 + mhpr.^2).^(3/2) ;
  moutr = [moutr, voutr * 337.5] ;
  
end ;

fmmean = reshape(vmean,19,7) ;
[mhr,mhpr,mhppr] = corhpp(fmmean',prgrid,ptgrid,r0) ;
moutm = -mhppr ./ (1 + mhpr.^2).^(3/2) ;
moutm = moutm * 337.5 ;

fmpc1 = reshape(ev1,19,7) ;
[mhr,mhpr,mhppr] = corhpp(fmpc1',prgrid,ptgrid,r0) ;
mout1 = -mhppr ./ (1 + mhpr.^2).^(3/2) ;
mout1 = mout1 * 337.5 ;

fmpc2 = reshape(ev2,19,7) ;
[mhr,mhpr,mhppr] = corhpp(fmpc2',prgrid,ptgrid,r0) ;
mout2 = -mhppr ./ (1 + mhpr.^2).^(3/2) ;
mout2 = mout2 * 337.5 ;

fmpc3 = reshape(ev3,19,7) ;
[mhr,mhpr,mhppr] = corhpp(fmpc3',prgrid,ptgrid,r0) ;
mout3 = -mhppr ./ (1 + mhpr.^2).^(3/2) ;
mout3 = mout3 * 337.5 ;

fmpc4 = reshape(ev4,19,7) ;
[mhr,mhpr,mhppr] = corhpp(fmpc4',prgrid,ptgrid,r0) ;
mout4 = -mhppr ./ (1 + mhpr.^2).^(3/2) ;
mout4 = mout4 * 337.5 ;




disp('Making Raw Plot') ;
subplot(3,4,1) ;
  plot(prgrid,moutr) ;
    title([imagstr zstr ' Raw']) ;
    axis(vaxis1) ;


disp('Making Mean Plot') ;
subplot(3,4,5) ;
  plot(prgrid,moutm) ;
    title([imagstr zstr ' Mean']) ;
      ssstr = ['% of SS,T = ' num2str(100 * ssm / sst)] ;
    text(1,45,ssstr) ;  
    axis(vaxis1) ;

disp('Making Mean Resid. Plot') ;
resid = moutr - vec2mat(moutm,size(moutr,2)) ;
subplot(3,4,9) ;
  plot(prgrid,resid) ;
    title([imagstr zstr ' Mean R']) ;
      ssstr = ['% of SS,T = ' num2str(100 * (sst - ssm) / sst)] ;
    text(1,6,ssstr) ;  
    axis(vaxis2) ;


disp('Making PC1 Plot') ;
mproj1 = mout1 * pc1' ;
subplot(3,4,2) ;
  plot(prgrid,mproj1) ;
    title([imagstr zstr ' PC1']) ;
      ssstr = ['% of SS,MR = ' num2str(100 * ss1 / ssmr)] ;
    text(1,6,ssstr) ;  
    axis(vaxis2) ;

disp('Making PC1 Resid. Plot') ;
resid = resid - mproj1 ;
subplot(3,4,3) ;
  plot(prgrid,resid) ;
    title([imagstr zstr ' Mean R']) ;
      ssstr = ['% of SS,MR = ' num2str(100 * (ssmr - ss1) / ssmr)] ;
    text(1,6,ssstr) ;  
    axis(vaxis2) ;


disp('Making PC2 Plot') ;
mproj2 = mout2 * pc2' ;
subplot(3,4,6) ;
  plot(prgrid,mproj2) ;
    title([imagstr zstr ' PC2']) ;
      ssstr = ['% of SS,1R = ' num2str(100 * ss2 / ss1r)] ;
    text(1,6,ssstr) ;  
    axis(vaxis2) ;

disp('Making PC2 Resid. Plot') ;
resid = resid - mproj2 ;
subplot(3,4,7) ;
  plot(prgrid,resid) ;
    title([imagstr zstr ' Mean R']) ;
      ssstr = ['% of SS,1R = ' num2str(100 * (ss1r - ss2) / ss1r)] ;
    text(1,6,ssstr) ;  
    axis(vaxis2) ;


disp('Making PC3 Plot') ;
mproj3 = mout3 * pc3' ;
subplot(3,4,10) ;
  plot(prgrid,mproj3) ;
    title([imagstr zstr ' PC3']) ;
      ssstr = ['% of SS,2R = ' num2str(100 * ss3 / ss2r)] ;
    text(1,6,ssstr) ;  
    axis(vaxis2) ;

disp('Making PC3 Resid. Plot') ;
resid = resid - mproj3 ;
subplot(3,4,11) ;
  plot(prgrid,resid) ;
    title([imagstr zstr ' Mean R']) ;
      ssstr = ['% of SS,2R = ' num2str(100 * (ss2r - ss3) / ss2r)] ;
    text(1,6,ssstr) ;  
    axis(vaxis2) ;


disp('Making PC4 Plot') ;
mproj4 = mout4 * pc4' ;
subplot(3,4,4) ;
  plot(prgrid,mproj4) ;
    title([imagstr zstr ' PC4']) ;
      ssstr = ['% of SS,3R = ' num2str(100 * ss4 / ss3r)] ;
    text(1,6,ssstr) ;  
    axis(vaxis2) ;

disp('Making PC4 Resid. Plot') ;
resid = resid - mproj4 ;
subplot(3,4,8) ;
  plot(prgrid,resid) ;
    title([imagstr zstr ' Mean R']) ;
      ssstr = ['% of SS,3R = ' num2str(100 * (ss3r - ss4) / ss3r)] ;
    text(1,6,ssstr) ;  
    axis(vaxis2) ;


disp('Making log power Plot') ;
resid = resid - mproj4 ;
subplot(3,4,12) ;
  bar(0:4,log10([ssm ss1 ss2 ss3 ss4])+5) ;
    title('log10 SS of PC''s + 5') ;
    axis([-.5, 4.5, 0, 7]) ;



if autop == 1 ;
  orient landscape ;
  if isurf == 2 ;
    eval(['print \matlab\steve\ps\co14' num2str(idif) 'h.ps -dps']) ;
  elseif isurf == 3 ;
    eval(['print \matlab\steve\ps\co14' num2str(idif) 'c.ps -dps']) ;
  end ;
end ;



