%
%
% Inhomogeneous Poisson process
%
% Fabrizio Gabbiani 11/10/08
%

n_rep = 10;

%isi data
isis_pp = [];

%constant rate PP
figure(1);

for i = 1:n_rep
  %loop over the n_rep iterations
  info_str = sprintf('iteration %i...',i);
  disp(info_str);
  
  %generate the poisson spike train
  [t_vect, spk_vect] = inhom_pp(40);
  subplot(11,1,i);
  plot(t_vect,spk_vect);
  set(gca,'XTick',[],'YTick',[]);
  
  %compute isis and add to the isis_pp vector
  spk_inds = find(spk_vect > 0.5);
  spk_times = t_vect(spk_inds);
  isis = spk_times(2:end) - spk_times(1:end-1);
  isis_pp = [isis_pp isis];
end;

subplot(11,1,11);
plot([0 1000],[40 40]);
xlabel('time (ms)');
ylabel('rate');
set(gca,'YTick',[]);

%isi data
isis_inpp = [];

%sinusoidal rate PP
figure(2);
for i = 1:n_rep
  %loop over the n_rep iterations
  info_str = sprintf('iteration %i...',i);
  disp(info_str);
  
  %generate the inhomogenous PP
  [t_vect, spk_vect, rho_vect] = inhom_pp2(40,40,10);
  subplot(11,1,i);
  plot(t_vect,spk_vect);
  set(gca,'XTick',[],'YTick',[]);
  
  %compute isis and add to the isis_pp vector
  spk_inds = find(spk_vect > 0.5);
  spk_times = t_vect(spk_inds);
  isis = spk_times(2:end) - spk_times(1:end-1);
  isis_inpp = [isis_inpp isis];
end;

subplot(11,1,11);
plot(t_vect,rho_vect);
xlabel('time (ms)');
ylabel('rate');
set(gca,'YTick',[]);

%plot the isi data
figure(3)

xout = 0:2:250;
n_inpp = hist(isis_inpp,xout);
n_inpp = n_inpp/sum(n_inpp);

n_pp = hist(isis_pp,xout);
n_pp = n_pp/sum(n_pp);

bar(xout,[n_pp; n_inpp]');
xlabel('isi (ms)');
ylabel('probability');