% ONE2   Exercise 2 from CHAPTER 1 of Dayan and Abbott

%   NEUR 1680
%   Last updated 02/08/2012

% Add a refractory period to the Poisson spike generator by allowing
% the firing rate to depend on time. Initially, set the firing rate to a
% constant value r_0. After every spike, set the rate to 0 and then
% allow it to recover exponentially back to r_0 with a time constant tau_ref
% that controls the refractory recovery rate.

% Plot the coefficient of variation
% as a function of tau_ref over the range 1 ms <= tau_ref <= 20 ms 
% and plot interspike interval histograms for a few different values of
% tau_ref in this range. Compute the Fano factor for spike counts
% obtained over counting intervals ranging from 1 to 100 ms for the
% case tau_ref = 10 ms.


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

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

%%%%%%  PARAMETERS

duration = 100;                     % duration of recording [sec]
dt = .0001;                         % discrete time step    [sec]
r_0 = 100;                          % recovery firing rate [Hz]
tau_ref = [.0002 .0005 .001 .005 .010 .015 .020 .050 .1];  % time constant of recovery from refractoriness [sec]
                        % NOTE: this is more than we are asked.
                        %   Some of these parameter values are biologically unrealistic
                        %   but might be theoretically interesting

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


N = round(duration/dt);             % number of time steps in recording
spike_train = zeros(1,N);
rates       = zeros(1,N);
n_tau_ref = length(tau_ref);
C_V = zeros(1,n_tau_ref);

for t_r = 1:n_tau_ref
    
    rate = r_0;
    tau = tau_ref(t_r);
    
    for n = 1:N
        p = rate*dt;            % probability of observing a spike in bin n
        spike = rand < p;
        spike_train(n) = spike;
        rates(n) = rate;
        % depending on presence or absence of spike, rate resets to 0 or recovers exponentially:
        rate = (1 - spike) * (rate + (dt/tau)*(r_0 - rate)); 
    end
    
    number_spikes = sum(spike_train);
    
    figure(t_r)
    subplot(2,1,1)
    N_display = 1000;                   % number of time steps to display
    plot(dt*(1:N_display),rates(1:N_display)/r_0,'r')   % plots the rate
    hold on
    h = stem(dt*(1:N_display),spike_train(1:N_display));
    hold off
    set(h(1),'Marker','none')           % removes the circles from the stem plot
    set(h(1),'LineWidth',.1)
    set(gca,'YLim',[-.1 1.5])
    title(['time constant of recovery from refractoriness: \tau_{ref} = ' num2str(tau_ref(t_r)) ' sec'])
    xlabel('time (sec)')
    disp(' ')
    disp(['time constant of recovery from refractoriness: ' num2str(tau) ' sec'])
    disp(['total duration of recording: ' num2str(duration) ' sec'])
    disp(['total number of spikes: ' num2str(number_spikes)])
    
    [c_v,ISI] = coef_var(spike_train);  % ISI's and coefficient of variation of ISI
    ISI = ISI*dt;                       % ISI's measured in seconds
    C_V(t_r) = c_v;

    % choose a bin length for the histogram:
    bin_l = 1*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, r_0*exp(-r_0*time_steps), 'r', 'linewidth',2) % displays the theoretical distribution without refractoriness
    hold off
    title('sample distribution of ISI''s')
    xlabel('interval (sec)')
    disp(['Coefficient of variation of ISI interval: ' num2str(c_v)])
    disp(['Predicted Fano Factor: ' num2str(c_v^2)])
    % NOTE: WE ARE DEALING WITH A RENEWAL PROCESS HERE
    %       HENCE: FANO FACTOR = (CV)^2 for large enough counting interval -- see book p.27
    
    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
    
end

figure(t_r+1)
plot(tau_ref,C_V,'x-')
title('coefficient of variation of ISI distribution as a function of \tau_{ref}')
grid
xlabel('\tau_{ref} (sec)')
ylabel('coefficient of variation')
