% HhatCE
% 11/2/18
function [pout] = HhatCE(z, x, p)
% Computes HhatCE probs pout from distribution Hhat defined by <x, p>, at z.  
%  CE = continuous and linearly extrapolated.
%  pout = 0 if extrapolated to negative value.
%  pout = 1 if extrapolated to value > 1.
% x, p ~ 1 X Nin.  p > 0.  Must sum to 1, be sorted by x.
% z, pout ~ 1 X Nout.  Nout may differ from Nin
% z must be increasing. 
% x and p may be pre-sorted by
%   [x, I] = sort(x);
%   p = p(I);
%  or 
%   xp = sortrows([x' p'])';
%   x = xp(1,:);
%   p = xp(2,:);
%  or 
%   [x; p] = sortrows([x' p']');
Nin = size(x, 2);  % must equal size(p, 2).
Nout = size(z, 2);
P = cumsum(p);     % 1 X Nin  P(0) implicitly 0, P(Nin) = 1
Q = P - p/2;
J = nan(1,Nout);
% J(i) = 2 if z(i) < x(2)
% J(i) = Nin if z(i) >= x(Nin)
% else J(i) = smallest j such that z(i) < x(j) 
j0 = 2;
for i = 1:Nout
  if z(i) >= x(Nin)
    J(i) = Nin;
    continue % to next i
  end % if z(i)
  for j = j0:Nin
    if z(i) < x(j)
      J(i) = j;
      j0 = j;
      break % out of j loop
    end % if q(i)
  end % for j
end % for i

theta = (z - x(J - 1)) ./ (x(J) - x(J-1));
pout = theta .* Q(J) + (1-theta) .* Q(J - 1);
pout = max(pout, 0);
pout = min(pout, 1);

return


  
  