function SpaSiZer1(data,vlambda,vtheta,vconpar,xgpar,kerstr,vvar,outstr,labelca) 
% SPASIZER, SPAtial smoothing version of SiZer
%     Assess significance of features in 1d Spatial Smoothing
%
%     Creates a series of images,
%     (as postscript files, or mpeg movies)
%     with family of smooths at top, colored SiZer map at bottom
%
%   Can use first 3, 4, 5, 6, 7, 8 or 9 arguments.
% Inputs:
%     data   - either n x 2 matrix of regression data
%                  or n x 1 column vector of time series data
%   vlambda  - vector of lambda tuning parameters
%                  1 entry:    use that as single lambda
%                  2 entries:  use as endpoints for logarthmic grid of:
%                                  5  for vconpar(1) = 1
%                                  25 for vconpar(1) = 2
%                  3 entries:  logarithmic grid, with
%                                  vlambda(1) = lower endpoint
%                                  vlambda(2) = upper endpoint
%                                  vlambda(3) = number of grid points
%                  >3 entries: use those lambdas
%   vtheta  - vector of theta tuning parameters
%                  1 entry:    use that as single theta
%                  2 entries:  use as endpoints for logarthmic grid of:
%                                  5  for vconpar(1) = 1
%                                  25 for vconpar(1) = 2
%                  3 entries:  logarithmic grid, with
%                                  vtheta(1) = lower endpoint
%                                  vtheta(2) = upper endpoint
%                                  vtheta(3) = number of grid points
%    vconpar  - vector of control parameters
%                  1st entry:  type of parametrization
%                       1 - images indexed by theta, SiZer by lambda
%                                   (default, when this is not specified)
%                       2 - images indexed by lambda, SiZer by theta
%                  2nd entry:  type of output
%                       1 - generate series of postscript files
%                                   (default, when this is not specified)
%                       2 - generate MPEG movie
%     ??? not supported yet  ???
%                  3nd entry:  pause parameter
%                       0 - generate everything with no pauses
%                                   (default, when vconpar(2) = 2)
%                       1 - pause between images
%                                   (default, when vconpar(2) = 1)
%                  4th entry:  variance parameter
%                       1 - use spatial conditional variance
%                                   (default, when this is not specified)
%                       2 - use linear operator variance
%                  5th entry:  alpha, significance level
%                                   (0.05 is default when not specified)
%                       1 - use spatial conditional variance
%                                   (default, when this is not specified)
%                       2 - use linear operator variance
%    xgpar  - xgrid parameters
%                  1 entry:    use the x-coordinates of the data 
%                                   (default, when this is not specified)
%                  2 entries:  generate equally spaced grid with:
%                                   xgpar(1) = left end
%                                   xgpar(2) = right end
%                                   nxgrid = 100 
%                  3 entries:    generate equally spaced grid with:
%                                   xgpar(1) = left end
%                                   xgpar(2) = right end
%                                   xgpar(3) = nxgrid
%   kerstr  - String with name of function for covariance kernel
%                  Currently available:
%                     'kergauss'  -  Gaussian covariance kernel
%                               (default, when this is not specified)
%                               (also get this using empty string, '')
%     vvar  - vector of relative variances of residuals,
%                  for handling heteroscedastic case.  This is
%                  used only "relatively", since overall variance is
%                  estimated from the data (and this is normalized to 
%                  have trace = n)
%     ???   not supported yet ???
%   outstr  - string with output prefix, including directory path
%                  for vconpar(1) = 1:   will add index, and ".ps"
%                  for vconpar(1) = 2:   will add ".mpg"
%                                (default is 'SpaSiZer1', which puts 
%                                 output in current directory)
%                  use '' to save no output
%  labelca  - label cell array, 2 x 3, for graphics labels
%                         (if not 2 x 3, or empty, use simple defaults)
%                         (for emtpy cells, use simple defaults)
%                         (create using command "cell")
%                  top row:     graphics labels for family plot
%                  bottom row:  graphics labels for SiZer plot
%                  first column:   title
%                  second column:  xlabel
%                  third column:   ylable
%
% Output:
%     Draws a succession of 2 x 1 SiZer type images,
%         and saves each as a postscript file
%     Or makes a movie, and saves as an MPEG file
%

%    Copyright (c) J. S. Marron 1999





%  Set defaults according to number of input arguments
%
if nargin <= 2 ;     %  not enough inputs, give error message
  disp(['!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!']) ;  
  disp(['!!!   Error from SpaSiZer1.m:             !!!']) ;  
  disp(['!!!   Must input a vlambda and a vtheta   !!!']) ;  
  disp(['!!!   Terminating Execution               !!!']) ;  
  disp(['!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!']) ;  
  return ;
end ;

if nargin == 3 ;    %  only 3 arguments input, use default vconpar
  partype = 1 ; 
  outtype = 1 ;
  pausepar = 1 ;
  varpar = 1 ;
  alpha = 0.05 ;
else ;              %  use input values
  if length(vconpar) == 1 ;   %  then need to add defaults
    partype = vconpar(1) ; 
    outtype = 1 ;
    pausepar = 1 ;
    varpar = 1 ;
    alpha = 0.05 ;
  elseif length(vconpar) == 2 ;   %  then need to add defaults
    partype = vconpar(1) ; 
    outtype = vconpar(2) ;
    if outtype == 1 ;    %  then are making .ps files
      pausepar = 1 ;
    elseif outtype == 2 ;    %  then are making a movie
      pausepar = 0 ;
    end ;
    varpar = 1 ;
    alpha = 0.05 ;
  elseif length(vconpar) == 3 ;   %  then need to add default
    partype = vconpar(1) ; 
    outtype = vconpar(2) ;
    pausepar = vconpar(3) ;
    varpar = 1 ;
    alpha = 0.05 ;
  elseif length(vconpar) == 4 ;   %  then need to add default
    partype = vconpar(1) ; 
    outtype = vconpar(2) ;
    pausepar = vconpar(3) ;
    varpar = vconpar(4) ;
    alpha = 0.05 ;
  else ;
    partype = vconpar(1) ; 
    outtype = vconpar(2) ;
    pausepar = vconpar(3) ;
    varpar = vconpar(4) ;
    alpha = vconpar(5) ;
  end ;
end ;


if nargin <= 4 ;    %  at most 4 arguments input, use default xgpar
  ixgpar = 0 ;    %  Default
else ;
  ixgpar = xgpar ;    %  value was specified, so use it
end ;


if nargin <= 5 ;    %  at most 5 arguments input, use default kerstr
  ikerstr = 'kergauss' ;    %  Default
else ;
  ikerstr = kerstr ;    %  value was specified, so use it
end ;


if nargin <= 6 ;    %  at most 6 arguments input, use default vvar
  ivvar = ones(size(data,1),1) ;    %  Default
else ;
  ivvar = vvar ;    %  value was specified, so use it
end ;


if nargin <= 7 ;    %  at most 7 arguments input, use default outstr
  ioutstr = 'SpaSiZer1' ;    %  Default
else ;
  ioutstr = outstr ;    %  value was specified, so use it
end ;


if nargin <= 8 ;    %  at most 8 arguments input, use default labelca
  ilabelca = cell(2,3) ;    %  Default
else ;
  ilabelca = labelca ;    %  value was specified, so use it
end ;





%  Set internal parameters
%
ndat = size(data,1) ;
          %  number of rows of input matrix
nlambda = length(vlambda) ;
ntheta = length(vtheta) ;




%  Do Checks of Inputs, and unpack
%
if size(data,2) == 1 ;     %   then have time series data,
                           %   so add column of time indices
  xdat = (1:ndat)' ;
  ydat = data ;
elseif size(data,2) ~= 2 ;    %  then have problem, so warn and quit
  disp(['!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!']) ;  
  disp(['!!!   Error from SpaSiZer1.m:     !!!']) ;  
  disp(['!!!   Input "data" is not valid   !!!']) ;  
  disp(['!!!   Terminating Execution       !!!']) ;  
  disp(['!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!']) ;  
  return ;
else ;
  xdat = data(:,1) ;
  ydat = data(:,2) ;
end ;


if length(vlambda) == 2 ;  %  then construct appropriate grid
  if outtype == 1 ;
    nlambda = 5 ;
  elseif outtype == 2 ;
    nlambda = 25 ;
  end ;
  ivlambda = logspace(log10(vlambda(1)),log10(vlambda(2)),nlambda) ;
elseif length(vlambda) == 3 ;  %  then construct appropriate grid
  nlambda = vlambda(3) ;
  ivlambda = logspace(log10(vlambda(1)),log10(vlambda(2)),nlambda) ;
else ;
  ivlambda = vlambda ;
end ;


if length(vtheta) == 2 ;  %  then construct appropriate grid
  if outtype == 1 ;
    ntheta = 5 ;
  elseif outtype == 2 ;
    ntheta = 25 ;
  end ;
  ivtheta = logspace(log10(vtheta(1)),log10(vtheta(2)),ntheta) ;
elseif length(vtheta) == 3 ;  %  then construct appropriate grid
  ntheta = vtheta(3) ;
  ivtheta = logspace(log10(vtheta(1)),log10(vtheta(2)),ntheta) ;
else ;
  ivtheta = vtheta ;
end ;


if length(ixgpar) == 1 ;  %  then use data for grid
  xgrid = sort(xdat) ;
          %  sort, so graphics look nice
  nxgrid = ndat ;
elseif length(ixgpar) == 2 ;  %  then generate grid
  nxgrid = 100 ;
  xgrid = linspace(ixgpar(1),ixgpar(2),nxgrid)' ;
elseif length(ixgpar) == 3 ;  %  then generate grid
  nxgrid = ixgpar(3) ;
  xgrid = linspace(ixgpar(1),ixgpar(2),nxgrid)' ;
elseif length(ixgpar) > 3 ;  %  then use input vector as grid
  xgrid = ixgpar ;
  nxgrid = length(xgrid) ;
end ;


if ~ischar(ikerstr) ;    %  Then this is not a string
  disp(['!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!']) ;  
  disp(['!!!   Warning from SpaSiZer1.m:            !!!']) ;  
  disp(['!!!   Input "kerstr" is not a string       !!!']) ;  
  disp(['!!!   Proceeding with Gaussian default     !!!']) ;  
  disp(['!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!']) ;  

  ikerstr = 'kergauss' ;
elseif isempty(ikerstr) ;
  ikerstr = 'kergauss' ;
end ;


if (size(ivvar,2) ~= 1) & (size(ivvar,1) ~= ndat) ;
                               %  Then size of this input is wrong
  disp(['!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!']) ;  
  disp(['!!!   Warning from SpaSiZer1.m:                 !!!']) ;  
  disp(['!!!   Input "vvar" is wrong size,               !!!']) ;  
  disp(['!!!   Proceeding with homoscedastic default     !!!']) ;  
  disp(['!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!']) ;  

  vvar = ones(ndat,1) ;
end ;


if ~ischar(ioutstr) ;    %  Then this is not a string
  disp(['!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!']) ;  
  disp(['!!!   Warning from SpaSiZer1.m:          !!!']) ;  
  disp(['!!!   Input "outstr" is not a string     !!!']) ;  
  disp(['!!!   Proceeding with default            !!!']) ;  
  disp(['!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!']) ;  

  ioutstr = 'SpaSiZer1' ; 
end ;


flagerror = logical(0) ;
if ~iscell(ilabelca) ;    %  then this is not a cell array
  flagerror = logical(1) ;
  ilabelca = cell(2,3) ;
end ;

toptitlebase = ilabelca{1,1} ;
if ~(ischar(toptitlebase) | isempty(toptitlebase)) ;    %  not a valid title
  flagerror = logical(1) ;
  toptitlebase = 'Spatial Smoothing Family, ' ; 
elseif isempty(toptitlebase) ;    %  then use default
  toptitlebase = 'Spatial Smoothing Family, ' ; 
end ;

bottomtitle = ilabelca{2,1} ;
if ~(ischar(bottomtitle) | isempty(bottomtitle)) ;  %  not a valid title
  flagerror = logical(1) ;
  bottomtitle = ['Spatial Smoothing Family, ' date] ; 
elseif isempty(bottomtitle) ;    %  then use default
  bottomtitle = ['Spatial Smoothing SiZer, ' date] ; 
end ;

topxlabel = ilabelca{1,2} ;
if ~(ischar(topxlabel) | isempty(topxlabel)) ;
  flagerror = logical(1) ;
  topxlabel = '' ; 
elseif isempty(topxlabel) ;    %  then use default
  topxlabel = '' ; 
end ;

bottomxlabel = ilabelca{2,2} ;
if ~(ischar(bottomxlabel) | isempty(bottomxlabel)) ;
  flagerror = logical(1) ;
  bottomxlabel = '' ; 
elseif isempty(bottomxlabel) ;    %  then use default
  bottomxlabel = '' ; 
end ;

topylabel = ilabelca{1,3} ;
if ~(ischar(topylabel) | isempty(topylabel)) ;
  flagerror = logical(1) ;
  topylabel = '' ; 
elseif isempty(topylabel) ;    %  then use default
  topylabel = '' ; 
end ;

bottomylabel = ilabelca{2,3} ;
if ~(ischar(bottomylabel) | isempty(bottomylabel)) ;
  flagerror = logical(1) ;
  if partype == 1 ;
    bottomylabel = 'log_{10}(\lambda)' ; 
  elseif partype == 2 ;
    bottomylabel = 'log_{10}(\theta)' ; 
  end ;
elseif isempty(bottomylabel) ;    %  then use default
  if partype == 1 ;
    bottomylabel = 'log_{10}(\lambda)' ; 
  elseif partype == 2 ;
    bottomylabel = 'log_{10}(\theta)' ; 
  end ;
end ;

if flagerror ;
  disp(['!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!']) ;  
  disp(['!!!   Warning from SpaSiZer1.m:   !!!']) ;  
  disp(['!!!   Input "labelca" is bad      !!!']) ;  
  disp(['!!!   Proceeding with default     !!!']) ;  
  disp(['!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!']) ;  
end ;





%  Do Common Calculations
%
mxddifxd = vec2mat(xdat,ndat) - vec2mat(xdat',ndat) ;
          %  matrix of pairwise differences, xdata - xdata
mxgdifxd = vec2mat(xgrid,ndat) - vec2mat(xdat',nxgrid) ;
          %  matrix of pairwise differences, xgrid - xdata






%  Do Main Loop through images
%
outeris = 1 ;
if partype == 1 ;
  outerie = ntheta ;
  scalemin = log10(vlambda(1)) ;
  scalemax = log10(vlambda(nlambda)) ;
elseif partype == 2 ;
  outerie = nlambda ;
  scalemin = log10(vtheta(1)) ;
  scalemax = log10(vtheta(ntheta)) ;
end ;

for outeri = outeris:outerie ;

  if partype == 1 ;

    theta = vtheta(outeri) ;
    disp(['  Working on theta = ' num2str(theta)]) ;

    toptitle = [toptitlebase '\theta = ' num2str(theta)] ;



    %  Do calculations that depend only on theta
    %
    eval(['[r,r0t,r0pt,dr00pp] = ' ikerstr ...
                       '(mxddifxd,mxgdifxd,theta) ;']) ;
          %  calls covariance kernel subroutine, e.g. kergauss.m
          %  Which gets:
          %  r - matrix sig11 / rho 
          %  r0t - matrix sig21 / rho = transpose of R0 matrix
          %                     = transpose of matrix sig12 / rho
          %                               (formerly called mbas)
          %  r0pt - matrix sig31 / rho = transpose of R0p matrix
          %                     = transpose of matrix sig13 / rho
          %                               (formerly called mdbas)
          %  dr00pp = ds33 - diagonal entries of sig33 matrix
          %                               (suitably rescaled)

    [u,md] = eig(r) ;
    vd = diag(md) ;
          %  gives eigen representation of r:   r = u D u'





  elseif partype == 2 ;

    lambda = vlambda(outeri) ;
    disp(['  Working on lambda = ' num2str(lambda)]) ;

    toptitle = [toptitlebase '\lambda = ' num2str(lambda)] ;


  end ;




  %  Setup inner loop through rows of Sizer map
  %
  mfhatgrid = [] ;
  mdfhatgrid = [] ;
  mvardfhg = [] ;
  vgq = [] ;

  inneris = 1 ;
  if partype == 1 ;
    innerie = nlambda ;
  elseif partype == 2 ;
    innerie = ntheta ;
  end ;


  for inneri = inneris:innerie ;


    if partype == 1 ;

      lambda = vlambda(inneri) ;

      disp(['    working on lambda = ' num2str(lambda)]) ;


    elseif partype == 2 ;

      theta = vtheta(inneri) ;

      disp(['    working on theta = ' num2str(theta)]) ;



      %  Do calculations that depend only on theta
      %
      eval(['[r,r0t,r0pt,dr00pp] = ' ikerstr ...
                       '(mxddifxd,mxgdifxd,theta) ;']) ;
          %  calls covariance kernel subroutine, e.g. kergauss.m
          %  Which gets:
          %  r - matrix sig11 / rho 
          %  r0t - matrix sig21 / rho = transpose of R0 matrix
          %                     = transpose of matrix sig12 / rho
          %                               (formerly called mbas)
          %  r0pt - matrix sig31 / rho = transpose of R0p matrix
          %                     = transpose of matrix sig13 / rho
          %                               (formerly called mdbas)
          %  dr00pp = ds33 - diagonal entries of sig33 matrix
          %                               (suitably rescaled)

      [u,md] = eig(r) ;
      vd = diag(md) ;
          %  gives eigen representation of r:   r = u D u'



    end ;






    %  calculate smooth, for family
    %
    dlii = diag(1 ./ (vd + lambda * ones(ndat,1))) ;
          %  inverse part of coefficients
    s11i = u * dlii * u' ;
          %  matrix sig11^(-1) (rescaled by rho)
          %              (formerly called cwoy)


    s11iy = s11i * ydat ;
          %  matrix sig11^(-1) * y
          %            (formerly called c)
    fhatgrid = r0t * s11iy ;
    fhatdat = r * s11iy ;

    mfhatgrid = [mfhatgrid, fhatgrid] ;




    %  get variance estimate
    %
    vresid = ydat - fhatdat ;    
    sig = std(vresid) ;



    %  calculate deriv, and var est, for SiZer
    %
    dfhatgrid = r0pt * s11iy ;

    mdfhatgrid = [mdfhatgrid, dfhatgrid] ;

    if varpar == 1 ;    %  Then do full conditional variance
      ds33 = dr00pp ;
          %  diagonal entries of sig33 matrix (suitably rescaled)
          %  (this matrix is actually exactly r00pp)
      ds31s11is13 = sum((r0pt * s11i) .* r0pt,2) ;
          %  diagonal entries of sig31*sig11^(-1)*sig13 (rescaled)
          %  exressed as sum along rows of product
      vardfhg = (sig^2 / lambda) * (ds33 - ds31s11is13) ;
          %  multiply by rho = sig^2 / lambda, to put on right scale

    elseif varpar == 2 ;    %  Then estimate variance as linear operator
      dhatmat = r0pt * s11i ;
          %  hat matrix for derivative estimation
      vardfhg = (dhatmat .* dhatmat) * (sig^2 * ones(ndat,1)) ;
          %  square of hat matrix times vector of local sig's

    end;

    mvardfhg = [mvardfhg, vardfhg] ;




    %  Set up SiZer quantiles
    %
    edf = sum(vd ./ (vd + lambda)) ;
          %  effective degrees of freedom

    beta = (1 - alpha)^(1/edf) ;
    gquant = -phiinv((1 - beta) / 2) ;

    vgq = [vgq gquant] ;





  end ;    %  of inneri loop  (through rows of SiZer map)




  %  Set up graphics
  %
  mdfhatgrid = mdfhatgrid' ;
  mvardfhg = mvardfhg' ;
  vgq = vgq' ;
          %  put on SiZer coordinates


  %  Put together SiZer map
  %
  if innerie > 1 ;     %  then have more than one row in SiZer plot
    mloci = mdfhatgrid - vec2mat(vgq,nxgrid) .* sqrt(mvardfhg) ;
          %  Lower confidence (simul.) surface for derivative
    mhici = mdfhatgrid + vec2mat(vgq,nxgrid) .* sqrt(mvardfhg) ;
          %  Upper confidence (simul.) surface for derivative
  else ;               %  then only one row in SiZer plot
    mloci = mdfhatgrid - vgq .* sqrt(mvardfhg) ;
          %  Lower confidence (simul.) surface for derivative
    mhici = mdfhatgrid + vgq .* sqrt(mvardfhg) ;
          %  Upper confidence (simul.) surface for derivative

  end ;



  %  Construct "gray level map", really assignment
  %    of pixels to integers, with idea:
  %          1 (very dark)    - Deriv. Sig. > 0 
  %          2 (darker gray)  - Eff. SS < 5
  %          3 (lighter gray) - Eff. SS >= 5, but CI contains 0
  %          4 (very light)   - Deriv. Sig. < 0 

  mapout = 3 * ones(size(mloci)) ;
            %  default is middle lighter gray

  flag = (mloci > 0) ;
            %  matrix of ones where lo ci above 0
  ssflag = sum(sum(flag)) ;
  if ssflag > 0 ;
    mapout(flag) = ones(ssflag,1) ;
            %  put dark grey where significantly positive
  end ;


  flag = (mhici < 0) ;
            %  matrix of ones where hi ci below 0
  ssflag = sum(sum(flag)) ;
  if ssflag > 0 ;
    mapout(flag) = 4 * ones(ssflag,1) ;
            %  put light grey where significantly negative
  end ;


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





  %  do graphics
  %

  %  family on top
  %
  subplot(2,1,1) ;
    plot(xdat,ydat,'go', ...
         xgrid,mfhatgrid,'c-' ...
                             ) ;
      title(toptitle) ;
      xlabel(topxlabel) ;
      ylabel(topylabel) ;





  %  SiZer map on bottom
  %
  subplot(2,1,2) ;
    image([xgrid(1),xgrid(nxgrid)],[scalemin,scalemax],mapout) ;
      set(gca,'YDir','normal') ;
      colormap(cocomap) ;
      xlabel('x') ;
      ylabel('log10(\lambda)') ;
      title(bottomtitle) ;
      xlabel(bottomxlabel) ;
      ylabel(bottomylabel) ;




  if partype == 1 ;      %  then write to a postscript file

    if ~isempty(ioutstr) ;    %  then want to do this

        orient tall ;
      prstr = ['print -dpsc ' ioutstr ...
                                   num2str(outeri) ...
                                            '.ps ;'] ;
      eval(prstr) ;

    end ;

  elseif partype == 2 ;     %  then save movie frame





  end  ;




  if (pausepar == 1) & (outeri ~= outerie) ;
    disp('Hit any key for next plot') ;
    pause ;
  end ;




end ;    %  of outeri loop (through images)





if (partype == 2) ;     %  then run movie, and save as mpeg



  if ~isempty(ioutstr) ;    %  then save as mpeg file


  end ;


end ;




