% DrawCEsample
% 5/17/18
function [xout, JS] = DrawCEsample(q, x, p)
% Draws CE sample xout from distribution defined by <x, p>, at prob. q.  
%  CE = continuous and linearly extrapolated.
% x, p ~ 1 X Nxp.  p > 0.  Must sum to 1, be sorted by x.
% q, xout ~ 1 X N
% q must be increasing, in (0, 1). 
% q need not be evenly spaced, but then must be derived from pxout > 0 
%  as q = Pxout - pxout/2, where Pxout = cumsum(pxout).
%  Calling program must then supply xout weights equal to pxout.
% JS gives Step Function indices so x(JS) is Step Fn draw.
%   If q reflects equal weights, DrawStratSample.m isfaster.
% JS(i) = smallest j such that q(i) <= P(j)
% 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']');
N = size(q, 2);
Nxp = size(x, 2);  % must equal size(p, 2).
P = cumsum(p);     % 1 X Nxp
Q = P - p/2;
J = nan(1,N);
% J(i) = 2 if q(i) < Q(2)
% J(i) = Nxp if q(i) >= Q(Nxp)
% else J(i) = smallest j such that q(i) < Q(j) 
j0 = 2;
for i = 1:N
  if q(i) >= Q(Nxp)
    J(i) = Nxp;
    continue % to next i
  end % if q(i)
  for j = j0:Nxp
    if q(i) < Q(j)
      J(i) = j;
      j0 = j;
      % Inactive bugtrap:
      if J(i) < 2
        iJ = i
        qi = q(i)
        Ji = J(i)
      end
      break % out of j loop
    end % if q(i)
  end % for j
  % Inactive bugtrap:
  if isnan(J(i))
    iJnan = i
    Jinan = J(i)
    qi = q(i)
  end
    
end % for i
% Inactive bugtrap:
Jnan = find(isnan(J));
sumJnan = sum(Jnan);
if sumJnan > 0
  sumJnan = sumJnan
  sumJnan = min([sumJnan; 10])
  Jnan = Jnan(1:sumJnan)
  Qnan = Q(Jnan(1:sumJnan))
end
theta = (q - Q(J - 1)) ./ (Q(J) - Q(J - 1));
xout = theta .* x(J) + (1-theta) .* x(J - 1);
JS = J - (q < P(J-1));
return


  
  