function tk = patjitRNG(tj,rk,alphak,ck,xk,dk,mk,nk,betak)
% function tk = patjitRNG(tj,rk,alphak,ck,xk,dk,mk,nk,betak)
% function tk = patjitRNG(tj,struct containing r,alpha,c,x,d,m,n,beta)
%
% used by patjitSample.m
%
% samples the next spike given the previous spike time using the
% pattern jitter sampling method
%
% If you encounter the following error, then the history dependence has
% gotten too complicated and rounding errors accumulated.  Even minor
% changes to the pattern size or the jitter length may fix this problem.
% Reverting to discrete jitter will always fix the problem.
%
% ???  In an assignment  A(I) = B, the number of elements in B and I must be the same.
%
% Error in ==> patjitRNG at 148
%    tk(h) = tkh;
    
% Copyright (C) 2008 Matthew T Harrison
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    See <http://www.gnu.org/licenses/> for a copy of the license.

% check for struct input
if isstruct(rk)
    % copy inputs
    alphak = rk.alpha;
    ck = rk.c;
    xk = rk.x;
    dk = rk.d;
    mk = rk.m;
    nk = rk.n;
    betak = rk.beta;
    rk = rk.r;
end

% check for a priori impossible
if isnan(alphak)
    tk = nan(size(tj));
    return
end

alphak1 = alphak(1);

% check for fixed spike time
if nk == 0
    tk = tj + rk;
    tk(tk >= alphak1) = NaN;
    tk(~isnan(tk)) = alphak1;
    return
end

% get the earliest next spike time (not inclusive)
tj = tj+rk;

% deal with impossible given this spike time
tj(tj >= alphak(end)) = NaN;

% store the uniform random numbers in the output
tk = rand(size(tj));

% check for uniform on a single interval (easy and frequent)
if nk == 1 && mk == 2
    % sample on this interval
    % tj = max(tj,alphak(1)); % this doesn't work because max ignores nan
    tj(tj < alphak1) = alphak1;
    tk = tk.*(alphak(2)-tj)+tj;
    return
end

% loop over spikes
for h = 1:numel(tj)

    tjh = tj(h);
    if isnan(tjh), tk(h) = NaN; continue, end
    tkh = tk(h);

    % rescale to account for the conditioning
    imin = 1;
    if tjh > alphak1
        for imin = 2:nk+1
            if tjh < alphak(imin)
                imin = imin - 1;
                break % always break here
            end
            % exit loop here should have been caught by impossible above
        end

        % copy the code for mypolyval to avoid overhead

        %tkh = tkh * mypolyval(ck(:,imin),1-(tjh-xk(imin))/dk(imin),mk(imin));

        %%%%%%%%%%%%%%%%%%%
        % start mypolyval
        %%%%%%%%%%%%%%%%%%%

        c = ck(:,imin);
        z = 1-(tjh-xk(imin))/dk(imin);
        m = mk(imin);
        p = c(1);

        for j = 2:m
            p = p.*z + c(j);
        end

        tkh = tkh * p;

        %%%%%%%%%%%%%%%%%%%
        % end mypolyval
        %%%%%%%%%%%%%%%%%%%

    end

    % figure out which segment
    for i = nk:-1:1
        if tkh <= betak(i)
            break
        end
    end
%     if i < imin
%         [h imin i]
%         error('algorithm error 1')
%     end

    % copy the relevant quantities
    m = mk(i);        
    a = max(tjh,alphak(i));
    b = alphak(i+1);
    d = dk(i);
    x = xk(i);
    % find the zero
    % deal with linear and quadratic separately (easy and frequent)
    if m == 2
        tkh = d*(1+(ck(2,i)-tkh)/ck(1,i))+x;
    elseif m == 3
        B = ck(2,i)/ck(1,i);
        D = sqrt(B*B-4*(ck(3,i)-tkh)/ck(1,i));
        tkh = d*(1-(D-B)/2)+x;
    else
        c = ck(1:m,i);
        c(m) = c(m)-tkh;
        tkh = d*(1-polyroot(c,1-(b-x)/d,1-(a-x)/d))+x;
    end
    
%     if isempty(tkh) || tkh < a || tkh > b
%         [h m a tkh b]
%         error('algorithm error 2')
%     end

    % store the sample
    tk(h) = tkh;
end
