function gpcheta(mcurv,mevents,itest,mwts) 
% GPCHETA, General Purpose Curve Heights at Event Times Analysis
%     This compares interaction between heights of a given curve,
%     and a set of "event times", by evaluating the curve at
%     the event times, and then comparing those numbers to 
%     similar values for some corresponding "background event
%     times" (e.g. randomly chosen times), to see if there is a
%     connection between the curve and the events.
% Inputs:
%     mcurv    -  2 x ng matrix defining the curve, 
%                     1st col is "time" values
%                     2nd col is "curve heights"
%     mevents  -  2 x n matrix of event times
%                     1st col is times of events of interest
%                     2nd col is times of "background events"
%     itest    -  scalar (optional) parameter controlling type of test
%                     1 - (default) test for event mean << background mean
%                     2 - test for event mean >> background mean
%                     3 - (not implemented yet) event spread << background spread
%                     4 - (not implemented yet) event spread >> background spread
%     mwts     -  either 0 for equal weights (default), or
%                 2 x n matrix (optional) of weights to use on event times
%                     1st col for events of interest
%                     2nd col for background events
%                     entries should be >= 0 , and should sum to 1
%                     along each column.
% Output:
%     graphics only, in a new figure window

%    Copyright (c) J. S. Marron 1999



%  Set input parameters and defaults according to number of input arguments
%
if nargin == 2 ;    %  only 2 arguments input, so use default
  iitest = 1 ;
else ;              %  then use input value
  iitest = itest ;
end ;

if nargin <= 3 ;    %  at most 3 arguments input, so use default
  imwts = 0 ;
else ;              %  then use input value
  imwts = mwts ;
end ;



%  set internal parameters
%
ng = size(mcurv,1) ;
          %  number of curve grid points
n = size(mevents,1) ;
          %  number of events in time
left = min(mcurv(:,1)) ;
right = max(mcurv(:,1)) ;
bottom = min(mcurv(:,2)) ;
top = max(mcurv(:,2)) ;
crange = top - bottom ;
bottom = bottom - 0.1 * crange ;
top = top + 0.1 * crange ;



%  Upper plot shows curve, and jitter plots of events
%
figure ;

subplot(2,1,1) ;
  plot(mcurv(:,1),mcurv(:,2),'-w') ;
    title('Curve Heights at Event Times Analysis') ;
    xlabel('time') ;
    ylabel('curve') ;
    axis([left, right, bottom, top]) ;

    hold on ;
      seed = 46586323 ;
      rand('seed',seed) ;

      randht = bottom + (.5 + .1 * rand(n,1)) * (top - bottom) ;
      plot(mevents(:,1),randht,'.r') ;
        tx = left + .01 * (right - left) ;
        ty = bottom + .55 * (top - bottom) ;
      text(tx,ty,'Event Times') ;

      randht = bottom + (.2 + .1 * rand(n,1)) * (top - bottom) ;
      plot(mevents(:,2),randht,'.b') ;
        ty = bottom + .25 * (top - bottom) ;
      text(tx,ty,'Background Times') ;
    hold off ;



%  Next evaluate curve at event times
%
flag = (mevents(:,1) < left) ;
nlo = sum(flag) ;
if nlo > 0 ;    %  then there are some event times below left end,
                      %  so truncate and send warning message
  evt = mevents(~flag,1) ;

  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
  disp('!!!    Warning from gpcheta:                    !!!') ;
  disp(['!!!       had ' num2str(nlo) ' main events below left end   !!!']) ;
  disp('!!!       will truncate them and proceed        !!!') ;
  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
else ;
  evt = mevents(:,1) ;
end ;

flag = (mevents(:,2) < left) ;
nlo = sum(flag) ;
if nlo > 0 ;    %  then there are some background times below left end,
                      %  so truncate and send warning message
  bgt = mevents(~flag,2) ;

  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
  disp('!!!    Warning from gpcheta:                    !!!') ;
  disp(['!!!       had ' num2str(nlo) ' bkgd events below left end   !!!']) ;
  disp('!!!       will truncate them and proceed        !!!') ;
  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
else ;
  bgt = mevents(:,2) ;
end ;

flag = (evt > right) ;
nhi = sum(flag) ;
if nhi > 0 ;    %  then there are some event times above right end,
                      %  so truncate and send warning message
  evt = evt(~flag) ;

  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
  disp('!!!    Warning from gpcheta:                    !!!') ;
  disp(['!!!       had ' num2str(nhi) ' main events above right end   !!!']) ;
  disp('!!!       will truncate them and proceed        !!!') ;
  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
end ;

flag = (bgt > right) ;
nhi = sum(flag) ;
if nhi > 0 ;    %  then there are some background times above right end,
                      %  so truncate and send warning message
  bgt = bgt(~flag) ;

  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
  disp('!!!    Warning from gpcheta:                    !!!') ;
  disp(['!!!       had ' num2str(nhi) ' bkgd events above right end   !!!']) ;
  disp('!!!       will truncate them and proceed        !!!') ;
  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
end ;

nev = length(evt) ;
nbg = length(bgt) ;


if imwts == 0 ;   %  then do equally weighted stuff (default)

  cevt = interp1(mcurv(:,1),mcurv(:,2),evt) ;
  cbgt = interp1(mcurv(:,1),mcurv(:,2),bgt) ;

  if (iitest == 1) | (iitest == 2) ;    %  then do test based on means
    evavg = mean(cevt) ;
    bgavg = mean(cbgt) ;
    sig = sqrt(std(cevt)^2 / nev + std(cbgt)^2 / nbg) ;
          %  estimated standard deviation for 2-sample test
    stat1str = ['Event Avg. = ' num2str(evavg)] ;
    stat2str = ['Bkgd Avg. = ' num2str(bgavg)] ;
    if iitest == 1 ;    %  then test for event mean small
      pval = 1 - phi((bgavg - evavg) / sig) ;
    elseif iitest == 2 ;    %  then test for event mean large
      pval = phi((bgavg - evavg) / sig) ;
    end ;
    pvalstr = ['p-val for diff. = ' num2str(pval)] ;
  elseif (iitest == 3) | (iitest == 4) ;    %  then do test based on spreads
    evstd = std(cevt) ;
    bgstd = std(cbgt) ;
    stat1str = ['Event s.d. = ' num2str(evstd)] ;
    stat2str = ['Bkgd s.d. = ' num2str(bgstd)] ;
  
    disp('!!! warning from gpcheta:  test not implemented yet !!!') ;
    disp('!!!     so returning p-value of 1/2                 !!!') ;
    pvalstr = 'p-val for diff. = 1/2' ;
  end ;

  [evkde,cgrid] = gpkde(cevt,0,[bottom,top]) ;
  bgkde = gpkde(cbgt,0,[bottom,top]) ;

else ;    %  do weighted version of everything
  disp('!!!  Warning from gpcheta:  weighted version not yet implemented !!!') ;

end ;



%  Lower plot analyzes distributions of curve heights
%
subplot(2,1,2) ;
  plot(cgrid,evkde,'-r',cgrid,bgkde,'--b') ;
    title(['CHETA, ' date]) ;
      maxkde = max([evkde; bgkde]) ;
    axis([bottom, top, 0, 1.1*maxkde]) ;
    xlabel('curve') ;
    ylabel('density') ;
    legend('Event Times','Background',0) ;
          %  0 for "automatic placement
      tx = bottom + .05 * (top - bottom) ;
      ty = 0 + .9 * (1.1*maxkde - 0) ;
    text(tx,ty,stat1str) ;
      ty = 0 + .8 * (1.1*maxkde - 0) ;
    text(tx,ty,stat2str) ;
      ty = 0 + .7 * (1.1*maxkde - 0) ;
    text(tx,ty,pvalstr) ;



