% ONE1
%   Exercise 1 from CHAPTER 1 of Dayan and Abbott

%   NEUR 1680
%   Last updated 02/08/2012

% Generate spikes for 10 s (or longer if you want better statistics) using
% a Poisson spike generator with a constant rate of 100 Hz, and record
% their times of occurrence.

% Compute the coefficient of variation of the
% interspike intervals, and the Fano factor for spike counts obtained
% over counting intervals ranging from 1 to 100 ms.

% Plot the interspike interval histogram.

%%%%%%%%%%%%%%%%%%


clear
close all
set(0,'DefaultFigureWindowStyle','docked')

%%%%%%  PARAMETERS   

duration = 10;                    % duration of recording [sec]
dt = .0001;                        % discrete time step    [sec]  (this might be overkill)
rate = 100;                        % firing rate           [Hz]

%%%%%%%%%%%%%%%%%%

N = round(duration/dt);             % number of time steps in recording

% Probability of observing a spike at any given discrete time:
p = rate*dt;                        % [dimensionless]


% The Poisson spike train is generated by N independent Bernoulli random variables,
%     i.e., independent "coin flips":
% spike_train = rand(1,N) < p;  

% Alternative mechanism: renewal process with exponential ISI
spike_train = zeros(1,N);  
current_time = dt;
while current_time < duration
    spike_train((round(current_time/dt))) = 1;
    ISI = exprnd(1/rate);
    current_time = current_time + ISI;
end

number_spikes = sum(spike_train);   % total number of spikes in recording

figure(1)
subplot(2,1,1)
N_display = 10000;                   % number of time steps to display

%  A convenient way to display a spike train is to use a "stem plot": 
h = stem(dt*(1:N_display),spike_train(1:N_display));
set(h(1),'Marker','none')           % removes the circles from the stem plot
set(h(1),'LineWidth',.1)
set(gca,'YLim',[-.1 1.5])
title('HOMOGENEOUS  POISSON  SPIKE  TRAIN')
xlabel('time (sec)')
disp(['duration of recording: ' num2str(duration) ' sec'])
disp(['total number of spikes: ' num2str(number_spikes)])

%  Since we will compute the coefficient of variation in other problems as well,
%  it is convenient to use a "function" (see coef_var.m): 
[c_v,ISI] = coef_var(spike_train);  % ISI's and coefficient of variation of ISI
ISI = ISI*dt;                       % ISI's measured in seconds

% choose a bin length for the ISI histogram:
bin_l = 10*dt;                     % should be at least equal to dt

max_interval = max(ISI);
bin_edges = 0:bin_l:max_interval;    
ISI_histogram = histc(ISI,bin_edges);                           % ISI histogram
ISI_distribution = ISI_histogram/(bin_l*sum(ISI_histogram));    % sample distribution = normalized histogram
max_value = max(ISI_distribution);                              % largest value

subplot(2,1,2)
bar(bin_edges + bin_l/2, ISI_distribution, 1)                   % displays the sample distribution
set(gca,'XLim',[0 max_interval])
set(gca,'YLim',[0 max_value])
hold on
time_steps = 0:dt:max_interval;
plot(time_steps, rate*exp(-rate*time_steps), 'r', 'linewidth',2) % displays the theoretical distribution
hold off
title('sample distribution of ISI''s')
xlabel('interval (sec)')
disp(' ')
disp(['Coefficient of variation of ISI interval: ' num2str(c_v)])

% These will be the counting-interval lengths for Fano factors:
counting_intervals = [1 10 100 1000 10000 100000 1000000];      % units of dt

for c_i = counting_intervals
    [fano, n_i] = fano_factor(spike_train,c_i);     
    disp(['Fano factor over ' num2str(n_i) ' counting intervals of length ' num2str(c_i*dt) ' sec: ' num2str(fano)])
end
