function gpcd1(data,vax,ncol,psfname) 
% GPCD1, General Purpose Curves as Data analysis
%     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"
%     vax     - vector of axis limits,  in order:
%                   main lower, main upper, Eigv lower, Eigv upper
%                   (use 0 to choose these from the data)
%     ncol    - 2  use 2 column format, eigvecs, resid
%               3, use 3 column format, eigvecs, mean +-, resid's
%               4, use 4 column format, add R2 histo, and kde's (default)
%     psfname - string with prefix for name of postscript file
%                    (will automatically add ".ps")
%                    (will be a color postscript file when icolor = 1)
% Outputs:
%     Only graphics, in current Figure, unless a postscript file is created
%
% Assumes path can find personal functions:
%    bwsjpib.m
%    gpkde.m

%    Copyright (c) J. S. Marron 1997



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

if nargin <= 2 ;    %  at most 2 arguments input, use default ncol
  incol = 4 ;
else ;              %  ncol was specified, use that
  incol = ncol ;
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 ;





%  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


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"
[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
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 ;



%  make output plot
%
npx = incol ;
if incol == 2 ;
  npy = nev + 2 ;
else ;
  npy = nev + 1 ;
end ;
xgrid = (1:d)' ;

subplot(npy,npx,1) ;
  plot(xgrid,data) ;
    title('Raw Data') ;
    mlvax = axis ;
    if length(ivax) == 1 ;     %  then use matlab's axes
      min1 = mlvax(3) ;
      max1 = mlvax(4) ;
    else ;                     %  then use input axis values
      min1 = ivax(1) ;
      max1 = ivax(2) ;
    end ;
    axis([1,d,min1,max1]) ;

if npx > 2 ;
  subplot(npy,npx,2) ;
else ;
  subplot(npy,npx,3) ;
end ;
  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']) ;


if npx > 2 ;
  subplot(npy,npx,3) ;
else ;
  subplot(npy,npx,4) ;
end ;
  plot(xgrid,mresid) ;
    title('Mean Resid.') ;
    mlvax = axis ;
    if length(ivax) == 1 ;     %  then use matlab's axes
      min2 = mlvax(3) ;
      max2 = mlvax(4) ;
    else ;                     %  then use input axis values
      min2 = ivax(3) ;
      max2 = ivax(4) ;
    end ;
    axis([1,d,min2,max2]) ;



if npx == 2  |  npx == 4 ;
  if npx == 2 ;
    subplot(npy,npx,2) ;
  else ;
    subplot(npy,npx,4) ;
  end ;
    bar(1:nev,vr2) ;
      title('PC R-squareds') ;
        xlabel('PC #') ;
        ylabel('R^2') ;
end ;


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

  istr = num2str(iev) ;

  %  do projection plot
  if npx == 2 ;
    subplot(npy,npx,2 * iev + 3) ;
  elseif npx == 3 ;
    subplot(npy,npx,3 * iev + 1) ;
  elseif npx == 4 ;
    subplot(npy,npx,4 * iev + 1) ;
  end ;
  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 = min1 + .9 * (max2 - min2) ;
      text(1.5,ty,[num2str(100 * vr2(iev)) '% of MR']) ;
    end ;


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


  %  do mean +- plot
  if npx > 2 ;
    if npx == 3 ;
      subplot(npy,npx,3 * iev + 2) ;
    elseif npx == 4 ;
      subplot(npy,npx,4 * iev + 2) ;
    end ;
    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']) ;
  end ;


  %  do kde of projections
  if npx == 4 ;
    subplot(npy,npx,4 * iev + 4) ;
    eval(['gpkde(pc' istr ') ;']) ;
      title(['kde, PC' istr]) ;
  end ;


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 ;



