function t = patjitSample(N,SS,t0)
%function t = patjitSample(N,SamplerStructure,t0)
%
% draws N samples from the description in SamplerStructure 
% t0 (default: -inf) are the initial spike times
%
% SamplerStructure can be created with PatternJitter.m
%
% uses: 
% patjitRNG.m (mypolyval.m, polyroot.m)
% patjitRNGint.m

% 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.

%%%%%%%%%%%%%%%%%%
% input checking %
%%%%%%%%%%%%%%%%%%

if N <= 0 || isempty(SS)
    t = [];
    return
end

if nargin < 3 || isempty(t0), t0 = -inf; end
if numel(t0) == 1
    t0 = repmat(t0,1,N);
elseif numel(t0) ~= N
    error('t0 incorrectly sized')
else
    t0 = t0(:).';
end

%%%%%%%%%%%%%%%%%%%%%
% precomputations
%%%%%%%%%%%%%%%%%%%%%

% number of patterns
m = numel(SS);

% extract some useful info
patn = zeros(m,1);
for k = 1:m
    patn(k) = SS(k).patn;
end

% where the patterns start and end
ndx = [0 ; cumsum(patn)];

% total number of spikes
nt = ndx(end);

if nt == 0, t = []; return, end

% initialize
t = zeros(nt,N);

%%%%%%%%%%%%%%%%%%%%
% sampling
%%%%%%%%%%%%%%%%%%%%

% loop over spikes
for k = 1:m
    if SS(k).Greal
        t0 = patjitRNG(t0,SS(k));
    else
        t0 = patjitRNGint(t0,SS(k));
    end
    % add the pattern
    pat = SS(k).pat;
    for j = 1:patn(k)
        t(ndx(k)+j,:) = t0+pat(j);
    end
end
