%% load data (run this cell from the crcns_ret-1 directory)
datafile = '20080516_R1.mat'; % data file name
load([pwd filesep 'Data' filesep datafile]); % load data

load([pwd filesep 'stimulus.mat']);

trial = 1; % trial number to analyze

%% set some parameters

w = 0.5; % window length in seconds
n = round(w/stimulus(trial).frame); % window length in frames

Nx = stimulus(trial).param.x/stimulus(trial).param.dx; % number of independent stimulus regions in horizontal direction 
Ny = stimulus(trial).param.y/stimulus(trial).param.dy; % number of independent stimulus regions in vertical direction
m = Nx*Ny; %number of stimulus dimensions (640 for trial 1)

time = (0:stimulus(trial).Nframes-1)'*stimulus(trial).frame+stimulus(trial).onset; % stimulus timings
T = length(time)-1; % number of stimulus frames

%% first, count the number of spikes during each frame

spikeCounts = zeros(T,datainfo.Ncell);
for j=1:datainfo.Ncell
    % we'll count the number of spikes in each frame period using a histogram with bins set by the time vector
    tmp = histc(spikes{j,trial},time);
    tmp(1:n-1) = 0; % ignore the first "n" frames
    spikeCounts(:,j) = tmp(1:end-1);
    spikeCountsAvg(:,j) = tmp(1:end-1)/sum(tmp(1:end-1)); %can divide by the number of spike for "averaging" purpose later
end

%% now perform spike-triggered averaging using matrix multiplication to avoid for-loops
% this is the neuron-by-neuron way to run the analysis: it produces an m-by-n
% matrix corresponding to the spike-triggered-average input to one cell
% (specified by cellnum)

cellnum = 1;
STA = zeros(m,n);
spikeCountsShift = spikeCountsAvg(:,cellnum);

for j=1:n % we loop over time offsets    
    
    STA(:,j) = r*spikeCountsShift;
    
% you should convince yourself that the inner product above is equivalent to the following for-loop:
%     for t = 1:T
%         STA(:,j) = STA(:,j) + r(:,t) * spikeCountsShift(t,:);
%     end

    spikeCountsShift = circshift(spikeCountsShift,[-1 0]); %shift things one time bin into the past and repeat
    
end
%% or, perform spike-triggered averaging for ALL cells at once
% we can save even more time by converting the STA matrix for each cell
% into an m*n-by-1 vector, like so:

STA = zeros(m*n,datainfo.Ncell);
spikeCountsShift = spikeCountsAvg; %spikeCountsShift is a T-by-7 matrix now

for j=1:n % we loop over time offsets as before
    
    STA((1:m)+(n-j)*m,:) = r*spikeCountsShift;
    % r is m x T and spikeCountsShift is T-by-7,
    % so their inner product is m-by-7, and gives the STA for each of the
    % seven cells at time offset j.
    
    spikeCountsShift = circshift(spikeCountsShift,[-1 0]); %shift things one time bin into the past and repeat
    
end

%% visualization of the STAs computed in the previous cell
timeAxis = (.5-n:0)*stimulus(trial).frame;
s = (1:max(Nx,Ny))*min(stimulus(trial).param.dx,stimulus(trial).param.dy)*stimulus(trial).pixelsize/1000;
figure(1);clf;
colormap gray;
for j=1:datainfo.Ncell
    
    subplot(1,datainfo.Ncell,j);
    
    % note that we can extract back out the STA matrix for any individual
    % neuron by calling Matlab's "reshape" function on STA:
    imagesc(timeAxis,s,reshape(STA(:,j),[],n));
    
    if(j==1)
        xlabel('time [s]');ylabel('space [mm]'); %apply axis labels
    end
end
