% generate_spike_trains.m
%
% Anthony Kellems, 6/11/08
% modifications done by Katie Ward, 9/22/08
%
% Generates N Poisson spike trains with rate lambda Hz for Tmax seconds
% from a homogeneous Poisson process.
%
% train = generate_spike_trains_vec(N,lambda,Tmax)
% N: number of neurons
% lambda: vector of spike rates of size Nx1 (spikes/sec)
% Tmax: maximum time (sec)
%
% trains: Matrix of Poisson spike trains
%   Col 1: times of spikes (sec)
%   Col 2: index of neuron spiking
%
% Ex: train = generate_spike_trains_vec(10,[1:10]',20)
%

function trains = generate_spike_trains_vec(N,lambda,Tmax)

[rows,cols] = size(lambda);
if rows~=N || cols~=1
    disp('Error in generage_spike_trains_vec: lambda must be Nx1')
    return
end

lambdaMax = max(lambda);

appx = ceil(1.003*Tmax*lambdaMax*N);
trains = zeros(appx,2);     % preallocate approximate # of spikes
                                     % The 1.01 is a factor to overestimate
                                     % spikes a bit, which gives a big
                                     % speedup in case we go over the
                                     % estimate Tmax*lambda*N. Change this
                                     % factor as you see fit.
                                                                          
ct = 0;
ind = [1:N]';

val = -log(rand(N,1)) ./ lambda;
trains(ct+ind,1) = val;      % first spike
trains(ct+ind,2) = ind;

ct = ct + N;

ind = find(trains(1:N,1) < Tmax);
neur = ind;
num = length(neur);

old_times = val(ind);

while num ~= 0  %stop when all neurons have a spike time > Tmax 
    
    new = -log(rand(num,1)) ./ lambda(neur);
    
    val = old_times + new;
    trains(ct+[1:num]',1) = val;     % Poisson process: generate spikes
    trains(ct+[1:num]',2) = neur;
    
    ct = ct + num;
    
    ind = find(val < Tmax);
    neur = neur(ind);
    num = length(neur);       

    old_times = val(ind);
    
end

trains(ct+1:end,:) = [];                % remove unused preallocated space

%Sort trains
[trains(:,1),IX] = sort(trains(:,1));
trains(:,2) = trains(IX,2);
                                    