% 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)
rand("seed",1234) % octave
		     
path(path,"./octave_functions");% where we will find some
% additional octave repository functions
% that do not come with the main package. Alternatively, you download them from
% http://octave.sourceforge.net/  (and install them properly on your system)


% -------------------------------------------------------------------------
% ------------- 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);
xlabel('time in seconds');
ylabel('spike count');
title('spike train histogramm');
bar(t, hist(tf,t))

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

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

% - 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);
xlabel('time in seconds');
ylabel('spike count');
title('spike train histogramm');
bar(t, hist(tf_r,t))

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



% ------------------------------------------------------------------
% ------------ exercise 14 -- Inhomogenous Poissonprocess ----------

 
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);
xlabel('Firing times of a homogeneous Poisson process')
ylabel('Firing times of an inhomogeneous Poisson process')
title('')
plot(t_hom,t_inh,'--o;;',[0 1],[0 1],'--;;')



% 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);

clear isi;

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


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

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

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