function makeplot = sizerdepSM(data,autocov,paramstruct) 
% SIZERSM, SIgnificance of ZERo crossings of derivatives, for DEPendent data
%   Steve Marron's matlab function
%     For determining which features in a smooth 
%     (density estimate, or nonparametric regression)
%     are statistically significant.
%     Recommended default version of SiZer, usable quite simply
%     with just a set of data (e.g. a first analysis of a new
%     set of data), but also very flexible with a wide range of 
%     options available.
% Inputs:
%   data        - either n x 1 column vector of time series data
%                      (times are assumed to be 1,2,...,n)
%                      (note n should be of reasonable length,
%                       100-600, since no binning will be done)
%   autocov    - nx1 column vector of autocovariances
%                      for lags 0,1,...,(n-1)
%                      i.e. var is in first entry
%   paramstruct - a Matlab structure of input parameters
%                    Use: "help struct" and "help datatypes" to
%                         learn about these.
%                    Create one, using commands of the form:
%
%       paramstruct = struct('field1',values1,...
%                            'field2',values2,...
%                            'field3',values3) ;
%
%                          where any of the following can be used,
%                          these are optional, misspecified values
%                          revert to defaults
%
%    fields            values
%
%
%    imovie           1  make movie version
%                     0  (default)  make a single still plot
%                            (will force reset to this for iout > 4)
%
%    icolor           1  (default)  full color version of SiZer
%                     0  fully black and white version
%
%    savestr          string controlling saving of output,
%                         either a full path, or a file prefix to
%                         save in matlab's current directory
%                     unspecified:  results only appear on screen
%                     movie version (imovie = 1):
%                         add .mpg and create MPEG file
%                     static version (imovie = 0):  
%                         add .ps, and save as either
%                              color postscript (icolor = 1)
%                         or
%                              black&white postscript (when icolor = 0)
%
%    xlabelstr        String for labeling x axes (default is none)
%
%    ylabelstr        String for labeling x axes (default is none)
%
%    labelfontsize    Font Size for labels (uses Matlab default)
%                                   (15 is "fairly large")
%
%    famoltitle       Title for family overlay plot
%                           (default is 'Family Overlay, date')
%
%    sizertitle       Title for slope SiZer plot
%                           (default is 'Slope SiZer Map')
%
%    titlefontsize    Font Size for titles (uses Matlab default)
%                                   (18 is "fairly large")
%
%    ndataoverlay     overlay raw data as a jitterplot on the family plot
%                                    (requires  imptyp >= 0)
%                                    (only has effect when plot is made here)
%                       0  -  no data plot
%                                    (use this to do just a single binning
%                                     of the data, e.g. for large data sets)
%                       1  -  (default) overlay up to 1000 points 
%                                    (random choice, when more)
%                       2  -  overlay full data set
%                       n > 2  -  overlay n random points
%
%    dolcolor         data overlay color
%                                    (only has effect when plot is made here)
%                           default is 'g'
%
%    ibigdot          0  (default)  use Matlab default for dot sizes
%                     1  force large dot size (useful since some
%                              postscript graphics leave dots too small)
%
%    dolhtseed        seed for random heights used in data overlay jitter plot
%                           default is [] (for using current Matlab seed) 
%                                    (should be an integer with <= 8 digits)
%
%    iscreenwrite     0  (default)  no screen writes
%                     1  write to screen to show progress
%
%    bpar             bin range boundary handling parameter
%                         0 - (default), move data to ends
%                         1 - truncate data outside ends
%
%    ibdryadj         index of boundary adjustment
%                         0 (or not specified)  -  no adjustment
%                         1  -  mirror image adjustment
%                         2  -  circular design
%                                 (only has effect for density estimation)
%
%    alpha            Usual "level of significance".  
%                         I.e. C.I.s have coverage probability 
%                         1 - alpha.  (0.05 when not specified)
%
%    nfh              number of h's for family
%                           (default, 41 for movies,  11 for static)
%
%    fhmin            left end of family h range (default = binwd * 2) 
%
%    fhmax            right end of family h range (default = range)
%
%    nsh              number of h's for SiZer
%                           (default, 41 for movies,  11 for static)
%
%    shmin            left end of SiZer h range (default = binwd * 2) 
%
%    shmax            right end of SiZer h range (default = range)
%
%    moviefps         movie speed, in frames per second (default = 5)
%
%    moviecstr        movie compression string, for type of AVI compression:
%                            most look bad with 256 color adapter,
%                            so use a higher one
%                       'MSVC'
%                            requires 256 color graphic adapter,
%                            streamlines and contours OK (on 1st run),
%                            but dots look bad
%                       'None'   (no compression)
%                            looks good but big file
%                       'Cinepak'   (default)
%                            looks good, small file
%                       'Indeo3'
%                            gives warning about "frame size"
%                            good color, but blurry, small file
%                       'Indeo5'
%                            gives warning about "frame size"
%                            good color, but blurry, small file
%
%    nrepeat          number of times to repeat movie (default = 2)
%
%    ishoweffwind    0  do not show effective window width
%                            (i.e. curves showing +- 2 h)
%                    1  (default) show effective window width
%
%    hhighlight      0  (default) don't highlight any h 
%                                (in static family output)
%                    -1 (default) highlight h closest to data driven
%                    h > 0  highlight h closest to this value
%                          Note: doesn't appear when no family plot
%                                is computed
%
% Outputs:
%     For iout = 1,2,3:   graphics in current Figure
%     For iout = 4,5,6,7:   graphics in current axes
%     When savestr exists,
%        For imovie = 1:  MPEG file saved in 'savestr'.mpg
%        For imovie = 0:  Postscript file saved in 'savestr'.ps
%                        (color postscript for icolor = 1)
%                        (B & W postscript for icolor = 0)
%    
%
% Assumes path can find personal functions:
%    bwsjpiSM.m
%    bwrswSM.m
%    kdeSM.m
%    nprSM.m
%    sz1SM.m
%    sc1SM.m
%    vec2matSM

%    Copyright (c) J. S. Marron, V. Rondonotti 2000-2003



%  First set all parameters to defaults
%
imovie = 0 ;
icolor = 1 ;
savestr = [] ;
xlabelstr = [] ;
ylabelstr = [] ;
labelfontsize = [] ;
famoltitle = ['Family Overlay, ' date] ;
sizertitle = ['Slope SiZer Map'] ;
titlefontsize = [] ;
ndataoverlay = 1 ;
dolcolor = 'g' ;
ibigdot = 0 ;
dolhtseed = [] ;
iscreenwrite = 0 ;
bpar = 0 ;
ibdryadj = 0 ;
alpha = 0.05 ;
nfh = [] ;
fhmin = [] ;
fhmax = [] ;
nsh = [] ;
shmin = [] ;
shmax = [] ;
moviefps = 5 ;
moviecstr = 'Cinepak' ;
nrepeat = 2 ;
ishoweffwind = 1 ;
hhighlight = 0 ;





%  Now update parameters as specified,
%  by parameter structure (if it is used)
%
if nargin > 2 ;   %  then paramstruct is an argument

  if isfield(paramstruct,'imovie') ;    %  then change to input value
    imovie = getfield(paramstruct,'imovie') ; 
  end ;

  if isfield(paramstruct,'icolor') ;    %  then change to input value
    icolor = getfield(paramstruct,'icolor') ; 
  end ;

  if isfield(paramstruct,'savestr') ;    %  then use input value
    savestr = getfield(paramstruct,'savestr') ; 
    if ~ischar(savestr) ;    %  then invalid input, so give warning
      disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
      disp('!!!   Warning from sizerSM.m:    !!!') ;
      disp('!!!   Invalid savestr,           !!!') ;
      disp('!!!   using default of no save   !!!') ;
      disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
      savestr = [] ;
    end ;
  end ;

  if isfield(paramstruct,'xlabelstr') ;    %  then change to input value
    xlabelstr = getfield(paramstruct,'xlabelstr') ; 
  end ;

  if isfield(paramstruct,'ylabelstr') ;    %  then change to input value
    ylabelstr = getfield(paramstruct,'ylabelstr') ; 
  end ;

  if isfield(paramstruct,'labelfontsize') ;    %  then change to input value
    labelfontsize = getfield(paramstruct,'labelfontsize') ; 
  end ;

  if isfield(paramstruct,'famoltitle') ;    %  then change to input value
    famoltitle = getfield(paramstruct,'famoltitle') ; 
  end ;

  if isfield(paramstruct,'sizertitle') ;    %  then change to input value
    sizertitle = getfield(paramstruct,'sizertitle') ; 
  end ;

  if isfield(paramstruct,'titlefontsize') ;    %  then change to input value
    titlefontsize = getfield(paramstruct,'titlefontsize') ; 
  end ;

  if isfield(paramstruct,'ndataoverlay') ;    %  then change to input value
    ndataoverlay = getfield(paramstruct,'ndataoverlay') ; 
  end ;

  if isfield(paramstruct,'dolcolor') ;    %  then change to input value
    dolcolor = getfield(paramstruct,'dolcolor') ; 
  end ;

  if isfield(paramstruct,'ibigdot') ;    %  then change to input value
    ibigdot = getfield(paramstruct,'ibigdot') ; 
  end ;

  if isfield(paramstruct,'dolhtseed') ;    %  then change to input value
    dolhtseed = getfield(paramstruct,'dolhtseed') ; 
  end ;

  if isfield(paramstruct,'iscreenwrite') ;    %  then change to input value
    iscreenwrite = getfield(paramstruct,'iscreenwrite') ; 
  end ;

  if isfield(paramstruct,'bpar') ;    %  then change to input value
    bpar = getfield(paramstruct,'bpar') ; 
  end ;

  if isfield(paramstruct,'ibdryadj') ;    %  then change to input value
    ibdryadj = getfield(paramstruct,'ibdryadj') ; 
  end ;

  if isfield(paramstruct,'alpha') ;    %  then change to input value
    alpha = getfield(paramstruct,'alpha') ; 
  end ;

  if isfield(paramstruct,'nfh') ;    %  then change to input value
    nfh = getfield(paramstruct,'nfh') ; 
  end ;

  if isfield(paramstruct,'fhmin') ;    %  then change to input value
    fhmin = getfield(paramstruct,'fhmin') ; 
  end ;

  if isfield(paramstruct,'fhmax') ;    %  then change to input value
    fhmax = getfield(paramstruct,'fhmax') ; 
  end ;

  if isfield(paramstruct,'nsh') ;    %  then change to input value
    nsh = getfield(paramstruct,'nsh') ; 
  end ;

  if isfield(paramstruct,'shmin') ;    %  then change to input value
    shmin = getfield(paramstruct,'shmin') ; 
  end ;

  if isfield(paramstruct,'shmax') ;    %  then change to input value
    shmax = getfield(paramstruct,'shmax') ; 
  end ;

  if isfield(paramstruct,'moviefps') ;    %  then change to input value
    moviefps = getfield(paramstruct,'moviefps') ; 
  end ;

  if isfield(paramstruct,'moviecstr') ;    %  then change to input value
    moviecstr = getfield(paramstruct,'moviecstr') ; 
  end ;

  if isfield(paramstruct,'nrepeat') ;    %  then change to input value
    nrepeat = getfield(paramstruct,'nrepeat') ; 
  end ;

  if isfield(paramstruct,'ishoweffwind') ;    %  then change to input value
    ishoweffwind = getfield(paramstruct,'ishoweffwind') ; 
  end ;
  
  if isfield(paramstruct,'hhighlight') ;    %  then change to input value
    hhighlight = getfield(paramstruct,'hhighlight') ; 
  end ;


end ;    %  of resetting of input parameters


%  Check that autocov could be an autocovariance
%
if  autocov(1) <= 0  ;
  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
  disp('!!!   Error from sizerdepSM.m:    !!!') ;
  disp('!!!   invalid autocovariance      !!!') ;
  disp('!!!   terminating execution       !!!') ;
  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
  return ;
end ;


fighand = gcf ;



%  Construct usual SiZer data matrix
%
ndat = length(data) ;
mdat = [[1:ndat]',data] ;

minx = 1 ;
maxx = ndat ;

centx = mean([minx;maxx]) ;




%  Determine bandwidth ranges
%
range = maxx - minx ;
binw = range / (ndat - 1) ;


if isempty(nfh) ;
  if imovie == 0 ;
    nfh = 11 ;
  else ;
    nfh = 41 ;
  end ;
end ;

if isempty(fhmin) ;
  fhmin = 2 * binw ;
end ;

if isempty(fhmax) ;
  fhmax = range ;
end ;

if nfh > 1 ;
  if fhmin >= fhmax ;
    disp('!!!  Warning from sizerSM:              !!!') ;
    disp('!!!  nfh > 1, and fhmin >= fhmax        !!!') ;
    disp('!!!  will reset nfh to 1, & use fhmax   !!!') ;
    nfh = 1 ;
    fvh = fhmax ;
  else ;
    fvh = logspace(log10(fhmin),log10(fhmax),nfh) ;
  end ;
else ;
  if fhmin == fhmax ;
    fvh = fhmax ;
  else ;
    disp('!!!  Warning from sizerSM:                 !!!') ;
    disp('!!!  nfh = 1, and fhmin, fhmax different   !!!') ;
    disp('!!!  will use fhmax as single h            !!!') ;
    fvh = fhmax ;
  end ;
end ;

if isempty(nsh) ;
  if imovie == 0 ;
    nsh = 11 ;
  else ;
    nsh = 41 ;
  end ;
end ;

if isempty(shmin) ;
  shmin = 2 * binw ;
end ;

if isempty(shmax) ;
  shmax = range ;
end ;


if nsh > 1 ;
  if shmin >= shmax ;
    disp('!!!  Warning from sizerSM:              !!!') ;
    disp('!!!  nsh > 1, and shmin >= shmax        !!!') ;
    disp('!!!  will reset nsh to 1, & use shmax   !!!') ;
    nsh = 1 ;
    svh = shmax ;
  else ;
    svh = logspace(log10(shmin),log10(shmax),nsh) ;
  end ;
else ;
  if shmin == shmax ;
    svh = shmax ;
  else ;
    disp('!!!  Warning from sizerSM:                 !!!') ;
    disp('!!!  nsh = 1, and shmin, shmax different   !!!') ;
    disp('!!!  will use shmax as single h            !!!') ;
    svh = shmax ;
  end ;
end ;



%  Set up colormap
%
if icolor ~= 0 ;     %  Then go for nice colors in slope and curvature sizer

  %  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) ;

  famcolorstr = 'b' ;
          %  Plot most curves in blue
  highlightcolorstr = 'k' ;
          %  Plot highlighted curve in black
  datcolorstr = dolcolor ;
          %  Plot overlaid data in chosen color
  cendatcolorstr = 'y' ;
          %  Plot overlaid censored data in yellow

else ;     %  Then use gray scale maps everywhere

  %  Set up gray level color map
  comap = [.2, .2, .2; ...
           .35, .35, .35; ...
           .5, .5, .5; ...
           .8, .8, .8] ;
  colormap(comap) ;

  famcolorstr = 'k' ;
          %  Plot most curves in black
  highlightcolorstr = 'k' ;
          %  Plot highlighted curve in black
  datcolorstr = 'k' ;
          %  Plot overlaid data in black
  cendatcolorstr = [.5, .5, .5] ;
          %  Plot overlaid censored data in gray

end ;




%  Bin the data
%
%bincts = lbinrSM(mdat,[minx; maxx; nbin],bpar) ;
%    %  eptflag = bpar
%bincts2 = lbinrSM([mdat(:,1), mdat(:,2).^2],[minx; maxx; nbin],bpar) ;
%    %  eptflag = bpar

bincts = [ones(ndat,1), mdat(:,2)] ;
    %  work with one data point per bin

bincts2 = 1 * (autocov(1) - (bincts(:,1) / 1).^2) ;
    %  kludge for sum of squares = n * (var - (mean)^2),
    %  in the case n = 1.

%  careful about how this interfaces with Vitaliana's code

bincts = [bincts, bincts2] ;



%  Calculate data driven bandwidths for highlighting (if needed)
%
if  imovie == 0   & ...
    hhighlight ~= 0  ;

  if isstr(hhighlight) ;

    disp('!!!   Warning from SiZerSM:  hhighlight must be numeric  !!!') ;
    disp('!!!       will not show highlighted bandwidth            !!!') ;

    hcflag = 0 ;

  else ;

    if  hhighlight < 0 ;    %  then do data driven 
      hcent = bwrswSM(mdat,0,[minx; maxx; ndat],bpar) ;
      hcflag = 1 ;
    else ;
      hcent = hhighlight ;
      hcflag = 1 ;
    end ;

    %  Check range and find bandwidth closest to data driven (or input)
    %      for highlighting (if needed)
    %
    if hcflag == 1 ;
      if hcent < fhmin ;    %  if data based h below range
        disp('!!!   Warning from SiZerSM: highlighted h below this range   !!!') ;
        hcflag = 0 ;
      elseif hcent > fhmax ;    %  if databased h above this range
        disp('!!!   Warning from SiZerSM: highlighted h above this range   !!!') ;
        hcflag = 0 ;
      else ;     %  if data based in range, then get its nearest index
        [temp, hcind] = min(abs(log10(hcent) - log10(fvh))) ;
        hcflag = 1 ; 
      end ;
    end ;

  end ;


else ;

  hcflag = 0 ;

end ;




%  Calculate and plot family of smooths (if needed)
%
if iscreenwrite == 1 ;
  disp('    SiZerSM Working on family') ;
end ;

famolh = subplot(2,1,1) ;

paramstruct = struct('vh',fvh,...
                     'vxgrid',[1; ndat; ndat], ...
                     'eptflag',bpar, ...
                     'titlestr',famoltitle, ...
                     'titlefontsize',titlefontsize, ...
                     'xlabelstr',xlabelstr, ...
                     'ylabelstr',ylabelstr, ...
                     'labelfontsize',labelfontsize, ...
                     'linecolor',famcolorstr, ...
                     'ndataoverlay',ndataoverlay, ...
                     'dolcolor',datcolorstr, ...
                     'ibigdot',ibigdot, ...
                     'iplot',1) ;

if ndataoverlay == 0   ;    %  then use bincounts, to avoid rebinning
  paramstruct = setfield(paramstruct,'imptyp',-1) ;
      %  indicate use of previously binned data
  [mfam, fxgrid] = nprSM(bincts(:,[1 2]),paramstruct) ;
      %  only need first two columns for this
else ;                      %  then use raw data, for overlay
  [mfam, fxgrid] = nprSM(mdat,paramstruct) ;
end ;


vax = axis ;
fbottom = vax(3) ;
ftop = vax(4) ;


%  get handles to lines in family plot
%
vaxchil = get(gca,'Children') ;
vfamh = [] ;
for i = length(vaxchil):-1:1 ;    %  since order is reversed
  if get(vaxchil(i),'LineStyle') == '-' ;
    vfamh = [vfamh; vaxchil(i)] ;
  end ;
end ;

%  Highlight data driven curve (if needed)
%
if hcflag == 1 ;
  set(vfamh(hcind),'LineWidth',2) ;
      %  use fatter line width for data based choice
  set(vfamh(hcind),'Color',highlightcolorstr) ;
      %  use given color for one in middle
end ;




%  Calculate SiZer and plot
%
if iscreenwrite == 1 ;
  disp('    SiZerSM Working on SiZer') ;
end ;

sizerh = subplot(2,1,2) ;

%paramstruct = struct('vxgp',[minx; maxx; nbin], ...
%                     'vhgp',[shmin; shmax; nsh], ...
%                     'eptflag',bpar, ...
%                     'ibdryadj',ibdryadj, ...
%                     'alpha',alpha, ...
%                     'imptyp',-1, ...
%                     'icolor',icolor, ...
%                     'titlestr',sizertitle, ...
%                     'titlefontsize',titlefontsize, ...
%                     'xlabelstr',xlabelstr, ...
%                     'ylabelstr','log_{10}(h)', ...
%                     'labelfontsize',labelfontsize, ...
%                     'iplot',1) ;
%
%sizermap = sz1SM(bincts,paramstruct) ;

%  replace above by lines adapted from sz1SM.m,
%  and from Vitaliana's gpsz1d.m:

%  First set parameters
%
vxgp = [1; ndat; ndat] ;
vhgp = [shmin; shmax; nsh] ;
imptyp = -1 ;
eptflag = bpar ;
titlestr = sizertitle ;
ylabelstr = 'log_{10}(h)' ;
iplot = 1 ;


%  Set x grid stuff
%
left = vxgp(1) ;
right = vxgp(2) ;
ngrid = vxgp(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
%
hmin = vhgp(1) ;
hmax = vhgp(2) ;
nh = vhgp(3) ;
if nh == 1 ;
  vh = hmax ;
else ;
  vh = logspace(log10(hmin),log10(hmax),nh) ;
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
ihprevcalc = 0 ;
    %  indicator of whether have done the full calculation
    %      (not done when all locations have ESS < 5)
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(bincts(:,1),kvec) ;
  ve = ve(k+1:k+ngrid) ;
        %  denominator of NW est.
        %    (same as sum for kde)


% 
%  Comment out my old lines, replaced by Vitaliana's
%
%  flag = (ve < 1) ;
%      %  locations with effective sample size < 1
%  ve(flag) = ones(sum(flag),1) ;
%      %  replace small sample sizes by 1 to avoid 0 divide
%      %  no problem below, since gets grayed out
%  s1 = conv(bincts(:,1) .* cxgrid , kvec) ;
%  s1 = s1(k+1:k+ngrid) ;
%  s2 = conv(bincts(:,1) .* cxgrid.^2 , kvec) ;
%  s2 = s2(k+1:k+ngrid) ;
%  t0 = conv(bincts(:,2),kvec) ;
%  t0 = t0(k+1:k+ngrid) ;
%      %  numerator of NW est.
%
%  xbar = conv(bincts(:,1) .* cxgrid , kvec) ;
%  xbar = xbar(k+1:k+ngrid) ;
%      %  Weighted sum of X_i
%  xbar = xbar ./ ve ;
%      %  Weighted avg of X_i
%
%  t1 = conv(bincts(:,2) .* cxgrid , kvec) ;
%  t1 = t1(k+1:k+ngrid) ;
%
%  numerd = t1 - t0 .* xbar ;
%      %  sum(Y_i * (X_i - xbar)) * W      (weighted cov(Y,X))
%  denomd = s2 - 2 * s1 .* xbar + ve .* xbar.^2 ;
%      %  sum((X_i - xbar)^2) * W      (weighted var(X))
%
%  flag2 = denomd < (10^(-10) * mean(denomd)) ;
%      %  for local linear, need to also flag locations where this
%      %  is small, because could have many observaitons piled up
%      %  at one location
%  ve(flag2) = ones(sum(flag2),1) ;
%      %  also reset these, because could have more than 5 piled up
%      %  at a single point, but nothing else in window
%  flag = flag | flag2 ;
%      %  logical "or", which flags both types of locations
%      %  to avoid denominator problems
%
%  denomd(flag) = ones(sum(flag),1) ;
%      %  this avoids zero divide problems, OK, since grayed out later
%
%  mhat = t0 ./ ve ;
%  vd = numerd ./ denomd ;
%      %  linear term from local linear fit (which est's deriv).
%      %       (sometimes called betahat)

%  Vitaliana's version of this (her mdhat = vd me):
  % kernel vector
  %
n = ngrid ;
kk=n-1;                % suggested maximum: lag fix(n/4), but it does not work 
%  arg=linspace(0,(n-1)/h,n)';         % vector n*1 [(i/h), i=0,...,n-1] 
%        %  checked this is same as my arg
  argh=arg*h;                         % vector n*1 [i, i=0,...,n]  

  kvec0=exp(-(arg.^2)/2)/sqrt(2*pi);  % vector n*1 [K(i/h), i=0,...,n-1, K~N]
  kvec1=kvec0.*argh;                  % vector n*1 [K(i/h)*(i), i=0,...,n-1]
  kvec2=kvec1.*argh;                  % vector n*1 [K(i/h)*(i^2), i=0,...,n-1]
 
  skvec0=kvec0.^2;                    % vector n*1 [(K(i/h))^2, i=0,...,n-1]
  skvec1=skvec0.*argh;                % vector n*1 [(K(i/h))^2*(i), i=0,...,n-1]
  skvec2=skvec1.*argh;                % vector n*1 [(K(i/h))^2*(i^2),i=0,..,n-1]


  % s0,s1,s2,t0,t1,u0,u1,u2
  %
  kvec0=[flipud(kvec0(2:n));kvec0];     % vector (2*n-1)*1 of kernel values

%  KVEC0=[KVEC0 kvec0];
    %  I don't need to keep these matrices

  skvec0=[flipud(skvec0(2:n));skvec0];  % vector (2*n-1)*1 of square kernels

  s0=conv(bincts(:,1),kvec0);         % vector [3*n-2]*1 
  s0=s0(n:2*n-1);                     % vector n*1
%  note: her s0 ~ ve me


%  S0=[S0 s0]; 
    %  I don't need to keep these matrices

  t0=conv(bincts(:,2),kvec0);         % vector [3*n-2]*1 
  t0=t0(n:2*n-1);                     % vector n*1 

  u0=conv(bincts(:,1),skvec0);        % vector [3*n-2]*1
  u0=u0(n:2*n-1);                     % vector n*1 

  kvec1=[-flipud(kvec1(2:n));kvec1];    % vector (2*n-1)*1 of kernel*argh

%  KVEC1=[KVEC1 kvec1];
    %  I don't need to keep these matrices

  skvec1=[-flipud(skvec1(2:n));skvec1]; % vector (2*n-1)*1 of kernel^2*argh

  s1=-conv(bincts(:,1),kvec1);        % vector [3*n-2]*1 
  s1=s1(n:2*n-1);                     % vector n*1 

%  S1=[S1 s1]; 
    %  I don't need to keep these matrices

  t1=-conv(bincts(:,2),kvec1);        % vector [3*n-2]*1 
  t1=t1(n:2*n-1);                     % vector n*1 

  u1=-conv(bincts(:,1),skvec1);       % vector [3*n-2]*1 
  u1=u1(n:2*n-1);                     % vector n*1

  kvec2=[flipud(kvec2(2:n));kvec2];     % vector (2*n-1)*1 of kernel*argh^2
  skvec2=[flipud(skvec2(2:n));skvec2];  % vector (2*n-1)*1 of kernel^2*argh^2

  s2=conv(bincts(:,1),kvec2);         % vector [3*n-2]*1
  s2=s2(n:2*n-1);                     % vector n*1 

  u2=conv(bincts(:,1),skvec2);        % vector [3*n-2]*1 
  u2=u2(n:2*n-1);                     % vector n*1


  % estimation of the regression function and its derivative
  %
  den=s2.*s0-s1.*s1; 

%  DEN=[DEN den];
    %  I don't need to keep these matrices

  mhat=(s2.*t0-s1.*t1)./den;          % vector n*1 of estimated values for m
  mdhat=(s0.*t1-s1.*t0)./den;         % vector n*1 of estimated values for m' 

  vd = mdhat ;
      %  put Vitaliana's quantity in my notation

%  MHAT=[MHAT mhat];                   % it will become a matrix n*nh 
%  MDHAT=[MDHAT mdhat];                % it will become a matrix n*nh 
    %  I don't need to keep these matrices




%  acov=ACOV(:,i);   % fixed gammahat function
  acov = autocov ;
      %  turn my input into Vitaliana's notation


  
  % estimation of part of the numerator of the variance
  %
  g0=acov(1);
  vnum1=g0*((s0.^2).*u2+(s1.^2).*u0-2*s0.*s1.*u1);
%  VNUM1=[VNUM1 vnum1];               % it will become a matrix n*nh
    %  I don't need to keep these matrices


 




%  sig2 = conv(bincts(:,3),kvec) ;
%  sig2 = conv(bincts2,kvec) ;
%  sig2 = sig2(k+1:k+ngrid) ;
%  sig2 = sig2 ./ ve - mhat.^ 2 ;

%sig2 = autocov(1) * ones(ndat,1) ;
%     %  my kludge for early dependent tests
%
%  flag2 = sig2 < (10^(-10) * mean(sig2)) ;
%  ve(flag2) = ones(sum(flag2),1) ;
%      %  also reset these
%  flag = flag | flag2 ;
%      %  logical "or", which flags both types of locations
%      %  to avoid denominator problems
%  sig2(flag) = ones(sum(flag),1) ;
%      %  this avoids zero divide problems, OK, since grayed out later
%
%  rho2 = vd.^2 .* (denomd ./ (sig2 .* ve)) ;
%  sig2res = (1 - rho2) .* sig2 ;
%      %  get the residual variance from local linear reg.
%
%  u0 = conv(bincts(:,1) .* sig2res , kvec.^2) ;
%  u0 = u0(k+1:k+ngrid) ;
%  u1 = conv(bincts(:,1) .* sig2res .* cxgrid , kvec.^2) ;
%  u1 = u1(k+1:k+ngrid) ;
%  u2 = conv(bincts(:,1) .* sig2res .* cxgrid.^2 , kvec.^2) ;
%  u2 = u2(k+1:k+ngrid) ;
%  vv = u2 - 2 * u1 .* xbar + u0 .* xbar.^2 ;
%  vv = vv ./ denomd.^2 ;
%      %  vector of variances of slope est. (from local linear)

%  Vitaliana's version of this (her vmdhat = vv me):
%  VMDHAT=[];        % matrix of estimated variance of the derivative function 
%  GQUANT=[];
    %  I don't need to keep these matrices





%  % loop through bandwidths: 
%  % continuation of the first loop, once gammahat is fixed
%  %
%
%  for j=1:nh;

%    kvec0=KVEC0(:,j);
%    kvec1=KVEC1(:,j);
    %  these are alredy defined above

    % estimation of the variance of the derivative of the regression function
    % 

    % w00, w01, w10, w11
    %
    ww00=[];
    ww11=[];
    ww01=[];
    ww10=[];

    SS0=[];
    SS1=[];

    kvec1f=flipud(kvec1);

    kkvec0=[0;kvec0(1:n-1)];                   % n*1
    kkvec1=[0;kvec1f(1:n-1)];                  % n*1

    KA=kkvec0*acov(3:kk+1)';                   % n*(kk-2)
    KA1=kkvec1*acov(3:kk+1)';                  % n*(kk-2)

    ss0=conv(acov(2:kk+1),kvec0(1:2*n-2));     % (3n-4)*1
    ss1=conv(acov(2:kk+1),kvec1f(1:2*n-2));    % (3n-4)*1


    for k=1:n-kk;

      ss00=ss0(n-k:2*n-(k+1));                
      ss11=ss1(n-k:2*n-(k+1));

      SS0=[SS0;ss00']; 
      SS1=[SS1;ss11'];

    end;   
 

    for k=n-kk+1:n-1; 

      ss00=[ss0(n-k);ss00(1:n-1)]-KA(:,n-k); 
      SS0=[SS0;ss00']; 

      ss11=[ss1(n-k);ss11(1:n-1)]-KA1(:,n-k); 
      SS1=[SS1;ss11']; 
    
    end;


    for k=1:n;

      kk0=kvec0(n-(k-1):2*n-(k+1));
      kk1=kvec1(n-(k-1):2*n-(k+1));

      w00=kk0'*SS0(:,k);
      w11=kk1'*SS1(:,k);
      w01=kk0'*SS1(:,k);
      w10=kk1'*SS0(:,k);

      ww00=[ww00;w00];
      ww11=[ww11;w11];
      ww01=[ww01;w01];
      ww10=[ww10;w10];

    end;

%    den=DEN(:,j);
%    vnum1=VNUM1(:,j);
%    s0=S0(:,j);
%    s1=S1(:,j);
    %  these were defined above

    vnum2=2*((s0.^2).*ww11+(s1.^2).*ww00-s0.*s1.*(ww01+ww10));
    vmdhat=(vnum1+vnum2)./(den.^2); 


    vv = vmdhat ;
      %  put Vitaliana's quantity in my notation





  % effective sample size without taking into account the dependence
  %
  ess=sqrt(2*pi)*s0;                 % n*1
%  ESS=[ESS ess];                     % it will become a matrix n*nh
    %  I don't need to keep these matrices


  % quantities for effective sample size taking into account the dependence
  %
  vybar = acov(1) / n + (2 / n) * ...
            sum((1 - (1:(n-1))' / n) .* acov(2:n)) ;
      %  known cov version, from equation (14) in paper

  nstar=acov(1)/vybar;
%  NSTAR=[NSTAR nstar];
    %  I don't need to keep these matrices

%  nstar=NSTAR(:,i);
       %  already defined above
  
  
  % effective sample size
  % 
%  ESTAR=(nstar/n)*ESS; 
  estar=(nstar/n)*ess; 
      %  my vector version
    




  %  Get quantiles, for CI's
  flag = (ve >= 5) ;
          %  locations where effective sample size >= 5
  if sum(flag) > 0 ;

%  my old version
%
%    nxbar = mean(ve(flag)) ;
%        %  Crude average effective sample size
%    numind = ndat / nxbar ;
%        %  Effective number of independent groups
%    beta = (1 - alpha)^(1/numind) ;
%    gquant = -norminv((1 - beta) / 2) ;

%  vitaliana's version of this:
%  lines 514-519
      mestar=mean(estar(flag));
      m=n/mestar;                  % effective number of independent groups
       
      beta=(1-alpha)^(1/m);
    
%      gquant=-phiinv((1-beta)/2);
      gquant = -norminv((1 - beta) / 2) ;






  else ;

    gquant = inf ;

  end ;


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

  vgq = [vgq gquant] ;
%  vitaliana's version of this:
%    GQUANT=[GQUANT gquant];

end ;



if nh > 1 ;    %  then have full matrices
  mloci = mdsurf - vec2matSM(vgq,ngrid) .* sqrt(mvsurf) ;
          %  Lower confidence (simul.) surface for derivative
  mhici = mdsurf + vec2matSM(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 ;
%  vitaliana's version of this:
%  LOWCI=MDHAT-vec2mat(GQUANT,n).*sqrt(VMDHAT);   % lower confidence surface
%  UPCI=MDHAT+vec2mat(GQUANT,n).*sqrt(VMDHAT);    % upper confidence surface



%  Construct "color 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 purple (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 blue (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 red (light gray) 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 gray where effective sample size < 5
end ;


%  Transpose for graphics purposes
mapout = mapout' ;         



%  Make plots if no numerical output requested
%
%if  nargout == 0  | ...
%      iplot == 1  ;  %  Then make a plot


  if icolor ~= 0 ;     %  Then go for nice colors in sizer and sicon

    %  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) ;

  else ;     %  Then use gray scale maps everywhere

    %  Set up gray level color map
    comap = [.2, .2, .2; ...
             .35, .35, .35; ...
             .5, .5, .5; ...
             .8, .8, .8] ;
    colormap(comap) ;

  end ;


  image([left,right],[log10(hmin),log10(hmax)],mapout) ;
    set(gca,'YDir','normal') ;


  if ~isempty(titlestr) ;
    if isempty(titlefontsize) ;
      title(titlestr) ;
    else ;
      title(titlestr,'FontSize',titlefontsize) ;
    end ;
  end ;


  if ~isempty(xlabelstr) ;
    if isempty(labelfontsize) ;
      xlabel(xlabelstr) ;
    else ;
      xlabel(xlabelstr,'FontSize',labelfontsize) ;
    end ;
  end ;


  if ~isempty(ylabelstr) ;
    if isempty(labelfontsize) ;
      ylabel(ylabelstr) ;
    else ;
      ylabel(ylabelstr,'FontSize',labelfontsize) ;
    end ;
  end ;

%
%end ;

%  end of lines from sz1SM.m










%  Highlight data driven curve
%
if hcflag == 1 ;
  hold on ;
    sizerlineh = plot([minx; maxx], ones(2,1)*log10(hcent), ...
                                       ['-' highlightcolorstr]) ;
  hold off ;
end ;


if ishoweffwind == 1 ;    %  then add curves showing effective window width

  hold on ;
    plot(centx + 2*svh, log10(svh),':w') ;
    plot(centx - 2*svh, log10(svh),':w') ;
  hold off ;

end ;




%  Do main movie construction
%
if imovie == 1 ;

  if iscreenwrite == 1 ;
    disp('    SiZerSM Working on Movie') ;
  end ;


  clear moviestruct ;


  %  Do highlight on family plot
  %
  axes(famolh) ;
  hold on ;
    famollineh = plot(fxgrid,mfam(:,1),highlightcolorstr) ;
    set(famollineh,'LineWidth',2) ;
        %  use fatter line width for data based choice
  hold off ;

  %  Add highlight to SiZer map
  %
  axes(sizerh) ;
      %  make SiZer axes current
  hold on ;
    sizerlineh = plot([minx; maxx], ones(2,1)*log10(fvh(1)), ...
                                         ['-' highlightcolorstr]) ;
  hold off ;


  moviestruct(1) = getframe(fighand) ;


  %  Loop through and change
  %
  for iframe = 2:nfh ;

    h = fvh(iframe) ;

    set(famollineh,'YData',mfam(:,iframe)) ;

    set(sizerlineh,'YData',ones(2,1)*log10(h)) ;

		moviestruct(iframe) = getframe(fighand) ;

  end ;    %  of iframe loop, to make movie



  %  Reorder frames, and replay
  %
  vmorder = [(1:nfh),((nfh-1):-1:2)] ;
  moviestruct = moviestruct(vmorder) ;
          %  reorder frames, to loop back to beginning


  movie(fighand,moviestruct,nrepeat,moviefps) ;
          %  Play movie on screen


end ;





%  Save results (if needed)
%
if ~isempty(savestr) ;     %  then save results

  if iscreenwrite == 1 ;
    disp('    SiZerSM saving results') ;
  end ;


  if imovie == 0 ;     %  then save as postscript file

    orient landscape ;

    if icolor ~= 0 ;     %  then make color postscript
      print('-dpsc',[savestr '.ps']) ;
    else ;                %  then make black and white
      print('-dps',[savestr '.ps']) ;
    end ;

  elseif imovie == 1 ;    %  then save as mpeg file

    movie2avi(moviestruct,savestr,'compression',moviecstr, ...
                          'keyframe',moviefps,'fps',moviefps) ;

  else ;

    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
    disp('!!!   Error from sizerSM.m:                !!!') ;
    disp('!!!   Invalid value of imovie,             !!!') ;
    disp('!!!   Terminating Execution without save   !!!') ;
    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
    return ;

  end ;


  if iscreenwrite == 1 ;
    disp('    SiZerSM finished save') ;
    disp('  ') ;
  end ;

end ;



