function [mapout,xgrid,vedf] = gpsz1ss(data,vxgp,vspargp,alpha) 
% GPSZ1, General Purpose SiZer for Smoothing Splines
%     Creates gray level map (function of location and bandwidth),
%     which is:
%         "dark"   at points where deriv sig > 0
%         "light"  at points where deriv sig < 0
%         "gray"   at points where deriv roughly 0
%         "darker gray" where "effective sample size" < 5
%   Can use first 1, 2, 3, 4, 5 or 6 arguments.
% Inputs:
%     data   - n x 2 matrix of regression data:
%                      X's in first column,  Y's in second
%     vxgp   - vector of x grid parameters:
%                  0 (or not specified)  -  use endpts of data and 101 bins
%                  [le; lr; nb] - le left, re right, and nb bins
%     vspargp   - vector of sparbda (smoothing) grid parameters:
%                  0 (or not specified)  -  use (2*binwidth) to (range),
%                                                and 21 spar's
%                  [sparmin; sparmax; nspar]  -  use sparmin to sparmax and
%                                             nspar smoothing parameter's.
%    alpha   -  Usual "level of significance".  I.e. C.I.s have coverage
%                  probability 1 - alpha.  (0.05 when not specified)
%
% Output:
%     (none)  -  Draws a gray level map (in the current axes)
%     mapout  -  output of gray level map
%     xgrid   -  col vector grid of points at which estimate(s) are 
%                    evaluted (useful for plotting), unless grid is input,
%                    can also get this from linspace(le,re,nb)'  
%     vedf    -   vector of global effective degrees of freedom, trace(A)
%
% Assumes path can find personal functions:
%    vec2mat.m
%    gplbinr.m
%    phiinv.m
%

%    copied from gpsz1.m
%    Copyright (c) J. S. Marron & Jin=Ting Zhang 2001


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

if nargin <= 2 ;    %  at most 2 arguments input, use default vspargp
  ivspargp = 0 ;
else ;              %  spargrid parameters were specified, use them
  ivspargp = vspargp ;
end ;

if nargin <= 3 ;    %  at most 3 arguments input, use default endpt trunc
  ialpha = 0.05 ;    %  Default
else ;
  ialpha = alpha ;    %  Endpt trunc was specified, so use it
end ;



%  Separate data matrix
%
xdat = data(:,1) ;
ydat = data(:,2) ;

%  Set x grid stuff
%
n = length(xdat) ;
if ivxgp == 0 ;   %  then use standard default x grid
  %ivxgp = [min(xdat),max(xdat),min([length(data(:,1)),401])] ;%% original 101
  ivxgp = [min(xdat),max(xdat),301]
end ;
left = ivxgp(1) ;
right = ivxgp(2) ;
ngrid = ivxgp(3) ;
xgrid = linspace(left,right,ngrid)' ;
          %  row vector to evaluate smooths at
cxgrid = xgrid - mean(xgrid) ;
          %  centered version, gives numerical stability

%  Set spar grid stuff
%
if ivspargp == 0 ;   %  then use standard default h grid
  if vernum >= 6.0 ;
    sparmin = 10^fzero('SSfamfitlim',1,[],nbin,0.8*nbin,[left,right]) ;
    sparmax = 10^fzero('SSfamfitlim',1,[],nbin,2.2,[left,right]) ;
  else ;
    sparmin = 10^fzero('SSfamfitlim',1,[],[],nbin,0.8*nbin,[left,right]) ;
    sparmax = 10^fzero('SSfamfitlim',1,[],[],nbin,2.2,[left,right]) ;
  end ;
  ivspargp = [sparmin,sparmax,21] ;
end ;
sparmin = ivspargp(1) ;
sparmax = ivspargp(2) ;
nspar = ivspargp(3) ;
if nspar == 1 ;
  vspar = sparmax ;
else ;
  vspar = logspace(log10(sparmin),log10(sparmax),nspar) ;
end ;


%  Bin the data
%
%bincts = gplbinr([xdat ydat],ivxgp,ieptflag) ;
%bincts2 = gplbinr([xdat, ydat.^2],ivxgp) ;
%n = sum(bincts(:,1)) ;
          %  put this here in case of truncations during binning



%  Construct Surfaces
%
[xfit,fhat,fsig,dhat,dsig,ledf,spargcv,vedf]=SSfamfit([xdat,ydat],vspar,ngrid,[left,right]) ;
mdsurf = dhat ;
          %  Derivative surface
mesurf = n ./ ledf ;
          %  Effective sample size surface
mvsurf = dsig.^2 ;
          %  Estimated Variance of Derivative
vgq = [] ;
          %  Vector of Gaussian Quantiles (for simultaneous CI's)



%  Loop through smoothing parameters
for ispar = 1:nspar ;

  ve = mesurf(:,ispar) ;

  %  Get Gaussian quantile, for CI's
  flag = (ve >= 5) ;
          %  locations where effective sample size >= 5
  if sum(flag) > 0 ;
    nxbar = mean(ve(flag)) ;
          %  Crude average effective sample size
    numind = n / nxbar ;
          %  Effective number of independent groups
    beta = (1 - ialpha)^(1/numind) ;
    gquant = -phiinv((1 - beta) / 2) ;
  else ;
    gquant = inf ;
  end ;

  vgq = [vgq, gquant] ;

end ;


%  Construct scale space CI surfaces
%
if length(vgq) > 1 ;    %  then have full matrices
     
  mloci = mdsurf - (ones(ngrid,1)*vgq).* sqrt(mvsurf) ;
          %  Lower confidence (simul.) surface for derivative
  mhici = mdsurf + (ones(ngrid,1)*vgq) .* sqrt(mvsurf) ;
          %  Upper confidence (simul.) surface for derivative
else ;       %  have only vectors (since only one h)
  mloci = mdsurf - vgq * sqrt(mvsurf) ;
          %  Lower confidence (simul.) surface for derivative
  mhici = mdsurf + vgq * sqrt(mvsurf) ;
          %  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 ;

flag = (mesurf <= 5) ;
            %  matrix of ones where effective np <= 5 ;
ssflag = sum(sum(flag)) ;
if ssflag > 0 ;

  mapout(flag) = 2 * ones(ssflag,1) ;
            %  put middle darker grey where effective sample size < 5
end ;


%  Transpose for graphics purposes
mapout = mapout' ;       

%  Make small h come out at the bottom
%%mapout = flipud(mapout) ;%% added by Zjt, 7/12/2001
%  (not needed when "specifying axes backwards", as below)

%  Make plots if no numerical output requested
%
if nargout == 0 ;  %  Then no numerical output, but make a plot
                   %  on the current axes

  %  Set up gray level color map
  comap = [.2, .2, .2; ...
           .35, .35, .35; ...
           .5, .5, .5; ...
           .8, .8, .8] ;
  if edflag==0,
    image([left,right],[log10(sparmin),log10(sparmax)],mapout) ;
  elseif edflag==1,  
    image([left,right],vedf,mapout) ;
  end  
   
    set(gca,'YDir','normal') ;
    colormap(comap) ;  
    xlabel('x') ;
    
    if edflag==0,
      ylabel('log10(lambda)') ;
    elseif vdflag==1,
      ylabel('EDF')
    end    


end ;