%% NeuralBox1 
%
% Odelia Schwartz
% Partly based on Cold Spring Harbor tutorial by Chichilnisky.
% This is a tutorial for understanding linear filter
% properties of neurons, and generating neural spikes
% based on a linear filter followed by a simple nonlinearity.
% This is part one of two tutorials. We are going to first show
% how to generate neural model responses, which we will
% later probe and "guess"/figure out the properties.
%
% The model neuron includes one linear and one nonlinear component.
% followed by generating a Poisson spike train. We have made several 
% versions of model neurons.
%
% ----------------------------------------------------------------------

%% clear all variables and close all figures
clear all; close all;

%% 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);
plot(stimulus(1:1000))
set(gca, 'fontsize', 14);
xlabel('Time (msec)')
ylabel('Stimulus strength')
title('First 1 second of stimulus');

%% We're now going to simulate a model neuron
% For purposes of this demo, we constructed some model
% neurons. The function getLinear1 computes the linear
% response of our model neuron. (This is a function we 
% constructed and so know the answer; in a usual
% experiment with real neurons one would not know this!)

% We've made several versions of model neurons.
% Toggle between
filterLen = 100;
[linearResp, filter] = getLinear1(stimulus, filterLen);
%[linearResp, filter] = getLinear2(stimulus, filterLen);
%[linearResp, filter] = getLinear3(stimulus, filterLen);

%% The linear response smooths (averages) the stimulus
% over time.  The top panel of the figure shows 
% the first second of the stimulus, the small middle 
% panel shows the linear filter (which we usually would
% not know), and the third panel shows the first second 
% of the linear response.
figure(2);
subplot(2,1,1);
plot(stimulus(1:1000))
set(gca, 'fontsize', 14);
title('Stimulus');
set(gca,'XTIck',[])
subplot(2,1,2);
plot(linearResp(1:1000))
set(gca, 'fontsize', 14);
title('Linear response');
xlabel('Time (ms)');
axes('position',[.13 .47 .2 .1])
plot(filter,'m-')
axis([0 length(filter) 1.1*min(filter) 1.1*max(filter)]);
set(gca,'XTick',[],'YTick',[],'Box','on','Units','Normalized')

%% We would like to unpack this and see what the linear
% filter is doing
% Choose a random starting point, and examine the stimulus
% starting from the random position and with a length equal
% to the filter length. We'll plot both this stimulus sequence
% and the filter. 
thestart = round(rand*100+1);
thelen = length(filter)
thestim = stimulus(thestart:thelen+thestart-1)
subplot(3,1,1); plot(thestim); 
set(gca, 'fontsize', 14); title('stimulus')
subplot(3,1,2); plot(filter); 
set(gca, 'fontsize', 14); title('filter')

%% Plot the point by point multiplication of the filter
% with the stimulus
subplot(3,1,3); plot(filter.*thestim); 
set(gca, 'fontsize', 14); title('stimulus times filter');

%% The linear filter response to the stimulus, is the sum
% of this point by point multiplication. This results in a
% single number. This is also known as inner product or dot product.
sum(filter.*thestim)
linearResp(thelen+thestart-1)

%% What if the input was exactly the linear filter
% or variant of (toggle these)
thestim = filter;
% thestim = -filter;
% thestim = flipud(filter);
subplot(3,1,1); plot(thestim); 
set(gca, 'fontsize', 14); title('stimulus')
subplot(3,1,2); plot(filter); 
set(gca, 'fontsize', 14); title('filter')
subplot(3,1,3); plot(filter.*thestim); 
set(gca, 'fontsize', 14); title('stimulus times filter');
sum(filter.*thestim)

%% Under a simple non-linear model, the firing rate of a neuron is
% a single-valued non-linear function of an underlying linear response.
% We can pick any such function we want -- this is done in getNonlinear.
% Here we're applying this non-linear transformation on the linear response of
% our simulated neuron. 
nonlinearResp = getNonlinear1(linearResp);
% Toggle between
%nonlinearResp = getNonlinear2(linearResp);

%% We can plot together the linear and nonlinear response
h = figure(3); clf;
h = subplot(3,1,1);
plot(linearResp(1:1000),'b-')
title('Linear response')
set(h,'XTick',[],'XTickLabel',[])
h = subplot(3,1,2);
plot(nonlinearResp(1:1000),'r-')
title('Nonlinear function of linear response')
set(h,'XTick',[],'XTickLabel',[])

%% To understand this better, we try just few values
% Try changing theval: Is there a pattern?
theval = 220
linearResp(theval)
nonlinearResp(theval)

%% Plot nonlinear function in the interval [-.3:.05:.3]
% Toggle between
modelNonlin = getNonlinear1([-.3:.05:.3]);
%modelNonlin = getNonlinear2([-.3:.05:.3]);
figure(4); clf;
plot([-.3:.05:.3], modelNonlin);
set(gca, 'fontsize', 14);
xlabel('Linear response'); ylabel('Nonlinear response');

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

%% So far, we constructed a model neuron and its response to experimental
% stimuli. We first constructed random Gaussian stimulus, we linearly 
% filtered it, put this linearly filtered signal through a non-linear 
% function to calculate an underlying firing rate of the cell, and used 
% this to simulate spikes coming out of the cell.  
% Here's the first second of each of these:
h = figure(5); clf;
h = subplot(3,1,1);
plot(linearResp(1:1000),'b-')
title('Linear response')
set(h,'XTick',[],'XTickLabel',[])
h = subplot(3,1,2);
plot(nonlinearResp(1:1000),'r-')
title('Nonlinear function of linear response')
set(h,'XTick',[],'XTickLabel',[])
h = subplot(3,1,3);
stem(neuralResponse(1:1000))
set(gca,'Ylim',[0 2]);
title('# of Spikes (1 ms bins)');
set(h,'XTick',[],'XTickLabel',[])

%% Things to do in class:
%
% 1. Try changing the linear function (choosing between getLinear1,
% getLinear2, getLinear3; see toggle comment). 
%
% 2. Try changing the nonlinear function (choose between getNonlinear1,
% getNonlinear2.m). Change it in the two places in the tutorial.
% What is the difference between getNonlinear1 and getNonlinear2?
