function tk = patjitRNGint(tj,rk,alphak,nk,betak)
% function tk = patjitRNGint(tj,rk,alphak,nk,betak)
% function tk = patjitRNGint(tj,struct containing r,alpha,n,beta)
%
% used by patjitSample.m
%
% samples the next spike given the previous spike time using the
% pattern jitter sampling method

% 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;
    nk = rk.n;
    betak = rk.beta;
    rk = rk.r;
end

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

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

% get the earliest next spike time
% map to indices
tj = tj+(rk+2-alphak);

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

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

% loop over spikes
alphak1 = alphak-1;
for h = 1:numel(tj)

    tjh = tj(h);
    if isnan(tjh), tk(h) = NaN; continue, end
    
    % rescale to account for previous spike
    if tjh > 1
        tkh = tk(h)*betak(tjh);
    else
        tkh = tk(h);
    end
    
    % figure out which segment
    for i = nk:-1:1
        if tkh <= betak(i)
            break
        end
    end
    if i < tjh
        [h tjh i]
        error('algorithm error 1')
    end
    
    % store the sample
    tk(h) = i + alphak1;
end
