disp('MATLAB Script File:  cocd5.m') ;
%
%   For the Cornea Study, this does a "curves as data" approach
%       to understanding the variability in the corneal groups.
%   Copied from cocd4.m, which it modifies by looking at
%       means and variances seperately
%
format compact ;

idel = 1 ;    %  1 to delete outliers (currently only for groups 7,9)

autop = 1 ;   %  1 to create a postscript file

ipause = 0 ;  %  0 to just make all pics, etc., with no pauses
              %  1 to pause between each plot

igroup = 9 ;  %  1 - Normal       (Central Image = "All K")
              %  2 - Ker Suspect  (Central Image = "All K")
              %  3 - Kerataconic  (Central Image = "All K")
              %  4 - Pre          (Central Image = "All P")
              %  5 - 1 Mo.        (Central Image = "All P")
              %  6 - 1 Yr.        (Central Image = "All P")
              %  7 - Pre - 1 Mo.  (Central Image not inter.)
              %  8 - Pre - 1 Yr.  (Central Image not inter.)
              %  9 - 1 Mo - 1 Yr  (Central Image not inter.)

ictype = 0 ;   %  0 - no centering
               %  1 - center about "central image" for this group
               %  2 - center about mean for this group

isurf = 1 ;   %  1 - for height map
              %  2 - for height map - sphere  (not working)
              %  3 - for curvature map


if igroup == 1 ;
  namestr = '\matlab\steve\cornea\fm\fnnor.asc' ;
  imagstr = 'Norm' ;
  icent = 119 ;    %  cofmap1.m index for Centering Feature Matrix
              %  119 - Mean All K
  centstr = '- all K' ;
  nincol = 1 ;
elseif igroup == 2 ;
  namestr = '\matlab\steve\cornea\fm\fnksus.asc' ;
  imagstr = 'Ksus' ;
  icent = 119 ;    %  cofmap1.m index for Centering Feature Matrix
              %  119 - Mean All K
  centstr = '- all K' ;
  nincol = 1 ;
elseif igroup == 3 ;
  namestr = '\matlab\steve\cornea\fm\fnker.asc' ;
  imagstr = 'Ker' ;
  icent = 119 ;    %  cofmap1.m index for Centering Feature Matrix
              %  119 - Mean All K
  centstr = '- all K' ;
  nincol = 1 ;
elseif igroup == 4 ;
  namestr = '\matlab\steve\cornea\fm\fnpre.asc' ;
  imagstr = 'Pre' ;
  icent = 122 ;    %  cofmap1.m index for Centering Feature Matrix
              %  122 - Mean all P
  centstr = '- all P' ;
  nincol = 1 ;
elseif igroup == 5 ;
  namestr = '\matlab\steve\cornea\fm\fn1mo.asc' ;
  imagstr = '1 Mo.' ;
  icent = 122 ;    %  cofmap1.m index for Centering Feature Matrix
              %  122 - Mean all P
  centstr = '- all P' ;
  nincol = 1 ;
elseif igroup == 6 ;
  namestr = '\matlab\steve\cornea\fm\fn1yr.asc' ;
  imagstr = '1 Yr.' ;
  icent = 122 ;    %  cofmap1.m index for Centering Feature Matrix
              %  122 - Mean all P
  centstr = '- all P' ;
  nincol = 1 ;
elseif igroup == 7 ;
  namestr = '\matlab\steve\cornea\fm\fnpre1mo.asc' ;
  imagstr = 'Pre - 1Mo' ;
  icent = 0 ;   %  no need to center differences
  centstr = '' ;
  nincol = 2 ;
elseif igroup == 8 ;
  namestr = '\matlab\steve\cornea\fm\fnpre1yr.asc' ;
  imagstr = 'Pre - 1Yr' ;
  icent = 0 ;   %  no need to center differences
  centstr = '' ;
  nincol = 2 ;
elseif igroup == 9 ;
  namestr = '\matlab\steve\cornea\fm\fn1mo1yr.asc' ;
  imagstr = '1Mo - 1Yr' ;
  icent = 0 ;   %  no need to center differences
  centstr = '' ;
  nincol = 2 ;
end ;


nr = 10 ;
nc = 7 ;



disp('Loading feature matrices for this group') ;
fid = fopen(namestr,'r') ;
if nincol == 1 ;
  mfindx = fscanf(fid,'%4d') 
elseif nincol == 2 ;
  mfindx = fscanf(fid,'%4d %4d') ;
  mfindx = reshape(mfindx,2,length(mfindx)/2)' ;
end ;

%  Used this for checking
%mfindx
%pause ;


prestr = '\matlab\steve\cornea\fm\fm42' ;
mfa = [] ;
for i = 1:size(mfindx,1) ;

  if mfindx(i,1) < 10 ;
    eval(['load ' prestr '000' num2str(mfindx(i,1))]) ;
  elseif mfindx(i,1) < 100 ;
    eval(['load ' prestr '00' num2str(mfindx(i,1))]) ;
  elseif mfindx(i,1) < 1000 ;
    eval(['load ' prestr '0' num2str(mfindx(i,1))]) ;
  else ;
    eval(['load ' prestr num2str(mfindx(i,1))]) ;
  end ;
          %  This creates a number of variables, especially fm,
          %  a 19 x 7 "feature matrix"

  if nincol == 2 ;    %  Then need to subtract second feature matrix
    fm1 = fm ;
    if mfindx(i,2) < 10 ;
      eval(['load ' prestr '000' num2str(mfindx(i,2))]) ;
    elseif mfindx(i,2) < 100 ;
      eval(['load ' prestr '00' num2str(mfindx(i,2))]) ;
    elseif mfindx(i,2) < 1000 ;
      eval(['load ' prestr '0' num2str(mfindx(i,2))]) ;
    else ;
      eval(['load ' prestr num2str(mfindx(i,2))]) ;
    end ;
          %  This creates a number of variables, especially fm,
          %  a 19 x 7 "feature matrix"
    fm = fm1 - fm ;
  end ;

  %  Convert c-s to 10 x 7 amplitude - phase form, 
  %        using block from cofmtst1.m
  mamp = abs(fm(1,:)) ;
  mpha = atan2(0,fm(1,:)) * 180 / pi ;
  for i = 1:9 ;
      vamp = sqrt(fm(2*i,:).^2 + fm(2*i+1,:).^2) ;
    mamp = [mamp; vamp] ;
      vpha = atan2(fm(2*i+1,:),fm(2*i,:)) * 180 / pi ;
        vind = find(vpha < 0) ;
      vpha(vind) = vpha(vind) + 360 ;
    mpha = [mpha; vpha] ;
  end ;

  %  Feature matrix form for output
  fmamp = fm(1,:) ;   %  To get signs right
  fmamp = [fmamp; mamp(2:10,:)] ;  %  Use amplitude for rest
  mfa = [mfa fmamp(:)] ;
          %  put matrix into a column vector and add to right
end ;


if igroup == 7 | igroup == 9 ;
  if idel == 1 ;
    if igroup == 7 ;
      mfa = [mfa(:,1:26) mfa(:,28:42)] ;
    elseif igroup == 9 ;
      mfa = [mfa(:,1:26) mfa(:,28)] ;
    end ;
  end ;
end ;


if  ictype > 0  &  nincol == 1  ;
  disp('Do Centering') ;
  if ictype == 1 ;   %  Then do "central image" centering
    if icent == 119 ;
    %  imagstr = 'Mean All K Study' ;
              %  All 3 groups in Keratoconus study
              %          - called kstud by NKT
              %          - called ker by JTZ
      fmistr = 'mn42ker.mat' ;
    elseif icent == 122 ;
    %  imagstr = 'Mean All P Study' ;
              %  All 4 groups in PRK study
              %          - called pstud by NKT
              %          - called pos by JTZ
              %                   nor by JTZ are not used (same as pre).
      fmistr = 'mn42pos.mat' ;
    end ;
    ['load \matlab\steve\cornea\fm\' fmistr]
    eval(['load \matlab\steve\cornea\fm\' fmistr]) ;
    mamp = mnamp ;
    mpha = mnpha ;
    %  Feature matrix form for output
    fmamp = mamp(1,:) .* cos(mpha(1,:) * pi / 180) ;
              %  Here Jinting used names mnamp & mnpha
    fmamp = [fmamp; mamp(2:10,:)] ;  %  Use amplitude for rest
    
    mfa = mfa - vec2mat(fmamp(:),size(mfindx,1)) ;
  elseif ictype == 2 ;  %  Then Center at Mean

    mfa = mfa - vec2mat(mean(mfa')',size(mfindx,1)) ;

    centstr = '- mean' ;
  end ;
else ;  %  No centering is done
  centstr = '' ;
end ;



d = size(mfa,1) ;
n = size(mfa,2) ;

disp(['joint feature matrix has size ' num2str(d) ' rows ']) ;
disp(['and ' num2str(n) ' cols']) ;


fclose('all') ;



disp('   Working on Sums of Squares') ;
vssf = zeros(10,1) ;
vssfm = zeros(10,1) ;
for ifour = 1:10 ;    %  Loop through Fourier levels

  for i = 1:n ;    %  Loop through individuals
    legcoef = mfa((0:10:60)+ifour,i) ;
    vssf(ifour) = vssf(ifour) + sum(legcoef .^ 2) ;
  end ;

  vssfm(ifour) = n * sum(mean( mfa((0:10:60)+ifour,:)' ) .^ 2) ;
end ;


sst = sum(sum(mfa.^2)) ;
sstm = n * sum(mean(mfa').^2) ;


mssfl = [] ;
mssflm = [] ;
for ifour = 1:10 ;    %  Loop through Fourier levels
  disp(['    Working on Fourier ' num2str(ifour - 1)]) ;

  vssfl = zeros(1,7) ;
  for i = 1:n ;    %  Loop through individuals
    legcoef = mfa((0:10:60)+ifour,i) ;
    vssfl = vssfl + (legcoef .^ 2)' ;
  end ;
  vssflm = n * mean(mfa((0:10:60)+ifour,:)') .^ 2 ;

  mssfl = [mssfl; vssfl] ;
  mssflm = [mssflm; vssflm] ;


end ;



%  Reset Graphics to get proper blanks
disp('   Working on Mean Histograms') ;
figure(1) ;
clf ;

subplot(3,4,1) ;
  bar(0:9,100 * vssfm / sstm) ;
    title([imagstr ', Raw ' centstr]) ;
    axis([-.5 9.5 0 100]) ;
    text(-.2,92,'% of SSM,T by F') ;
    if idel == 1 ;
    text(-.2,84,'Outlier Deleted') ;
    end ;

for ifour = 1:10 ;    %  Loop through Fourier levels
  subplot(3,4,ifour+1) ;
    bar(0:6,100 * mssflm(ifour,:) / vssfm(ifour)) ;
      title([imagstr ', F' num2str(ifour-1) ' ' centstr]) ;
      axis([-.5 6.5 0 100]) ;
      text(-.2,92,'% of SSM,F by L') ;
      if idel == 1 ;
      text(-.2,84,'Outlier Deleted') ;
      end ;
end ;

if autop == 1 ;
    orient landscape ;
  eval(['print -dps \matlab\steve\ps\cd5' num2str(igroup)...
               num2str(isurf) num2str(ictype) 'm.ps']) ;
end ;



if ipause == 1 ;
  disp('Hit any key to do Residual Histograms') ;
  pause ;
end ;
clf ;


subplot(3,4,1) ;
  bar(0:9,100 * (vssf - vssfm) ./ (sst - sstm)) ;
    title([imagstr ', Raw ' centstr]) ;
    axis([-.5 9.5 0 100]) ;
    text(-.2,92,'% of SSR,T by F') ;

for ifour = 1:10 ;    %  Loop through Fourier levels
  subplot(3,4,ifour+1) ;
    bar(0:6,100 * (mssfl(ifour,:) - mssflm(ifour,:)) ...
                             ./ (vssf(ifour) - vssfm(ifour))) ;
      title([imagstr ', F' num2str(ifour-1) ' ' centstr]) ;
      axis([-.5 6.5 0 100]) ;
      text(-.2,92,'% of SSR,F by L') ;
end ;

if autop == 1 ;
    orient landscape ;
  eval(['print -dps \matlab\steve\ps\cd5' num2str(igroup)...
               num2str(isurf) num2str(ictype) 'r.ps']) ;
end ;



if ipause == 1 ;
  disp('Hit any key to do ratio bars') ;
  pause ;
end ;
clf ;


for ifour = 1:10 ;    %  Loop through Fourier levels
  subplot(3,4,ifour+1) ;
    bar(0:6,sqrt(mssflm(ifour,:)) ./ sqrt(mssfl(ifour,:))) ;
      title([imagstr ', F' num2str(ifour-1) ' ' centstr]) ;
      axis([-.5 6.5 0 1]) ;
      text(-.2,.92,'M/T Ratio, each L') ;
      if idel == 1 ;
        text(-.2,.84,'Outlier Deleted') ;
      end ;
end ;

if autop == 1 ;
    orient landscape ;
  eval(['print -dps \matlab\steve\ps\cd5' num2str(igroup)...
               num2str(isurf) num2str(ictype) 's.ps']) ;
end ;



if ipause == 1 ;
  disp('Hit any key to do Z-score bars') ;
  pause ;
end ;
clf ;


shat = mean(mean(sqrt(mssfl - mssflm))) ;
for ifour = 1:10 ;    %  Loop through Fourier levels
  subplot(3,4,ifour+1) ;
    bar(0:6, (sqrt(mssflm(ifour,:)) ./ shat) ) ;
      title([imagstr ', F' num2str(ifour-1) ' ' centstr]) ;
      axis([-.5 6.5 0 10]) ;
      text(-.2,9.2,'Z-scores, each L') ;
      if idel == 1 ;
        text(-.2,8.4,'Outlier Deleted') ;
      end ;
end ;

if autop == 1 ;
    orient landscape ;
  eval(['print -dps \matlab\steve\ps\cd5' num2str(igroup)...
               num2str(isurf) num2str(ictype) 'z.ps']) ;
end ;



%  Used these for Checking
%disp('Check these are 0') ;
%abs(sum(sum(mssfl)) - sst) 
%sum(abs(vssf - sum(mssfl')')) 
%pause ;

vssl = sum(mssfl) ;
vsslm = sum(mssflm) ;


%  Write results to a file
  foutstr = ['\matlab\steve\ps\cd4' num2str(igroup) ...
                       num2str(ictype) '.out'] ;
fid = fopen(foutstr,'wt') ;
    titlstr1 = 'Output from the MATLAB Script cocd5.m,   ' ;
    titlstr1 = [titlstr1,date] ;
    titlstr2 = '  Studies Mean and Resid Sums of Squares for Cornea Groups' ;
  cntbytes = fprintf(fid,'%1s\n',titlstr1) ;
  cntbytes = fprintf(fid,'%1s\n\n\n',titlstr2) ;

  if idel == 1 ;
    titlstr3 = '      When Outlier was deleted:' ;
  cntbytes = fprintf(fid,'%1s\n\n',titlstr3) ;
  end ;

  hdrstr1 = ['For the Group: ' imagstr ',    Centered by: ' centstr] ;
  cntbytes = fprintf(fid,'%1s\n\n\n',hdrstr1) ;

  hdrstr2 = ['                    % of Total Sum of Squares, Mean Only'] ;
  hdrstr3 = ['                L0      L1      L2      L3      L4'] ;
  hdrstr3 = [hdrstr3, '      L5      L6'] ;
  cntbytes = fprintf(fid,'%1s\n',hdrstr2) ;
  cntbytes = fprintf(fid,'%1s\n',hdrstr3) ;

    fstr =   '       all    %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n\n' ;
    vout = 100 * vsslm / sstm ;
  cntbytes = fprintf(fid,fstr,vout) ;
  for iout = 0:9 ;
      fstr = ['F' num2str(iout)] ;
      fstr = [fstr, '  %7.3f   %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n'] ;
      vout = 100 * [vssfm(iout+1) mssflm(iout+1,:)] / sstm ;
    cntbytes = fprintf(fid,fstr,vout) ;
  end ;


  hdrstr2 = ['                    % of Total Sum of Squares, Resid Only'] ;
  hdrstr3 = ['                L0      L1      L2      L3      L4'] ;
  hdrstr3 = [hdrstr3, '      L5      L6'] ;
  cntbytes = fprintf(fid,'\n\n%1s\n',hdrstr2) ;
  cntbytes = fprintf(fid,'%1s\n',hdrstr3) ;

    fstr =   '       all    %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n\n' ;
    vout = 100 * (vssl - vsslm) / (sst - sstm) ;
  cntbytes = fprintf(fid,fstr,vout) ;
  for iout = 0:9 ;
      fstr = ['F' num2str(iout)] ;
      fstr = [fstr, '  %7.3f   %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n'] ;
      vout = 100 * [(vssf(iout+1) - vssfm(iout+1)) ...
                     (mssfl(iout+1,:) - mssflm(iout+1,:))] / (sst - sstm) ;
    cntbytes = fprintf(fid,fstr,vout) ;
  end ;


  hdrstr2 = ['                                  Z-scores'] ;
  hdrstr3 = ['                L0      L1      L2      L3      L4'] ;
  hdrstr3 = [hdrstr3, '      L5      L6'] ;
  cntbytes = fprintf(fid,'\n\n%1s\n',hdrstr2) ;
  cntbytes = fprintf(fid,'%1s\n',hdrstr3) ;

    fstr =   '       all    %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n\n' ;
    vout = 100 * vsslm / sstm ;
  cntbytes = fprintf(fid,fstr,vout) ;
  for iout = 0:9 ;
      fstr = ['F' num2str(iout)] ;
      fstr = [fstr, '  %7.3f   %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n'] ;
      vout = sqrt(n) * [
       (sqrt(vssfm(iout+1)) ./ sqrt(vssf(iout+1) - vssfm(iout+1))) ...
       (sqrt(mssflm(iout+1,:)) ./ sqrt(mssfl(iout+1,:) - mssflm(iout+1,:)))] ;
    cntbytes = fprintf(fid,fstr,vout) ;
  end ;



fclose('all') ;


