function makeplot = sizer2(data,paramstruct) 
% SIZER2, SIgnificance of ZERo crossings of derivatives, version 2
%     Much more flexible version of gpanal2, which uses a "Structure"
%     input form of parameters, to allow more choices.  Important
%     new feature is the movie version of SiZer.  Intended for a
%     quick first look at data sets. 
% Inputs:
%   data        - either n x 1 column vector of density estimation data
%                     or n x 2 matrix of regression data:
%                             X's in first column,  Y's in second
%   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,...
%                                             ) ;
%
%                          where any of the following can be used,
%                          these are optional, misspecified values
%                          revert to defaults
%
%    fields            values
%
%    iout             1  (default)  use 4 panels:  family overlay & 
%                                          surface, SiZer, SiCon
%                     2  use 3 panels: family overlay, SiZer, SiCon
%                     3  use 2 panels: family overlay, SiZer
%                     4  family overlay only
%                     5  family surface only
%                     6  SiZer only
%                     7  SiCon only
%
%    ihazard          0  (default)  give ordinary density or 
%                            regression estimates
%                     1  give hazard rate estimates
%                            (note, no SiCon, or data-driven
%                                 bandwidths allowed)
%
%    icensor          0  (default)  assume data are uncensored
%                     1  assume data are censored, and have a
%                            2nd column of censoring indicators 
%                            (works only for density or hazard
%                                 estimation)
%                            (note, no SiCon, or data-driven
%                                 bandwidths allowed)
%
%    imovie           1  (default)  make movie version
%                     0  make a single still plot
%                            (reset to this for iout >= 4)
%
%    icolor           1  (default)  full color version
%                     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')
%
%    famsurtitle      Title for family surface plot
%                           (default is 'SiZer colored Scale Space')
%
%    sizertitle       Title for SiZer plot
%                           (default is 'SiZer Map')
%
%    sicontitle       Title for SiCon plot
%                           (default is 'SiCon Map')
%
%    titlefontsize    Font Size for titles (uses Matlab default)
%                                   (18 is "fairly large")
%
%    viewangle        Angle for Viewing of 3d surface plot, 
%                         in degrees azimuth and elevation, 
%                         recommend  [165,30] (default)
%                                or  [195,30]
%
%    ndataoverlay     0  no data plot
%                     1  (default) plot 1000 points (random choice)
%                     2  plot full data set
%                     n > 2  plot n random points
%
%    iscreenwrite     0  (default)  no screen writes
%                     1  write to screen to show progress
%
%    nrepeat          number of times to repeat movie (default = 5)
%
%    nbin             number of bins (default = 401)
%
%    minx             left end of bin range (default is min of data)
%
%    maxx             right end of bin range (default is max of data)
%
%    bpar             bin range boundary handling parameter
%                           0 - (default), move data to ends
%                           1 - truncate data outside ends
%
%    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/SiCon
%                           (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)
%
%    ihhighlight      0 don't highlight data driven h (in static output)
%                     1 (default) highlight data driven h
%
% 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:
%    bwsjpib.m
%    bwrswb.m
%    gpkde.m
%    gpnpr.m
%    gpsz1.m
%    gpsc1.m

%    Copyright (c) J. S. Marron 2000



%  First set all parameter to defaults
%
iout = 1 ;
ihazard = 0 ;
icensor = 0 ;
imovie = 1 ;
icolor = 1 ;
savestr = [] ;
ndataoverlay = 1 ;
iscreenwrite = 0 ;
nbin = 401 ;
minx = [] ;
maxx = [] ;
bpar = 0 ;
nfh = [] ;
fhmin = [] ;
fhmax = [] ;
nsh = [] ;
shmin = [] ;
shmax = [] ;
ihhighlight = 1 ;
nrepeat = 5 ;
viewangle = [165,30] ;
xlabelstr = [] ;
ylabelstr = [] ;
labelfontsize = [] ;
famoltitle = ['Family Overlay, ' date] ;
famsurtitle = ['SiZer colored Scale Space'] ;
sizertitle = ['SiZer Map'] ;
sicontitle = ['SiCon Map'] ;
titlefontsize = [] ;





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

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

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

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

  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 sizer2.m:     !!!') ;
      disp('!!!   Invalid savestr,           !!!') ;
      disp('!!!   using default of no save   !!!') ;
      disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
      savestr = [] ;
    end ;
  end ;

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

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

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

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

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

  if isfield(paramstruct,'bpar') ;    %  then change to input value
    bpar = getfield(paramstruct,'bpar') ; 
  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,'ihhighlight') ;    %  then change to input value
    ihhighlight = getfield(paramstruct,'ihhighlight') ; 
  end ;

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

  if isfield(paramstruct,'nrepeat') ;    %  then change to input value
    nrepeat = getfield(paramstruct,'nrepeat') ; 
  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,'famsurtitle') ;    %  then change to input value
    famsurtitle = getfield(paramstruct,'famsurtitle') ; 
  end ;

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

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

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


end ;  %  of resetting of input parameters






%  Turn iout into control parameters
%
if iout == 1 ;
  viplot = [1; 1; 1; 1] ;
         %  indicators for family overlay, family surface, SiZer, SiCon
elseif iout == 2 ;
  viplot = [1; 0; 1; 1] ;
elseif iout == 3 ;
  viplot = [1; 0; 1; 0] ;
elseif iout == 4 ;
  viplot = [1; 0; 0; 0] ;
elseif iout == 5 ;
  viplot = [0; 1; 0; 0] ;
elseif iout == 6 ;
  viplot = [0; 0; 1; 0] ;
elseif iout == 7 ;
  viplot = [0; 0; 0; 1] ;
end ;
nplot = sum(viplot) ;

if  nplot == 1  &  imovie == 1  ;
          %  then reset imovie, since it focuses on figures, not axes 
  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
  disp('!!!   Warning from sizer2.m:     !!!') ;
  disp('!!!   resetting imovie to 0      !!!') ;
  disp('!!!   (use iout < 4 for movie)   !!!') ;
  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
  imovie = 0
end ;

ihazcen = (ihazard == 1  |  icensor == 1) ;
          %  one if are doing either hazard estimation, or censoring
if  ihazcen == 1   &   viplot(4) == 1 ;
  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
  disp('!!!   Error from sizer2.m:       !!!') ;
  disp('!!!   SiCon not supported for    !!!') ;
  disp('!!!   Hazard Est. or Censoring   !!!') ;
  disp('!!!   Choose another iout        !!!') ;
  disp('!!!   Terminating Execution      !!!') ;
  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
  return ;  
end ;


%  set movie parameters
%
nrif = 4 ;
          %  number of times to Repeat Individual Movie Frames
          %    (controls speed of .mpeg version)
movopt = 1 ;
%  Movie options, currently using same as Cornean
%
%  This is MATLAB "help" to use for choosing parameters:
%	MPGWRITE(M, map, 'filename', options) Encodes M in MPEG format
% 	using the specified colormap and writes the result to the
% 	specified file.  The options argument is an optional vector of
% 	8 or fewer options where each value has the following meaning:
% 	1. REPEAT:
% 		An integer number of times to repeat the movie
% 		(default is 1).
% 	2. P-SEARCH ALGORITHM:
% 		0 = logarithmic	(fastest, default value)
% 		1 = subsample
% 		2 = exhaustive	(better, but slow)
% 	3. B-SEARCH ALGORITHM:
% 		0 = simple	(fastest)
% 		1 = cross2	(slightly slower, default value)
% 		2 = exhaustive	(very slow)
% 	4. REFERENCE FRAME:
% 		0 = original	(faster, default)
% 		1 = decoded	(slower, but results in better quality)
% 	5. RANGE IN PIXELS:
% 		An integer search radius.  Default is 10.
% 	6. I-FRAME Q-SCALE:
% 		An integer between 1 and 31.  Default is 8.
% 	7. P-FRAME Q-SCALE:
% 		An integer between 1 and 31.  Default is 10.
% 	8. B-FRAME Q-SCALE:
% 		An integer between 1 and 31.  Default is 25.






%  detect whether density or regression data
%
if icensor == 1 ;    %  then assume data are censored

  xdat = data(:,1) ;
  vdel = data(:,2) ;
  idatyp = 1 ;
  mdat = [xdat, vdel] ;

else ;    %  then assume data are uncensored

  if size(data,2) == 1 ;    %  Then is density or hazard estimation
    xdat = data(:,1) ;
    idatyp = 1 ;
    mdat = xdat ;
  else ;                    %  Then assume regression ;
    xdat = data(:,1) ;
    ydat = data(:,2) ;
    idatyp = 2 ;
    mdat = [xdat,ydat] ;
  end ;

end ;



%  check not all data points are same
%
if std(xdat) == 0 ;     %  then all data points are same, warn and quit
  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
  disp('!!!   Error from sizer2.m:     !!!') ;
  disp('!!!   All x values are same,   !!!') ;
  disp('!!!   Terminating Execution    !!!') ;
  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
  return ;
end ;




%  Determine data ranges
%
if isempty(minx) ;
  minx = min(xdat) ;
end ;

if isempty(maxx) ;
  maxx = max(xdat) ;
end ;

ndat = length(xdat) ;
centx = mean([minx;maxx]) ;




%  Check there are differing data points within range
%
flag = (xdat < maxx) & (xdat > minx) ;
if sum(flag) == 0 ;    %  then no data in range
  errflag = 1 ;
else ;                 %  there are data in range, so check are different
  if std(xdat(flag)) == 0 ;
    errflag = 1 ;
  else ;
    errflag = 0 ;
  end ;
end ;

if errflag == 1 ;     %  write error message, and quit
  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
  disp('!!!   Error from sizer2.m:           !!!') ;
  disp('!!!   no x values in chosen range,   !!!') ;
  disp('!!!   or all x values are the same   !!!') ;
  disp('!!!   Terminating Execution          !!!') ;
  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
  return ;
end;





%  Determine bandwidth ranges
%
range = maxx - minx ;
binw = range / (nbin - 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 SiZer2:               !!!') ;
    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 SiZer2:                  !!!') ;
    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 SiZer2:               !!!') ;
    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 SiZer2:                  !!!') ;
    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 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) ;

  famcolorstr = 'b' ;
          %  Plot most curves in blue
  highlightcolorstr = 'w' ;
          %  Plot highlighted curve in white
  datcolorstr = 'g' ;
          %  Plot overlaid data in green
  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 = 'w' ;
          %  Plot most curves in white
  highlightcolorstr = 'w' ;
          %  Plot highlighted curve in white
  datcolorstr = 'w' ;
          %  Plot overlaid data in white
  cendatcolorstr = [.5, .5, .5] ;
          %  Plot overlaid censored data in gray

end ;





%  Calculate family of smooths (if needed)
%
if  viplot(1) == 1  |  viplot(2) == 1  ;    %  Then will show a family plot

  if iscreenwrite == 1 ;
    disp('    SiZer 2 Working on family') ;
  end ;

  if idatyp == 1 ;    %  doing density or hazard estimation

    if ihazcen == 1 ;      %  then are doing hazard or censoring est.
      [mfam, fxgrid] = CHkde(mdat,fvh,ihazard+1,[minx; maxx; nbin],bpar) ;
          %  argument itype (of CHkde) is (ihazard + 1)
      fbottom = 0 ;
    else ;      %  do ordinary density estimation
      [mfam, fxgrid] = gpkde(mdat,fvh,[minx; maxx; nbin],0,bpar) ;
      fbottom = 0 ;
    end ;

  else ;              %  doing regression

    [mfam, fxgrid] = gpnpr(mdat,fvh,[minx; maxx; nbin],0,1,bpar) ;
    fbottom = min(min(mfam)) ;

  end ;
  ftop = max(max(mfam)) ;
  frange = ftop - fbottom ;
  fbottom = fbottom - 0.05 * frange ;
  ftop = ftop + 0.05 * frange ;


end ;    %  of family plot construction if block





%  Calculate data driven bandwidths for highlighting (if needed)
%
if  imovie == 0   & ...
    ihhighlight == 1   & ...
    ihazcen == 0  ;

  if idatyp == 1 ;    %  doing density estimation
    hcent = bwsjpib(mdat,[minx; maxx; nbin],0,bpar) ;
  else ;              %  doing regression
    hcent = bwrswb(mdat,0,[minx; maxx; nbin],bpar) ;
  end ;
  if hcent < fhmin ;    %  if data based h below range
    disp('!!!   Warning from SiZer2: databased h below this range   !!!') ;
    hcflag = 0 ;
  elseif hcent > fhmax ;    %  if databased h above this range
    disp('!!!   Warning from SiZer2: databased 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 ;

else ;

  hcflag = 0 ;

end ;




%  Calculate SiZer (if needed)
%
if  viplot(2) == 1  |  viplot(3) == 1  ; 

  if iscreenwrite == 1 ;
    disp('    SiZer 2 Working on SiZer') ;
  end ;

  if ihazcen == 0 ;
    sizermap = gpsz1(mdat,[minx,maxx,nbin],[shmin,shmax,nsh],bpar) ;
  else ;
    if ihazard == 1 ;
      itype = 2 ;
          %  do hazard estimation
    else ;
      itype = 1 ;
          %  do density estimation
    end ;
    sizermap = CHsz1(mdat,itype,[minx,maxx,nbin],[shmin,shmax,nsh],bpar) ;
  end ;

end ;




%  Do family overlay graphics (if needed) 
%
if viplot(1) == 1 ;

  if nplot > 1 ;    %  then doing multiple graphics
    clf ;   
    fighand = gcf ;
    if nplot == 4 ;
      famolh = subplot(2,2,1) ;
    else ;
      famolh = subplot(nplot,1,1) ;
    end ;
  end ;


  %  Plot family
  %
  vfamh = plot(fxgrid,mfam,famcolorstr) ;
          %  vfamh = vector of family curve handles     
    vax = [minx,maxx,fbottom,ftop] ;
    axis(vax) ;
    th = title(famoltitle) ;
    if ~isempty(titlefontsize) ;
      set(th,'FontSize',titlefontsize) ;
    end ;
    lxh = xlabel(xlabelstr) ;
    lyh = ylabel(ylabelstr) ;
    if ~isempty(labelfontsize) ;
      set(lxh,'FontSize',labelfontsize) ;
      set(lyh,'FontSize',labelfontsize) ;
    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 red color for one in middle
  end ;


  %  Overlay data (if needed)
  %
  if ndataoverlay > 0 ;    %  then overlay a plot of the raw data

    seed = 34593874 ;
    rand('seed',seed) ;

    if ndataoverlay == 1 ;
      ndo = min(ndat,1000) ;
    elseif ndataoverlay == 2 ;
      ndo = ndat ;
    else ;
      ndo = min(ndat,ndataoverlay) ;
    end ;


    if ndo < ndat ;    %  then need to subsample
      [temp,randperm] = sort(rand(ndat,1)) ;
            %  randperm is a random permutation of 1,2,...,ndat
      vkeep = randperm(1:ndo) ;
            %  indices of elements to keep for display
      mdatol = mdat(vkeep,:) ;
      if icensor == 1 ;
        vdelol = logical(vdel(vkeep)) ;
      end ;
    else ;    %  overlay full data set
      mdatol = mdat ;
      if icensor == 1 ;
        vdelol = logical(vdel) ;
      end ;
    end ;


    %  overlay selected data
    %
    hold on ;
      if idatyp == 1 ;    %  doing density estimation
        yrand = fbottom + (0.7 + 0.2 * rand(ndo,1)) ...
                                              * (ftop - fbottom) ;
               %  y coords for jittering
        if icensor == 1 ;
          plot(mdatol(vdelol,1),yrand(vdelol),[datcolorstr '.']) ;
          plot(mdatol(~vdelol,1),yrand(~vdelol),[cendatcolorstr '.']) ;
        else ;
          plot(mdatol,yrand,[datcolorstr '.']) ;
        end ;
      else ;              %  doing regression
        plot(mdatol(:,1),mdatol(:,2),[datcolorstr '.']) ;
      end ;
    hold off ;

  end ;


end ;








%  Do family surface graphics (if needed) 
%
if viplot(2) == 1 ;

  if nsh == nfh ;    % then go ahead with surface graphics

    if nplot > 1 ;    %  then doing multiple graphics
      famsurh = subplot(2,2,2) ;
            %  this only appears in the 4 panel plot
    end ;

    %  First do decimation if needed
    %
    if nbin > 41 ;

      decimfact = ceil((nbin - 1) / 40) ;
            %  factor to decimate by
      keepflag = (1:decimfact:nbin)' ;

      mfamdecim = mfam(keepflag,:) ;
      fxgriddecim = fxgrid(keepflag) ;
      sizermapdecim = sizermap(:,keepflag) ;

    else ;

      mfamdecim = mfam ;
      fxgriddecim = fxgrid ;
      sizermapdecim = sizermap ;

    end ;


    l10hdecim = log10(fvh) ;
    if nfh > 21 ;

      decimfact = ceil((nfh - 1) / 20) ;
            %  factor to decimate by
      keepflag = (1:decimfact:nfh)' ;

      mfamdecim = mfamdecim(:,keepflag) ;
      sizermapdecim = sizermapdecim(keepflag,:) ;
      l10hdecim = l10hdecim(keepflag) ;

    end ;


    nrowsmd = size(sizermapdecim,1) ;
    ncolsmd = size(sizermapdecim,2) ;
    sizermapdecim = sizermapdecim(2:nrowsmd,2:ncolsmd) ;
          %  cut off one row and one column, to give number
          %  of panels

    %  Make surface plot
    %
    surf(fxgriddecim,l10hdecim,mfamdecim',sizermapdecim) ;
          %  surface plot, using SiZer colors
      th = title(famsurtitle) ;
      if ~isempty(titlefontsize) ;
        set(th,'FontSize',titlefontsize) ;
      end ;
      lxh = xlabel(xlabelstr) ;
      lyh = ylabel('log10(h)') ;
      lzh = zlabel(ylabelstr) ;
      if ~isempty(labelfontsize) ;
        set(lxh,'FontSize',labelfontsize) ;
        set(lyh,'FontSize',labelfontsize) ;
        set(lzh,'FontSize',labelfontsize) ;
      end ;

      vax = [minx,maxx,log10(shmin),log10(shmax),fbottom,ftop] ;
      axis(vax) ;
      set(gca,'Xgrid','off') ;
      set(gca,'Ygrid','off') ;
      set(gca,'Zgrid','off') ;
      set(gca,'Xdir','reverse')
      view(viewangle) ; 

      caxis([1, 8]) ;
      if icolor ~= 0 ;     %  Then go for nice colors in sizer and sicom
        colormap(cocomap) ;  
      else ;     %  Then use gray scale maps everywhere
        colormap(comap) ;  
      end ;

    else ;

      disp(['!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!']) ;
      disp(['!!!   Error from SiZer2.m:    !!!']) ;
      disp(['!!!   For surface plot,       !!!']) ;
      disp(['!!!   must have nsh = nfh     !!!']) ;
      disp(['!!!   Terminating execution   !!!']) ;
      disp(['!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!']) ;
      return ;

    end ;


  end ;    %  of family surface plot construction if block






%  Plot SiZer (if needed)
%
if viplot(3) == 1 ;    %  Then will show a SiZer map

  if nplot > 1 ;    %  then doing multiple graphics
    if nplot == 4 ;
      sizerh = subplot(2,2,3) ;
    else ;
      sizerh = subplot(nplot,1,2) ;
    end ;
  end ;


  image([minx,maxx],[log10(shmin),log10(shmax)],sizermap) ;
    set(gca,'YDir','normal') ;
    th = title(sizertitle) ;
    if ~isempty(titlefontsize) ;
      set(th,'FontSize',titlefontsize) ;
    end ;
    lxh = xlabel(xlabelstr) ;
    lyh = ylabel('log10(h)') ;
    if ~isempty(labelfontsize) ;
      set(lxh,'FontSize',labelfontsize) ;
      set(lyh,'FontSize',labelfontsize) ;
    end ;
          %  next add marking lines
    if icolor ~= 0 ;     %  Then go for nice colors in sizer and sicom
      colormap(cocomap) ;  
    else ;     %  Then use gray scale maps everywhere
      colormap(comap) ;  
    end ;
    hold on ;
      if hcflag == 1 ;
        plot([minx; maxx], ones(2,1)*log10(hcent),['-' highlightcolorstr]) ;
      end ;
      plot(centx + 2*svh, log10(svh),':k') ;
      plot(centx - 2*svh, log10(svh),':k') ;
    hold off ;


end ;    %  of SiZer plot construction if block






%  Calculate and plot SiCon (if needed)
%
if viplot(4) == 1 ;    %  Then will show a SiCon map

  if iscreenwrite == 1 ;
    disp('    SiZer 2 Working on SiCon') ;
  end ;

  if nplot > 1 ;    %  then doing multiple graphics
    if nplot == 4 ;
      siconh = subplot(2,2,4) ;
    else ;
      siconh = subplot(nplot,1,3) ;
    end ;
  end ;


  siconmap = gpsc1(mdat,[minx,maxx,nbin],[shmin,shmax,nsh],bpar) ;
          %  get matrix of map values, so can handle both color sets at once
  if icolor ~= 0 ;     %  Then need to adjust colors in sicom, to fit
                        %  joint color range
    siconmap = siconmap + 4 ;
  end ;
  image([minx,maxx],[log10(shmin),log10(shmax)],siconmap) ;
    set(gca,'YDir','normal') ;
    if icolor ~= 0 ;     %  Then go for nice colors in sizer and sicom
      colormap(cocomap) ;  
    else ;     %  Then use gray scale maps everywhere
      colormap(comap) ;  
    end ;
    th = title(sicontitle) ;
    if ~isempty(titlefontsize) ;
      set(th,'FontSize',titlefontsize) ;
    end ;
    lxh = xlabel(xlabelstr) ;
    lyh = ylabel('log10(h)') ;
    if ~isempty(labelfontsize) ;
      set(lxh,'FontSize',labelfontsize) ;
      set(lyh,'FontSize',labelfontsize) ;
    end ;
          %  next add marking lines
    hold on ;
      if hcflag == 1 ;
        plot([minx; maxx], ones(2,1)*log10(hcent),['-' highlightcolorstr]) ;
      end ;
      plot(centx + 2*svh, log10(svh),':k') ;
      plot(centx - 2*svh, log10(svh),':k') ;
    hold off ;


end ;    %  of SiCon plot construction if block








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

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



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

  %  Do highlight on family surface (if needed)
  %
  if viplot(2) == 1 ;
    axes(famsurh) ;
    hold on ;
      famsurlineh = plot3(fxgrid,log10(fvh(1))*ones(nbin,1), ...
                                   mfam(:,1),highlightcolorstr) ;
      set(famsurlineh,'LineWidth',2) ;
          %  use fatter line width for data based choice
    hold off ;
  end ;

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

  %  Add highlight to SiCon map (if needed)
  %
  if viplot(4) == 1 ;
    axes(siconh) ;
          %  make SiCon axes current
    hold on ;
      siconlineh = plot([minx; maxx], ones(2,1)*log10(fvh(1)), ...
                                           ['-'highlightcolorstr]) ;
    hold off ;
  end ;


  mmovie(:,1) = getframe(fighand) ;


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

    h = fvh(iframe) ;

    if viplot(1) == 1 ;
      set(famollineh,'YData',mfam(:,iframe)) ;
    end ;


    if viplot(2) == 1 ;
      set(famsurlineh,'YData',log10(fvh(iframe))*ones(nbin,1)) ;
      set(famsurlineh,'ZData',mfam(:,iframe)) ;
    end ;


    if viplot(3) == 1 ;
      set(sizerlineh,'YData',ones(2,1)*log10(h)) ;
    end ;


    if viplot(4) == 1 ;
      set(siconlineh,'YData',ones(2,1)*log10(h)) ;
    end ;


    mmovie(:,iframe) = getframe(fighand) ;

  end ;    %  of iframe loop, to make movie



  %  Reorder frames, and replay
  %
  vmorder = [(1:nfh),((nfh-1):-1:2)] ;
  vmorder = vec2mat(vmorder,nrif) ;
  vmorder = reshape(vmorder,1,size(vmorder,2)*nrif) ;
          %  replicate frame number, NRIF times

  mmovie = mmovie(:,vmorder') ;
          %  reorder frames, to loop back to beginning

  movie(fighand,mmovie,nrepeat) ;
          %  Play movie on screen



end ;





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

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


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

    if  nplot == 1  |  nplot == 4 ;
      orient landscape ;
    else ;
      orient tall ;
    end ;

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

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

    if icolor ~= 0 ;     %  then save in full color
      mpgwrite(mmovie,cocomap,savestr,movopt) ;
    else ;                %  then save in black and white
      mpgwrite(mmovie,comap,savestr,movopt) ;
    end ;

  else ;

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

  end ;


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

end ;



