% Computational Neuroscience II: Foundations of Neural Coding, SS 2003
% Solutions for Assignment 5 (22.05.03)
% (c) written by Christian Leibold and Richard Kempter

clear; % --- forget
rand('state',1234) % initialize random number generator at state 1234
                   % (omit the last line if you want to get differet
                   % random numbers in each run)

% -------------------------------------------------------------------------
% ------------- exercise 13-- Random Distributions ------------------------

% Algorithm X = -\lambda_0^{-1} log(1-u)) implemented in function Poisson

% --------------  10 second Poisson train

n=1;         % Spikecounter
lambda0=10.; % rate in units of Hz

% - initialize firing time vector tf
tf(n) = Poisson(lambda0); % first spike at the end of the first ISI 
while tf(n)<10,
  isi(n) = Poisson(lambda0);	    
  n = n+1;
  tf(n) = tf(n-1) + isi(n-1);
end

% - define grid:
DT=0.05;
t=0:DT:11;

figure(1);
subplot(2,1,1);
bar(t, hist(tf,t))
xlabel('time in seconds');
ylabel('spike count');
title('spike train histogramm');

% - define isigrid
DTISI=0.002;
tisi=0:DTISI:.1;

subplot(2,1,2);
bar(tisi, hist(isi,tisi))
xlabel('time in seconds');
ylabel('ISI count');
title('ISI distribution');
% Interpret the peak at 0.1 ms in the ISI distribution (hint: >> help hist)!
% What is the difference between 'hist' and 'histc' ?

% - refractoriness
tref = 0.02; 
n=1;
tf_r(n) = Poisson(lambda0);
while tf_r(n)<10,
  isi_r(n) = tref + Poisson(lambda0);	    
  n = n+1;
  tf_r(n) = tf_r(n-1) + isi_r(n-1);
end


% - define grid:
DT=0.05;
t=0:DT:11;

figure(2);
subplot(2,1,1);
bar(t, hist(tf_r,t))
xlabel('time in seconds');
ylabel('spike count');
title('spike train histogramm');

%define isigrid
DTISI=0.002;
tisi=0:DTISI:.1;
subplot(2,1,2);
bar(tisi, histc(isi_r,tisi)); % try also 'hist' here
xlabel('time in seconds');
ylabel('ISI count');
title('ISI distribution');



% ------------------------------------------------------------------
% ------------ exercise 14 -- Inhomogeneous Poisson process ----------

 
T =1;  % period in units of seconds

% homogeneous Poisson train
n=1;         % Spikecounter
lambda0=100.; % rate in units of Hz
t_hom(n) = Poisson(lambda0); % first spike at the end of the first ISI 
while t_hom(n)<T,
  isi_hom(n) = Poisson(lambda0);	    
  n = n+1;
  t_hom(n) = t_hom(n-1) + isi_hom(n-1);
end

% inhomogeneous Poisson train
t_inh = timewarp(t_hom);

figure(3)
plot(t_hom,t_inh,'--o',[0 1],[0 1],'--')
xlabel('Firing times of a homogeneous Poisson process')
ylabel('Firing times of an inhomogeneous Poisson process')



% 20 more 1-second trains of an inhomogeneous Poisson process :
N = 20;
r=10;  % rate in units of Hz 
f=3.3; % frequency 3.3 Hz-- see lambda_t.m 

DT=0.002;  % Integration stepsize: 2 ms
T=1; % Duration: 1 second

%Grids
t=0:DT:T;
tcorr=-T:DT:T;

%stimulus for spike stimulus correlation ssavg
stim=lambda_t(r,f,t);

%Pad averages with zeros
ravg=zeros(1,length(t));
xavg=zeros(1,2*length(t)-1);
ssavg=zeros(1,2*length(t)-1);

for i=1:N,
	[spiketrain,isi]=inhPoisson(DT,r,f,T);
        spikehist=histc(spiketrain,t);
        
        ssavg = ssavg + xcorr(spikehist,stim,'unbiased');
        xavg = xavg + xcorr(spikehist,spikehist,'unbiased');
        ravg = ravg + spikehist;
end


figure(4);
subplot(3,1,1);
bar(t,ravg/DT/N);
xlabel('time in seconds');
ylabel('rate in Hz');
title('average rate');

subplot(3,1,2);
bar(tcorr,xavg/N/(length(xavg)*DT)^2);
xlabel('time in seconds');
ylabel('correlation in Hz^2');
title('average spike autocorrelation');

subplot(3,1,3);
bar(tcorr,ssavg/N/(length(xavg)*DT));
xlabel('time in seconds');
ylabel('correlation in Hz^2');
title('average spike-stimulus correlation');