%% Neural Box 
% Odelia Schwartz
% Partly based on Cold Spring Harbor white noise tutorial by Chichilnisky.
% This is a tutorial for figuring out the response properties
% of a neuron by presenting many random stimuli.
% The neuron here is a model simulation, but the
% same approach is used for understanding real neurons.
% We probe the model neuron as if it is a black box and
% we do not know its properties.
% This tutorial specifically focuses on spike-triggered average,
% estimating the linear filter, and estimating the nonlinearity.
%
% At the end of the tutorial, there are some questions.
% ----------------------------------------------------------------------

%% We want to first choose experimental stimuli that are random
% At each frame, the intensity of the uniform screen changes:
% it is drawn randomly, here from a Gaussian distribution.
numsamples = 25000
stimulus=(1/3*randn(numsamples,1))

%% Plot the stimulus
figure(1);
thelen = min(1000, numsamples);
plot(stimulus(1:thelen))
set(gca, 'fontsize', 14);
xlabel('Time (msec)')
ylabel('Stimulus strength')
title('First 1 second of stimulus');

%% Let's check that it is Gaussian
figure(2);
[H, x] = hist(stimulus, 20);
bar(x, H./sum(H));
set(gca, 'fontsize', 14);
xlabel('Stimulus strength');
ylabel('Histogram');

%% We're now going to simulate a model neuron
% For purposes of this demo, we are constructing the model
% neurons and so know their filters and nonlinearity
% (in an experiment with real neurons, we would be handed 
% the spike trains and would not know this!)

% We've made several versions of model neurons.
% We have 3 possible linear filters.
% Toggle between
kernelSize = 60;
[linearResp, filter] = getLinear1(stimulus, kernelSize);
%[linearResp, filter] = getLinear2(stimulus, kernelSize);
%[linearResp, filter] = getLinear3(stimulus, kernelSize);

%% Let's look at the filter
figure(3);
h=plot(filter, 'o-');
set(h, 'linewidth', 2);
set(gca, 'fontsize', 14);
title('actual filter')
xlabel('time');

%% We also have two versions of nonlinearities.
nonlinearResp = getNonlinear1(linearResp);
% Toggle between
%nonlinearResp = getNonlinear2(linearResp);

%% We can use this non-linear response to simulate a 
% Poisson-ish spike train... 
xr=rand(size(nonlinearResp));
neuralResponse = nonlinearResp>.05*xr;
spikeCounts = neuralResponse;

%% So far, we constructed a model neuron and its response to experimental
% stimuli. Here's the first second of each step:
thelen = min(numsamples, 1000);
h = figure(4); clf;
h = subplot(3,1,1);
plot(linearResp(1:thelen),'b-')
title('Linear response')
set(h,'XTick',[],'XTickLabel',[])
h = subplot(3,1,2);
plot(nonlinearResp(1:thelen),'r-')
title('Nonlinear function of linear response')
set(h,'XTick',[],'XTickLabel',[])
h = subplot(3,1,3);
stem(neuralResponse(1:thelen))
set(gca,'Ylim',[0 2]);
title('# of Spikes (1 ms bins)');
set(h,'XTick',[],'XTickLabel',[])

%% Now we compute the spike-triggered average stimulus.  This is accomplished
% by taking the 60 milliseconds of stimulus immediately preceding each spike
% and adding them together.  This sum is then divided by the total number 
% of spikes fired over the course of the entire experiment to determine the 
% average stimulus preceding a spike.
% This spike-triggered average is, in a sense, a template for what the neuron
% is "looking for" in the stimulus.

totalCount=sum(spikeCounts)
a=zeros(kernelSize, 1);
for i=kernelSize+1:length(spikeCounts)
  a = a + spikeCounts(i)*stimulus((i-kernelSize+1):i);
end
a=a/totalCount;

%% We'll first look at the answer; then unpack what we did
figure(5);
clf;
h=plot(a, 'o-');
set(h, 'linewidth', 2);
set(gca, 'fontsize', 14);
xlabel('time'); ylabel('estimated linear filter')

%% Because this is a tutorial, we *know* exactly what filtering 
% was done on the stimulus to get the linear response ("linearResp"). 
% Below, we compare the spike-triggered average to the filter we used.
% They are similar shape up to a constant multiplication factor
figure(6);
clf;
h=plot([a/abs(sum(a)),filter/abs(sum(filter))], 'o-');
set(h, 'linewidth', 2);
set(gca, 'fontsize', 14);
legend('estimated filter', 'actual filter', 'Location', 'NorthWest')

%% Remember we summed together stimuli that led to a spike
% We can look at individual such stimuli and the average 
% as we have more samples. This code is slower, so we'll
% do steps of 100
atmp = 0; totalTmp = 0;
figure(7)
thelen = min(numsamples, 10000);
for i=kernelSize+1:thelen
  atmp = atmp + spikeCounts(i)*stimulus((i-kernelSize+1):i);
  totalTmp = totalTmp+spikeCounts(i);
  if spikeCounts(i)==1 & mod(totalTmp, 100)==0
    subplot(2,1,1);
    h=plot(stimulus((i-kernelSize+1):i));
    set(gca, 'fontsize', 14); set(h, 'linewidth', 2);
    title('stimulus that resulted in spike')
    atmp = atmp + stimulus((i-kernelSize+1):i);
    subplot(2,1,2);
    h=plot(atmp); set(h, 'linewidth', 2);
    set(gca, 'fontsize', 14);
    title('average stimulus that resulted in spike')
    pause(.5);
  end
end

%% Extra intuition
% To get intuition about why averaging the spiked stimuli
% works, we can look at how the linear response relates to
% spikes versus no spikes. This allows us to differentiate 
% between stimuli that lead to spikes or no spikes.
figure(8)
linearEst=zeros(size(spikeCounts));
for i=kernelSize+1:length(spikeCounts)
  linearEst(i) = a' * stimulus((i-kernelSize+1):i);
end
plot(linearEst,spikeCounts,'o');
set(gca, 'fontsize', 12);
xlabel('Linear response');
ylabel('Spikes');
title('Relation between estimated linear responses and spikes')

%% Extra
%% We can also actually estimate the nonlinearity.
% We can plot the *average* number of spikes fired in response
% to similar linear responses.  First we decide on linear
% response bins...
figure(9);
x=[-.2:.05:.3];
% ...and then calculate the averages.
N=zeros(size(x));
se=zeros(size(x));
for i=1:length(x)
  L = logical((linearEst>x(i)) & (linearEst<(x(i)+.05)));
  N(i)=mean(spikeCounts(L));
  se(i)=sqrt(var(spikeCounts((linearEst>x(i)) & (linearEst<(x(i)+.05))))/sum(L));
end
x=x+.025;
errorbar(x,N,se,'ro')
set(gca,'xLim',[min(x), max(x)]);
xlabel('Linear response component');
ylabel('Mean spike count');

% We can compare this to the nonlinearity that we know because we constructed
% the simulation - but usually would not know.
% That is, we can superimpose the non-linear function that we actually used to
% determine the spike firing probabilities.  The plot of linear response versus
% mean spike count should have the same shape as this function, but, there
% is an arbitrary scale factor relating these two quantities.  Below, we estimate
% this scale factor using least-squares.
xx = getNonlinear1(x);           % this we know
gamma = inv(xx*xx') * xx * N';   % find the scale factor
vals=[-.4:.05:.4];
Nth =  getNonlinear1(vals);
hold on
plot(vals,Nth*gamma)
hold off

%% Questions
% Each time make only one change to the code:
%
% 1. Try lowering the numsamples in the code above
% (example, to 200, 700, and 2000)? Plot the 
% filter estimates for these. When does the filter 
% estimate break down? Why does it break down?
%
% 2. Try changing the linear function (choosing between getLinear1,
% getLinear2, getLinear3; see toggle comment). Can we recover each of
% the linear filters properly?
%
% 3. Try changing the nonlinear function (choose between 
% getNonlinear1, getNonlinear2.m) and keeping getLinear1. 
% Can we recover the linear filter of the neuron for each 
% of these nonlinearities? If not, then why?
%
% Extra:
% 4. What is special about the stimuli that lead to a spike?
% What would happen if we tried to average random stimuli
% without knowing which ones lead to a spike? To do so, you can
% try to create a false spike train each time (for instance do help on
% randperm, and try to permute spikeCounts) and compute the spike-triggered
% average. 

