function gpcd2(data,itype,iout,psfname) 
% GPCD2, General Purpose Curves as Data analysis, 2nd generation
%     Makes a single page with eigenvectors and projections
%   Can use first 1, 2, 3 or 4 arguments.
% Inputs:
%     data    - d x n matrix of data, each column vector is 
%                      a "d-dim digitized curve"
%     itype   - scalar, with analysis type as:
%                 1 - (default) Standard Principle Components
%                 2 - Correlation Principle Components
%                 3 - Spearman's Correlation PCA
%                 4 - Shrink to sphere PCA
%                 5 - map to ellipse (determined by MAD) PCA
%     iout    - plots to make:
%                 1 - just simple page of eignevectors and projections
%                 2 - add a SiZer and Q-Q analysis
%                 3 - add a draftsman's plot
%                 4 - both SiZer and draftman's plots
%     psfname - string with prefix for name of postscript file
%                    (will automatically add ".ps")
%                    (will be a color postscript file for SiZer-Sicon maps)
% Outputs:
%     Only graphics, in Figures 1 - 3, unless postscript files are created
%
% Assumes path can find personal functions:
%    bwsjpib.m
%    gpkde.m

%    Copyright (c) J. S. Marron 1998



%  Set parameters and defaults according to number of input arguments
%
if nargin == 1 ;    %  only 1 argument input, use default iitype
  iitype = 1 ;
else ;              %  itype was specified, use that
  iitype = itype ;
end ;

if nargin <= 2 ;    %  at most 2 arguments input, use default iiout
  iiout = 1 ;
else ;              %  iout was specified, use that
  iiout = iout ;
end ;

if nargin <= 3 ;    %  at most 3 arguments input, create no postscript file
  ipsfname = [] ;    %  Default
else ;
  ipsfname = psfname ;    %  Create Postscript file, with given prefix
end ;


%  set preliminary stuff
%
figure(1) ;
clf ;
figcount = 1 ;
rand('seed',20983743) ;



%  Get statistics and do eigenvectors
%
d = size(data,1) ;
         %  dimension of each data curve
n = size(data,2) ;
         %  number of data curves
nev = min([d,n,4]) ;
         %  number of eigenvectors to display


%  get mean and subtract, and get covariance matrix
%
if iitype == 1 ;    %  Standard PCA
  vmean = mean(data')' ;
  ssm = n * sum(vmean.^2) ;
  mresid = data - vec2mat(vmean,n) ;
  ssmr = sum(sum(mresid .^ 2)) ;

  mcov = cov(mresid') ;
          %  Get covariance matrix, transpose, since want "coordinates
          %                      as variables"
  typestr = 'Standard PCA' ;

elseif iitype == 2 ;    %  Correlation PCA
  vmean = mean(data')' ;
  ssm = n * sum(vmean.^2) ;
  mresid = data - vec2mat(vmean,n) ;
  ssmr = sum(sum(mresid .^ 2)) ;

  mcov = corrcoef(mresid') ;
          %  Get correlation matrix, transpose, since want "coordinates
          %                      as variables"
  typestr = 'Correlation PCA' ;

elseif iitype == 3 ;    %  Spearman Correlation PCA
  vmean = median(data')' ;
  ssm = n * sum(vmean.^2) ;
  mresid = data - vec2mat(vmean,n) ;
  ssmr = sum(sum(mresid .^ 2)) ;

  [temp,msind] = sort(mresid') ;
          %  Transpose, since want "coordinates as variables"
          %  Sort gives indices for doing sort
  mrank = [] ;
  vind = (1:n)' ;
  for i = 1:d ;
    mrank = [mrank, vind(msind(:,i))] ;
          %  matrix of ranks
  end ;
  mcov = corrcoef(mresid') ;
          %  Get covariance matrix, transpose, since want "coordinates
          %                      as variables"
  typestr = 'Spearman Corr PCA' ;

elseif iitype == 4 ;    %  Shrink to sphere PCA
  vmean = rmean(data',(10^-6),20,0)' ;
          %  accuracy parameters, and 0 for no screen writes
  ssm = n * sum(vmean.^2) ;
  mresid = data - vec2mat(vmean,n) ;
  ssmr = sum(sum(mresid .^ 2)) ;

  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
  typestr = 'Spherical PCA' ;

elseif iitype == 5 ;    %  Shrink to Ellipse PCA

  vmad = gpmad(data')' ;
          %  column vector of MADs of each row

  flag0mad = (vmad == 0) ;
  n0mad = sum(flag0mad) ;
          %  one where vmad = 0 
  if n0mad > 0 ;    %  then have some zero mad's
    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
    disp('!!!   Warning from gpcd2.m:     !!!') ;
    disp(['!!!   Have 0 mad in ' num2str(n0mad) ' coordinates']) ;
    disp('!!!   will replace with zeros   !!!') ;
    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
    vmad(flag0mad) = ones(n0mad,1) ;
  end ;

  sdata = data .* vec2mat(1 ./ vmad, n) ;
          %  make each coordinate have MAD 1
  if n0mad > 0 ;
    sdata(flag0mad) = zeros(n0mad,1) ;
  end ;

  svmean = rmean(sdata',(10^-6),20,0)' ;
          %  accuracy parameters, and 0 for no screen writes
  vmean = svmean .* vmad ;
  ssm = n * sum(vmean.^2) ;
  smresid = sdata - vec2mat(svmean,n) ;
  mresid = smresid .* vec2mat(vmad, n) ;
  ssmr = sum(sum(mresid .^ 2)) ;

  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
  typestr = 'Elliptical PCA' ;

end ;



[meigvec,meigval] = eig(mcov) ;
veigval = (diag(meigval)) ;
          %  Get eigenvalues and eigenvectors
[veigval,vind] = sort(veigval) ;
meigvec = meigvec(:,vind) ;
          %  Sort to avoid possible ordering problems
vev = flipud(veigval) ;
          %  put largest one first

vev = vev / sum(vev) ;
          %  put on scale of R^2
cev = cumsum(vev) ;
          %  get cumulative sum


ev1 = meigvec(:,d) ;
          %  eigenvector with largest eigenvalue
pc1 = mresid' * ev1 ;
          %  vector of projections of data onto ev1
mpc1 = ev1 * pc1' ;
          %  matrix for display
ss1 = sum(sum(mpc1.^2)) ;
ss1r = ssmr - ss1 ;
mresid1 = (mresid - mpc1) ;
vr2 = ss1 / ssmr ;
vflast = ss1r / ssmr ;

ev2 = meigvec(:,d-1) ;
          %  eigenvector with 2nd largest eigenvalue
pc2 = mresid1' * ev2 ;
          %  vector of projections of data onto ev2
mpc2 = ev2 * pc2' ;
          %  matrix for display
ss2 = sum(sum(mpc2.^2)) ;
ss2r = ss1r - ss2 ;
mresid2 = (mresid1 - mpc2) ;
vr2 = [vr2, ss2 / ssmr] ;
vflast = [vflast, ss2r / ss1r] ;


if nev > 2 ;

  ev3 = meigvec(:,d-2) ;
          %  eigenvector with 3rd largest eigenvalue
  pc3 = mresid2' * ev3 ;
          %  vector of projections of data onto ev3
  mpc3 = ev3 * pc3' ;
          %  matrix for display
  ss3 = sum(sum(mpc3.^2)) ;
  ss3r = ss2r - ss3 ;
  mresid3 = (mresid2 - mpc3) ;
  vr2 = [vr2, ss3 / ssmr] ;
  vflast = [vflast, ss3r / ss2r] ;

  if nev > 3 ;

    ev4 = meigvec(:,d-3) ;
          %  eigenvector with 3rd largest eigenvalue
    pc4 = mresid3' * ev4 ;
          %  vector of projections of data onto ev4
    mpc4 = ev4 * pc4' ;
          %  matrix for display
    ss4 = sum(sum(mpc4.^2)) ;
    ss4r = ss3r - ss4 ;
    mresid4 = (mresid3 - mpc4) ;
    vr2 = [vr2, ss4 / ssmr] ;
    vflast = [vflast, ss4r / ss3r] ;

  end ;

end ;

cr2 = cumsum(vr2) ;


%  make output plot
%
npx = 4 ;
          %  always do 4 column output
npy = nev + 1 ;
          %  number of rows output
xgrid = (1:d)' ;
min1 = min(min(data)) ;
max1 = max(max(data)) ;
range1 = max1 - min1 ;
min1 = min1 - 0.05 * range1 ;
max1 = max1 + 0.05 * range1 ;
min2 = min(min(mresid)) ;
max2 = max(max(mresid)) ;
range2 = max2 - min2 ;
min2 = min([min2; -max2]) ;
max2 = max([max2; -min2]) ;
min2 = min2 - 0.05 * range2 ;
max2 = max2 + 0.05 * range2 ;


subplot(npy,npx,1) ;
  plot(xgrid,data) ;
    title('Raw Data') ;
    axis([1,d,min1,max1]) ;


subplot(npy,npx,2) ;
  plot(xgrid,vmean) ;
    title('Mean') ;
    axis([1,d,min1,max1]) ;
      ty = min1 + .9 * (max1 - min1) ;
    text(1.5,ty,[num2str(100 * ssm / (ssm + ssmr)) '% of Tot']) ;


subplot(npy,npx,3) ;
  plot(xgrid,mresid) ;
    title('Mean Resid.') ;
    axis([1,d,min2,max2]) ;


subplot(npy,npx,4) ;
  plot(1:d,vev,'m-',1:d,cev,'c--',1:nev,vr2,'mo',1:nev,cr2,'c+') ;
    title('PC R-squareds') ;
    xlabel('PC #') ;
    ylabel('R^2') ;
    minev = min([0; vev]) ;
    if minev < -10^(-3) ;
      negevflag = (vev < 0) ;
      negevind = 1:d ;
      negevind = negevind(negevflag) ;
      negev = vev(negevflag) ;
      hold on ;
        plot(negevind,negev,'yo') ;
      hold off ;
      minev = minev - 0.05 ;
    end ;
    axis([0,d+1,minev,1]) ;
    text(1.5,.5,typestr) ;


%  loop through principal components
%
for iev = 1:nev ;

  istr = num2str(iev) ;

  subplot(npy,npx,4 * iev + 1) ;
  eval(['plot(xgrid,mpc' istr ')']) ;
    title(['PC' num2str(iev)]) ;
    axis([1,d,min2,max2]) ;
    if npx == 2 ;    %  then put this here, else put on mean +- plot
        ty = min2 + .9 * range2 ;
      text(1.5,ty,[num2str(100 * vr2(iev)) '% of MR']) ;
    end ;


  %  do residual plot
  subplot(npy,npx,4 * iev + 3) ;
  eval(['plot(xgrid,mresid' istr ')']) ;
    title(['resid, PC' num2str(iev)]) ;
    axis([1,d,min2,max2]) ;
      ty = min2 + .9 * range2 ;
    text(1.5,ty,[num2str(100 * vflast(iev)) '% of last']) ;


  %  do mean +- plot
  subplot(npy,npx,4 * iev + 2) ;
    plot(xgrid,vmean) ;
      title(['Mean +- PC' istr]) ;
      axis([1,d,min1,max1]) ;
      hold on ;
        eval(['plot(xgrid,vmean + max(pc' istr ') * ev' istr ',''--'') ;']) ;
        eval(['plot(xgrid,vmean + min(pc' istr ') * ev' istr ','':'') ;']) ;
      hold off ;
        ty = min1 + .9 * (max1 - min1) ;
      text(1.5,ty,[num2str(100 * vr2(iev)) '% of MR']) ;


  %  do kde of projections
  eval(['min3 = min(pc' istr ') ;']) ;
  eval(['max3 = max(pc' istr ') ;']) ;
  range3 = max3 - min3 ;
  min3 = min3 - 0.05 * range3 ;
  max3 = max3 + 0.05 * range3 ;
  subplot(npy,npx,4 * iev + 4) ;


    eval(['[kde,kdexgrid] = gpkde(pc' istr ',0,[min3; max3]) ;']) ;
    plot(kdexgrid,kde) ;
      title(['kde, PC' istr]) ;
      top3 = 1.05 * max(kde) ;
      axis([min3,max3,0,top3]) ;
      hold on ;
        rhts = (.5 + .1 * rand(n,1)) * top3 ;
          %  random heights
        eval(['plot(pc' istr ',rhts,''r.'') ;']) ;
          % do jitter plot
      hold off ;


end ;


if length(ipsfname) > 0 ;     %  then create postscript file with given name
    if npy > npx ;
      orient tall ;
    else ;
      orient landscape ;
    end ;
  eval(['print -dps \matlab\steve\ps\' ipsfname '.ps']) ;
end ;



if iiout == 2 | iiout == 4 ;   %  then add a Sizer - SiCon analysis

  figure(2) ;
  figcount = figcount + 1 ;

  for i = 1:nev ;
    istr = num2str(i) ;

    eval(['pcdat = pc' istr ' ;']) ;
    min3 = min(pcdat) ;
    max3 = max(pcdat) ;
    range3 = max3 - min3 ;
    min3 = min3 - 0.05 * range3 ;
    max3 = max3 + 0.05 * range3 ;

    disp(['    Working on family, pc' istr]) ;
    subplot(4,nev,i) ;

    inhp = 11 ;
    ngrid = 401 ;
    centd = mean([min3;max3]) ;

    %  Set h grid stuff, as in SiZer1
    range = max3 - min3 ;
    binw = range / (ngrid - 1) ;
    hmin = 2 * binw ;
    hmax = range ;
    vh = logspace(log10(hmin),log10(hmax),inhp) ;

    hcent = bwsjpib(pcdat) ;
    if hcent < min(vh) ;    %  if data based h below range
      disp('!!!   Warning from gpanal: databased h below this range   !!!') ;
      hcflag = 0 ;
    elseif hcent > max(vh) ;    %  if databased h above this range
      disp('!!!   Warning from gpanal: databased h above this range   !!!') ;
      hcflag = 0 ;
    else ;     %  if data based in range, then get its nearest index
      [temp, hcind] = min(abs(log10(hcent) - log10(vh))) ;
      hcind = inhp + 1 - hcind ;
                        %  since handles appear in reverse order in plot
      hcflag = 1 ; 
    end ;
    [msmoo, xgrid] = gpkde(pcdat,vh,[min3; max3; ngrid]) ;
    bottom = 0 ;
    top = 1.05 * max(max(msmoo)) ;

          %  Next lines from gpfam.m
    plot(xgrid,msmoo,'c') ;
          %  Plot most curves in cyan
      vcurvh = get(gca,'Children') ;
          %  Vector of handles for curves
      if hcflag == 1 ;
        set(vcurvh(hcind),'LineWidth',2) ;
          %  use fatter line width for data based choice
        set(vcurvh(hcind),'Color',[1 0 0]) ;
          %  use red color for one in middle
      end ;

      vax = [min3,max3,bottom,top] ;
      axis(vax) ;
      title(['Family Plot, pc' istr]) ;

      hold on ;
        yrand = vax(3) + (0.8 + 0.1 * rand(size(pcdat))) ...
                                               * (vax(4) - vax(3)) ;
               %  y coords for jittering
        plot(pcdat,yrand,'r.') ;
      hold off ;


    disp(['    Working on SiZer, pc' istr]) ;
    subplot(4,nev,4+i) ;
    gpsz1(pcdat,[min3,max3,ngrid]) ;
      title(['SiZer, pc' istr]) ;
          %  next add marking lines
      hold on ;
        if hcflag == 1 ;
          plot([min3; max3], ones(2,1)*log10(hcent),'-k') ;
        end ;
        plot(centd + 2*vh, log10(vh),':k') ;
        plot(centd - 2*vh, log10(vh),':k') ;
      hold off ;

      %  Set up colorful color map
      cocomap = [0,    0,   1; ...
                .35, .35, .35; ...
                .5,    0,  .5; ...
                 1,    0,   0; ...
                 1,   .5,   0; ...
               .35,  .35, .35; ...
                 0,    1,   0; ...
                 0,    1,   1] ; 
      colormap(cocomap) ;


    disp(['    Working on SiCon, pc' istr]) ;
    subplot(4,nev,8+i) ;
      mapout = gpsc1(pcdat,[min3,max3,ngrid]) ;
          %  get matrix of map values, so can handle both color sets at once
      mapout = mapout + 4 ;
                        %  Adjust colors in sicom, to fit
                        %  joint color range
      image([min3,max3],[log10(hmin),log10(hmax)],mapout) ;
        set(gca,'YDir','normal') ;
          colormap(cocomap) ;  
        ylabel('log10(h)') ;
        title(['SiCon, pc' istr]) ;
          %  next add marking lines
        hold on ;
          if hcflag == 1 ;
            plot([min3; max3], ones(2,1)*log10(hcent),'-k') ;
          end ;
          plot(centd + 2*vh, log10(vh),':k') ;
          plot(centd - 2*vh, log10(vh),':k') ;
        hold off ;


    disp(['    Working on Q-Q plot, pc' istr]) ;
    subplot(4,nev,12+i) ;
    xbar = mean(pcdat) ;
    s = std(pcdat) ;
    gq = xbar + s * phiinv((.5:1:n)' / n ) ;
          %  corresponding gaussian quantiles
    plot(sort(pcdat),gq,'r-',[min3;max3],[min3;max3],'c--') ;
      title(['Gaussian Q-Q, pc' istr]) ;
      eval(['xlabel(''pc' istr ''') ;']) ;
      ylabel('Normal Q') ;
      axis([min3,max3,min3,max3]) ;

  end ;


  if length(ipsfname) > 0 ;   %  then create postscript file with given name
      orient landscape ;
    eval(['print -dpsc \matlab\steve\ps\' ipsfname 's.ps']) ;
  end ;


end ;



if iiout == 3 | iiout == 4 ;    %  then add a draftsman's plot

  figure(figcount + 1) ;


  %  First work on diagonal, and put down same kde stuff
  %
  for i = 1:nev ;
    istr = num2str(i) ;

    eval(['min3 = min(pc' istr ') ;']) ;
    eval(['max3 = max(pc' istr ') ;']) ;
    range3 = max3 - min3 ;
    min3 = min3 - 0.05 * range3 ;
    max3 = max3 + 0.05 * range3 ;
    subplot(4,nev,5*(i-1)+1) ;
      eval(['[kde,kdexgrid] = gpkde(pc' istr ',0,[min3; max3]) ;']) ;
      plot(kdexgrid,kde) ;
        title(['kde, PC' istr]) ;
        top3 = 1.05 * max(kde) ;
        axis([min3,max3,0,top3]) ;
        hold on ;
          rhts = (.5 + .1 * rand(n,1)) * top3 ;
            %  random heights
          eval(['plot(pc' istr ',rhts,''r.'') ;']) ;
            % do jitter plot
        hold off ;

  end ;


  
  %  Now do off diagonals
  %
  for i = 1:nev ;

    for j = 1:nev ;


      if i ~= j ;   %  then not on diagonal, so make a plot

        istr = num2str(i) ;
        jstr = num2str(j) ;

        eval(['min3i = min(pc' istr ') ;']) ;
        eval(['max3i = max(pc' istr ') ;']) ;
        range3i = max3i - min3i ;
        min3i = min3i - 0.05 * range3i ;
        max3i = max3i + 0.05 * range3i ;

        eval(['min3j = min(pc' jstr ') ;']) ;
        eval(['max3j = max(pc' jstr ') ;']) ;
        range3j = max3j - min3j ;
        min3j = min3j - 0.05 * range3j ;
        max3j = max3j + 0.05 * range3j ;

        subplot(nev,nev,nev*(i-1)+j) ;
        eval(['plot(pc' jstr ',pc' istr ',''r+'') ;']) ;
          title('2D Loading plot') ;
          eval(['xlabel(''pc' jstr ''') ;']) ;
          eval(['ylabel(''pc' istr ''') ;']) ;
          axis([min3j, max3j, min3i, max3i]) ;

      end ;


    end ;

  end ;




  if length(ipsfname) > 0 ;     %  then create postscript file with given name
      orient landscape ;
    eval(['print -dpsc \matlab\steve\ps\' ipsfname 'd.ps']) ;
  end ;


end ;



