function [mapout,xgrid] = CHsz1(data,itype,vxgp,vhgp,eptflag,alpha,simflag,llflag) 
% CHSZ1, Censored, Hazard, SIgnificant derivative ZERo crossings,
%     A variation of gpsz1.m
%     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, 6 or 7 arguments.
%
% Inputs:
%     data   - either n x 1 column vector of uncensored data
%                  or n x 2 matrix of censored data, with:
%                      X's in first column,
%                      delta's in second column, with values:
%                          1 - when X is the actual value
%                          0 - when X is the (right) censoring time
%                                (i.e. the actual value is only
%                                      known to be larger)
%    itype   - index of estimation type:
%                  1  density estimation (default)
%                  2  hazard rate estimation
%     vxgp   - vector of x grid parameters:
%                  0 (or not specified)  -  use endpts of data and 401 bins
%                  [le; lr; nb] - le left, re right, and nb bins
%     vhgp   - vector of h (bandwidth) grid parameters:
%                  0 (or not specified)  -  use (2*binwidth) to (range),
%                                                and 21 h's
%                  [hmin; hmax; nh]  -  use hmin to hmax and nh h's.
%    eptflag - endpoint truncation flag (only has effect when imptyp = 0):
%                  0 (or not specified)  -  move data outside range to
%                                   nearest endpoint
%                  1  -  truncate data outside range
%    alpha   -  Usual "level of significance".  I.e. C.I.s have coverage
%                  probability 1 - alpha.  (0.05 when not specified)
%    simflag -  Confidence Interval type (simultaneous vs. ptwise)
%                  0  -  Use Pointwise C.I.'s
%                  1 (or not specified)  -  Use Simultaneous C.I.'s
%    llflag  -  Type of Regression Estimator (LL  vs. NW)
%                  0  -  Use Nadaraya Watson est. (in regression problems)
%                  1 (or not specified)  -  Use Local Linear est.
%                            (has no effect for density est)
%
% Outputs:
%     (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)'  
%
% Assumes path can find personal functions:
%    vec2mat.m
%    CHlbinr.m
%    phiinv.m
%

%    Copyright (c) J. S. Marron 2000


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

if nargin <= 2 ;    %  at most 2 arguments input, use default vxgp
  ivxgp = 0 ;
else ;              %  xgrid was specified, use that
  ivxgp = vxgp ;
end ;

if nargin <= 3 ;    %  at most 3 arguments input, use default vhgp
  ivhgp = 0 ;
else ;              %  hgrid parameters were specified, use them
  ivhgp = vhgp ;
end ;

if nargin <= 4 ;    %  at most 4 arguments input, use default endpt trunc
  ieptflag = 0 ;    %  Default
else ;
  ieptflag = eptflag ;    %  Endpt trunc was specified, so use it
end ;

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

if nargin <= 6 ;    %  at most 6 arguments input, use default simul CI's
  isimflag = 1 ;    %  Default
else ;
  isimflag = simflag ;    %  Endpt trunc was specified, so use it
end ;

if nargin <= 7 ;    %  at most 7 arguments input, use default LL in reg
  illflag = 1 ;    %  Default
else ;
  illflag = llflag ;    %  Endpt trunc was specified, so use it
end ;



%  detect whether censored or uncensored data
%
if size(data,2) == 1 ;    %  Then is uncensored
  xdat = data(:,1) ;
  vdel = ones(length(xdat),1) ;
  icen = 0 ;
else ;                    %  Then assume censored
  xdat = data(:,1) ;
  vdel = data(:,2) ;
  icen = 1 ;
end ;


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

%  Set h grid stuff
%
range = right - left ;
binw = range / (ngrid - 1) ;
if ivhgp == 0 ;   %  then use standard default h grid
  ivhgp = [2 * binw,range,21] ;
end ;
hmin = ivhgp(1) ;
hmax = ivhgp(2) ;
nh = ivhgp(3) ;
if nh == 1 ;
  vh = hmax ;
else ;
  vh = logspace(log10(hmin),log10(hmax),nh) ;
end ;



%  Sort and bin the data
%
[xdat,vsortind] = sort(xdat) ;
if icen == 1 ;    %  then there is some censoring
  vdel = vdel(vsortind) ;
end ;

if iitype == 1 ;    %  then are doing density estimation
  fbar = ones(n,1) ;
else ;    %  then are doing hazard estimation
  fbar = 1 - KMcdf(xdat,vdel,1) ;
end ;

if icen == 0 ;     %  then there is no censoring
  gbar = ones(n,1) ;
else ;     %  then there is censoring
  gbar = 1 - KMcdf(xdat,1 - vdel,1) ;
end ;

lbar = fbar .* gbar ;

bdat = [xdat, vdel, lbar] ;
bincts = CHlbinr(bdat,ivxgp,ieptflag) ;


%  Also get "unadjusted bin counts", to use in ESS
%
if icen == 0 ;
  uabincts = gplbinr(xdat,ivxgp,ieptflag) ;
else ;
  if sum(vdel) > 0 ;    %  then there are some uncensored obs's
    uabincts = gplbinr(xdat(logical(vdel)),ivxgp,ieptflag) ;
  else ;
    uabincts = zeros(ngrid,1) ;
  end ;
end ;


if ieptflag == 1 ;
  n = sum(xdat >= left  &  xdat <= right) ;
          %  put this here in case of truncations during binning
end ;



%  Construct Surfaces
%
mdsurf = [] ;
          %  Derivative surface
mesurf = [] ;
          %  Effective sample size surface
mvsurf = [] ;
          %  Estimated Variance of Derivative
vgq = [] ;
          %  Vector of Gaussian Quantiles (for simultaneous CI's)

%  Create grid for kernel values
delta = (right - left) / (ngrid - 1) ;    %  binwidth
k = ngrid - 1 ;    %  index of last nonzero entry of kernel vector

%  Loop through bandwidths
for ih = 1:nh ;
  h = vh(ih) ;

  %  Set common values
  arg = linspace(0,k * delta / h,k + 1)' ;
  kvec = exp(-(arg.^2) / 2) ;
  kvec = [flipud(kvec(2:k+1)); kvec] ;
        %  construct symmetric kernel


  %  Vector of Effective sample sizes 
  %          (notation "s0" is consistent with below)
  ve = conv(uabincts(:,1),kvec) ;
          %  use unadjusted bin counts here, to count
          %  actual pieces of data
  ve = ve(k+1:k+ngrid) ;
          %  denominator of NW est.
          %    (same as sum for kde)


  kvecd = -arg .* exp(-(arg.^2) / 2) / sqrt(2 * pi) ;
  kvecd = [-flipud(kvecd(2:k+1)); kvecd] ;
        %  construct symmetric kernel


  flag0 = (uabincts == 0) ;
  n0 = sum(flag0) ;
  if n0 > 0 ;
    denom = uabincts ;
    denom(flag0) = ones(n0,1) ;
    twiddlefactor = bincts ./ denom ;
    twiddlefactor = (1 - flag0) .* twiddlefactor ;
  else ;
    twiddlefactor = bincts ./ uabincts ;
  end ;
  cabincts = bincts .* twiddlefactor ;

  vd = conv(bincts,kvecd) ;
  vv = conv(cabincts,kvecd.^2) ;
  vd = vd(k+1:k+ngrid) / (n * h^2) ;
  vv = vv(k+1:k+ngrid) / (n * h^4) ;


  vv = vv - vd.^2 ;
  vv = vv / n ;



  %  Get Gaussian quantile, for CI's
  flag = (ve >= 5) ;
          %  locations where effective sample size >= 5
  if sum(flag) > 0 ;
    if isimflag == 0 ;         %  do pt'wise CI's
      gquant = phiinv(1 - (ialpha / 2)) ;
    else ;                     %  do simultaneous CI's
      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) ;
    end ;
  else ;
    gquant = inf ;
  end ;


  mdsurf = [mdsurf vd] ;
  mesurf = [mesurf ve] ;
  mvsurf = [mvsurf vv] ;
  vgq = [vgq gquant] ;

end ;    %  of ih loop through bandwidths


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

  image([left,right],[log10(hmin),log10(hmax)],mapout) ;
    set(gca,'YDir','normal') ;
    colormap(comap) ;  
    xlabel('x') ;
    ylabel('log10(h)') ;


end ;

