disp('Running MATLAB script file MarronPerouList1FigA.m') ;
%
%    This makes Figure A for the Marron-Perou List paper.
%    which shows the example graphic illustrating the Cluster Index
%
%    Copied from GeneArray6.m
%    which was for analyzing Chuck Perou's two Gene Array data sets,
%    this was essentially ipart == 15, in that program
%


ipart = 22  ;    %  0 - Read combined ASCII file, 
                %  15 - Examples Illustrating Cluster Index


viscale = [3] ;     %  Index of scaling type
                    %  1 - raw ratios
                    %  2 - log2(raw ratios)
                    %  3 - row median normalized log2 ratios (Chuck's choice)
                    %  4 - row median & MAD standardized
                    %  5 - row median & max standardized



vidisc = [1] ;    %  Index of dicrimination method (when needed)
                  %  1 - Cluster Index




vntrim = [2; 5; 10; 20; 30; 50; 70; 100; 200; 459] ;
          %  numbers to trim to for discrimination
nntrim = length(vntrim) ;




asciifilestr = 'GA6OnlyInputs.txt' ;
filestr = 'GA6DataAll' ;

%  set up class labels
%
classlabelcellstr = {'Basal', ...
                     'ERB', ...
                     'Normal', ...
                     'Lum A', ...
                     'Lum BC'} ;
nclust = size(classlabelcellstr,2) ;
classcolorcellstr = {'m', ...
                     'c', ...
                     'g', ...
                     'b', ...
                     'r'} ;



%   Use preloaded version of main data

eval(['load ' filestr ' ;']) ;


mdataint = mdata(:,flagintrinsic) ;
mdatafumi = mdata(:,flagfumi) ;

flagclass2int = flagclass2(flagintrinsic) ;
flagclass2fumi = flagclass2(flagfumi) ;

d = size(mdata,1) ;
nint = size(mdataint,2) ;
nfumi = size(mdatafumi,2) ;

figh = figure(1) ;
clf ;




%  Set up color map
%
genecomap = [zeros(32,1) linspace(1,0,32)' zeros(32,1)] ;
genecomap = [genecomap; [linspace(0,1,32)' zeros(32,1) zeros(32,1)]] ;
genecomap = [genecomap; [0.5 0.5 0.5]] ;
        %  add a gray color
colormap(genecomap) ;


bottomi = -0.5 ;
topi = d + 0.5 ;
left = -0.5 ;
rightint = nint + 0.5 ;




%  do scaling
%
for iscale = viscale ;

  if iscale == 1 ;    %  raw ratios

    mdataintrs = mdataint ;
    mdatafumirs = mdatafumi ;

    ccenter = 1 ;
    cradius = (log(8) + 1) / 2 ;
        %  color scale, center and radius

    scaletitstr = 'Raw Ratios' ;
    scaleoutstr = 'RawR' ;

  elseif iscale == 2 ;    %  log2(raw ratios)

    mdataintrs = log(mdataint) / log(2) ;
    mdatafumirs = log(mdatafumi) / log(2) ;
        %  use log, not log2, for better handling of NaNs

    ccenter = 0 ;
    cradius = log2(8) ;
        %  color scale, center and radius

    scaletitstr = 'Log2 Ratios' ;
    scaleoutstr = 'Log2R' ;

  elseif iscale == 3 ;    %  row median normalized log2 ratios (Chuck's choice)

    mdataintrs = log(mdataint) / log(2) ;
    mdatafumirs = log(mdatafumi) / log(2) ;

    vmedian = (nanmedian(mdataintrs'))' ;
        %  double transpose, since nanmedian works down columns only
        %  nanmedian correctly handles NaNs

    mdataintrs = mdataintrs - vec2matSM(vmedian,nint) ;
    mdatafumirs = mdatafumirs - vec2matSM(vmedian,nfumi) ;

    ccenter = 0 ;
    cradius = log(8) ;
        %  color scale, center and radius

    scaletitstr = 'Log2 Ratios - med' ;
    scaleoutstr = 'L2Rmm' ;

  elseif iscale == 4 ;    %  row median & MAD standardized

    mdataintrs = log(mdataint) / log(2) ;
    mdatafumirs = log(mdatafumi) / log(2) ;

    vmedian = (nanmedian(mdataintrs'))' ;
        %  double transpose, since nanmedian works down columns only
        %  nanmedian correctly handles NaNs

    mdataintrs = mdataintrs - vec2matSM(vmedian,nint) ;
    mdatafumirs = mdatafumirs - vec2matSM(vmedian,nfumi) ;
    
    vmad = (nanmedian(abs(mdataint')))' ;
        %  double transpose, since nanmedian works down columns only
        %  nanmedian correctly handles NaNs
    vmad = vmad / (norminv(0.75) - norminv(0.25)) ;
        %  adjust to scale of standard deviation

    mdataintrs = mdataintrs ./ vec2matSM(vmad,nint) ;
    mdatafumirs = mdatafumirs ./ vec2matSM(vmad,nfumi) ;
    
    ccenter = 0 ;
    cradius = 2 ;
        %  color scale, center and radius

    scaletitstr = 'Log2(R - med) / mad' ;
    scaleoutstr = 'L2Rmmom' ;

  elseif iscale == 5 ;    %  row median & max standardized

    mdataintrs = log(mdataint) / log(2) ;
    mdatafumirs = log(mdatafumi) / log(2) ;

    vmedian = (nanmedian(mdataintrs'))' ;
        %  double transpose, since nanmedian works down columns only
        %  nanmedian correctly handles NaNs

    mdataintrs = mdataintrs - vec2matSM(vmedian,nint) ;
    mdatafumirs = mdatafumirs - vec2matSM(vmedian,nfumi) ;
    
    vmax = (nanmax(abs(mdataintrs')))' ;
        %  double transpose, since nanmedian works down columns only
        %  nanmax correctly handles NaNs

    mdataintrs = mdataintrs ./ vec2matSM(vmax,nint) ;
    mdatafumirs = mdatafumirs ./ vec2matSM(vmax,nfumi) ;

    ccenter = 0 ;
    cradius = 2 / 3 ;
        %  color scale, center and radius

    scaletitstr = 'Log2(R - med) / max' ;
    scaleoutstr = 'L2Romax' ;

  end ;



  %  First appropriately reorder data
  %
  %  sort columns according to classes, row 4
  %
  flagclass4int = flagclass4(flagintrinsic) ;
  [sortflag4int,vind] = sort(flagclass4int) ;
  
  vdif = sortflag4int(2:nint) - sortflag4int(1:(nint-1)) ;
  vbdry = 1.5:1:nint ;
  vbdry = vbdry(vdif > 0.5) ;


  mdataintrsco = mdataintrs(:,vind) ;
        %  use "case ordering", depending on labels

  sort4flag2int = flagclass2int(vind) ;
      %  corresponding version of row 2 labels


  %  order rows using the modified Chuck's ordering
  %  (with 3 new genes at end)
  %
  sortkey = load('GA6RowSortKey.txt') ;
  [temp,vind] = sort(sortkey) ;

  mdataintro = mdataintrsco(vind,:) ;
  mdatafumiro = mdatafumirs(vind,:) ;
      %  re-ordered versions  of rescaled data




  vig = [329, 453] ;
      %  vector of gene indices, 
      %  chosen with help of GA6PareListClustInd.txt
      %  first has the best cluster index
      %  second has the worst cluster index
  vci = [0.83879, 0.00143] ;
      %  vector of corresponding cluster indices
  commentcellstr = {'Best','Worst'} ;

  nonmissgenecomap = genecomap(1:(size(genecomap,1)-1),:) ;
  colormap(nonmissgenecomap) ;


  figure(1) ;
  clf ;
  for i = 1:2 ;    %  loop through examples

    vdata = mdataintro(vig(i),:) ;
    
    left = nanmin(vdata) ;
    right = nanmax(vdata) ;
    range = max(abs(left),abs(right)) ;
    left = -1.02 * range ;
    right = 1.02 * range ;
    
    xgrid = linspace(left,right,200)' ;


    subplot(2,1,i) ;
    
    hold on ;
    %  loop through classes
    %
    bottom = 0 ;
    vtop = [] ;
    for iclass = 1:nclust ;
    
      cdata = vdata(sortflag4int == iclass) ;
      missflag = isnan(cdata) ;
      cdata = cdata(~missflag) ;
      ncdata = length(cdata) ;

      mu = mean(cdata) ;
      sig = std(cdata) ;
      
      pdfplot = normpdf(xgrid,mu,sig) ;
      
      plot(xgrid,pdfplot,[classcolorcellstr{iclass} '-'], ...
                         'linewidth',3) ;
      
      vtop = [vtop; 1.05 * max(pdfplot)] ;

    end ;    %  of iclass loop through classes
    
    
    top = max(vtop) ;


    %  loop through classes again to add other parts
    %
    vtop = [] ;
    for iclass = 1:nclust ;

      cdata = vdata(sortflag4int == iclass) ;
      missflag = isnan(cdata) ;
      cdata = cdata(~missflag) ;
      ncdata = length(cdata) ;

        bottomi = bottom + ((iclass - 1) / nclust) * 0.6 * (top - bottom) ;
        topi = bottom + (iclass / nclust) * 0.6 * (top - bottom) ;
      plot([cdata; cdata],[bottomi * ones(1,ncdata) ; ...
                           topi * ones(1,ncdata)], ...
                 [classcolorcellstr{iclass} '-'], ...
                 'linewidth',1.5) ;

        tx = left + ((iclass - 0.7) / nclust) * (right - left) ;
        ty = bottom + 0.8 * (top - bottom) ;
      text(tx,ty,classlabelcellstr{iclass}, ...
              'Color',classcolorcellstr{iclass}, ...
              'FontSize',15) ;

    end ;    %  of iclass loop through classes

    hold off ;

    title([commentcellstr{i} ...
           ' Gene ' num2str(vig(i)) ...
           ', Cluster Index = ' num2str(vci(i))], ...
           'FontSize',15) ;
    ylabel('Gaussian fit','FontSize',12) ;
    axis([left,right,0,top]) ;

    if i == 1 ;
      cbh = colorbar('horiz') ;
      set(cbh,'XTickLabel',[]) ;
      set(cbh,'XTick',[]) ;
      vpos = get(cbh,'Position') ;
          %  position vector of color bar, in pixel units?
      cbrad = vpos(3) / 2 ;
          %  color bar radius
      cbcent = vpos(1) + cbrad ;
          %  color bar center
      cbrad = cbrad * 3 / (1.02 * range) ;
      vpos = [(cbcent - cbrad) vpos(2) (2 * cbrad) vpos(4)] ;
      set(cbh,'Position',vpos) ;
    elseif i == 2 ;
      xlabel('normalized log2 color ratio','FontSize',12) ;
    end ;


  end ;    %  of i-loop through examples



    outstr = 'MarronPerouList1FigA' ;
    orient tall ;
  print('-dpsc2',outstr) ;
  print('-depsc2',outstr) ;
  print('-djpeg90',outstr) ;
  print('-dtiff',outstr) ;



end ;    %  of iscale loop



