function makeplot = gpszj(data,icolor,vrange,nhp,psfname,idatap,vlhp) 
% GPszj, General Purpose SiZer for Jumps
%     Makes a single page with family analysis, SiZer 1, and a
%     linear bandwidth scale SiZer
%     for studying data with suspected jumps
%   Can use first 1, 2, 3, 4, 5 or 6 arguments.
% 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
%     icolor  - 0  use gray levels
%               1 (or not specified) do red, blue purple
%     vrange  - 0 (or not specified) use the min and max of the data
%               2 vector:  use minx as vrange(1) and maxx as vrange(2)
%               3 vector:  use minx as vrange(1) and maxx as vrange(2)
%                               and number of grid points as vrange(3)
%     nhp     - 0 use default value of 11
%               otherwise, value is number in family, and rows in SiZer
%     psfname - string with prefix for name of postscript file
%                    (will automatically add ".ps")
%                    (will be a color postscript file when icolor = 1)
%     idatap  - 0  don't add data to family plot
%               1  (or not specified) add data to family plot, as
%                    jitter plot (density estimation) or scatterplot
%     vlhp    - vector of log h parameters
%               0  (or not specified) use lower half of SiZer range
%               2 vector: [lhmin, lhmax]
%
% Outputs:
%     Only graphics, in current Figure, unless a postscript file is created
%
% Assumes path can find personal functions:
%    bwsjpib.m
%    bwrswb.m
%    gpkde.m
%    gpnpr.m
%    gpsz1.m

%    Copyright (c) J. S. Marron 1998




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

if nargin <= 2 ;    %  at most 2 arguments input, use default range
  ivrange = 0 ;
else ;              %  vrange was specified, use that
  ivrange = vrange ;
end ;

if nargin <= 3 ;    %  at most 3 arguments input, use default nhp
  inhp = 11 ;
else ;              %  vrange was specified, use that
  if nhp == 0 ;
    inhp = 11 ;
  else ;
    inhp = nhp ;
  end ;
end ;

if nargin <= 4 ;    %  at most 4 arguments input, create no postscript file
  ipsfname = [] ;    %  Default
else ;
  ipsfname = psfname ;    %  Create Postscript file, with given prefix
end ;

if nargin <= 5 ;    %  at most 5 arguments input, so plot data
  iidatap = 1 ;    %  Default
else ;
  iidatap = idatap ;    %  Plot or not, depending on input value
end ;

if nargin <= 6 ;    %  at most 6 arguments input, so use default vhlp
  ivlhp = 0 ;    %  Default
else ;
  ivlhp = vlhp ;
end ;






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


disp('    Working on family') ;
subplot(3,1,1) ;
  if length(ivrange) == 1 ;
    mind = min(xdat) ;
    maxd = max(xdat) ;
    ngrid = 401 ;
  else ;
    mind = vrange(1) ;
    maxd = vrange(2) ;
    if length(ivrange) == 3 ;
      ngrid = vrange(3) ;
    else ;
      ngrid = 401 ;
    end ;
  end ;
  centd = mean([mind;maxd]) ;

  %  Set h grid stuff, as in SiZer1
  range = maxd - mind ;
  binw = range / (ngrid - 1) ;
  hmin = 2 * binw ;
  hmax = range ;
  vh = logspace(log10(hmin),log10(hmax),inhp) ;

  if idatyp == 1 ;    %  doing density estimation
    hcent = bwsjpib(data) ;
  else ;              %  doing regression
    hcent = bwrswb(data) ;
  end ;
  if hcent < min(vh) ;    %  if data based h below range
    disp('!!!   Warning from gpanal: databased h below this range   !!!') ;
    hcflag = 0 ;
  elseif hcent > max(vh) ;    %  if databased h above this range
    disp('!!!   Warning from gpanal: 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(vh))) ;
    hcind = inhp + 1 - hcind ;
                        %  since handles appear in reverse order in plot
    hcflag = 1 ; 
  end ;
  if idatyp == 1 ;    %  doing density estimation
    [msmoo, xgrid] = gpkde(data,vh,[mind; maxd; ngrid]) ;
    bottom = 0 ;
  else ;              %  doing regression
    [msmoo, xgrid] = gpnpr(data,vh,[mind; maxd; ngrid]) ;
    bottom = min(min(msmoo)) ;
  end ;
  top = max(max(msmoo)) ;
  range = top - bottom ;
  bottom = bottom - 0.05 * range ;
  top = top + 0.05 * range ;

           %  Next lines from gpfam.m
  plot(xgrid,msmoo,'c') ;
          %  Plot most curves in cyan
    vcurvh = get(gca,'Children') ;
          %  Vector of handles for curves
    if hcflag == 1 ;
      set(vcurvh(hcind),'LineWidth',2) ;
          %  use fatter line width for data based choice
      set(vcurvh(hcind),'Color',[1 0 0]) ;
          %  use red color for one in middle
    end ;

    vax = [mind,maxd,bottom,top] ;
    axis(vax) ;
    title(['Family Plot, from gpszj, ' date]) ;

    if iidatap == 1 ;    %  then overlay a plot of the raw data
      hold on ;
        if idatyp == 1 ;    %  doing density estimation
          yrand = vax(3) + (0.8 + 0.1 * rand(size(data))) ...
                                                 * (vax(4) - vax(3)) ;
               %  y coords for jittering
          plot(data,yrand,'g.') ;
        else ;              %  doing regression
          plot(data(:,1),data(:,2),'g.') ;
        end ;
      hold off ;
    end ;



disp('    Working on SiZer') ;
subplot(3,1,2) ;
  gpsz1(data,[mind,maxd,ngrid]) ;
    title(['SiZer 1 Plot, from gpszj, ' date]) ;
          %  next add marking lines
    hold on ;
      if hcflag == 1 ;
        plot([mind; maxd], ones(2,1)*log10(hcent),'-k') ;
      end ;
      plot(centd + 2*vh, log10(vh),':k') ;
      plot(centd - 2*vh, log10(vh),':k') ;
    hold off ;
    %  Note, one more bar gets added further down



if iicolor ~= 0 ;     %  Then go for nice colors in sizer and sicom
  %  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 ;




disp('    Working on linear bandwidths') ;
subplot(3,1,3) ;

  %  Here are lines from:  gpsz1
  %
  %  detect whether density or regression data
  %
  ialpha = 0.05 ;

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

  %  Set x grid stuff
  %

%     vrange  - 0 (or not specified) use the min and max of the data
%               2 vector:  use minx as vrange(1) and maxx as vrange(2)
%               3 vector:  use minx as vrange(1) and maxx as vrange(2)
%                               and number of grid points as vrange(3)

  n = length(xdat) ;
  if ivrange == 0 ;   %  then use standard default x grid
    ivxgp = [min(xdat),max(xdat),401] ;
  elseif length(ivrange) == 2 ;   %  then use standard default number of grid points
    ivxgp = [ivrange(1),ivrange(2),401] ;
  else ;     %  then use input values
    ivxgp = ivrange ;
  end ;
  left = ivxgp(1) ;
  right = ivxgp(2) ;
  ngrid = ivxgp(3) ;
  xgrid = linspace(left,right,ngrid)' ;
          %  row vector to evaluate smooths at

  %  Set h grid stuff
  %
  range = right - left ;
  binw = range / (ngrid - 1) ;
  if ivlhp == 0 ;    %  use default linear h grid
    hmin = 2 * binw ;
    hmax = range ;
    hmax = exp(mean([log(hmin); log(hmax)])) ;
          %  replace maximum bandwidth, by geometric mean
          %       (this way linear scale shows more interesting bandwidth)
  else ;     %  use input values
    hmin = ivlhp(1) ;
    hmax = ivlhp(2) ;
  end ;
  vlh = linspace(hmin,hmax,inhp) ;
          %  linearly spaced rid of bandwidths



  %  Bin the data
  %
  if idatyp == 1 ;        %  Treating as density estimation
    bincts = gplbinr(xdat,ivxgp) ;
  elseif idatyp == 2 ;    %  Treating as regression
    bincts = gplbinr([xdat ydat],ivxgp) ;
    bincts2 = gplbinr([xdat, ydat.^2],ivxgp) ;
  end ;
  n = sum(bincts(:,1)) ;
          %  put this here in case of truncations during binning



  %  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:inhp ;
    h = vlh(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)


    if idatyp == 1 ;       %  Doing density estimation
          %  main lines from gpkde.m, via nmsur5.m

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

      vd = conv(bincts,kvecd) ;
      vv = conv(bincts,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 ;
          %  did this outsidea loop in nmsur5.m

    else ;    %    Doing regression
          %  using modification of lines from gpnpr.m

        %  main lines from gpnpr.m, via szeg4.m

      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) .* xgrid , kvec) ;
      s1 = s1(k+1:k+ngrid) ;
      s2 = conv(bincts(:,1) .* xgrid.^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) .* xgrid , kvec) ;
      xbar = xbar(k+1:k+ngrid) ;
          %  Weighted sum of X_i
      xbar = xbar ./ ve ;
          %  Weighted avg of X_i

      t1 = conv(bincts(:,2) .* xgrid , 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)

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

      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 .* xgrid , kvec.^2) ;
      u1 = u1(k+1:k+ngrid) ;
      u2 = conv(bincts(:,1) .* sig2res .* xgrid.^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)

    end ;


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


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

  end ;


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



  image([left,right],[hmin,hmax],mapout) ;
    set(gca,'YDir','normal') ;
%    colormap(cocomap) ;  
    xlabel('x') ;
    ylabel('h') ;
    title(['Linear Bandwidth Scale, from gpszj, ' date]) ;
          %  next add marking lines
    hold on ;
      if hcflag == 1 ;
        plot([mind; maxd], ones(2,1)*hcent,'-k') ;
      end ;
      plot(centd + 2*vlh, vlh,':k') ;
      plot(centd - 2*vlh, vlh,':k') ;
    hold off ;




%  Now go back and add bar to 2nd plot, showing top of first
%
subplot(3,1,2) ;
  hold on ;
    plot([mind; maxd], ones(2,1)*log10(hmax),'--k') ;
  hold off ;




if length(ipsfname) > 0 ;     %  then create postscript file with given name
    orient tall ;

  if iicolor ~= 0 ;     %  then make color postscript
    eval(['print -dpsc \matlab\steve\ps\' ipsfname '.ps']) ;
  else ;                %  then make black and white
    eval(['print -dps \matlab\steve\ps\' ipsfname '.ps']) ;
  end ;

end ;

