
%%% This is a tutorial for linear visual filters (receptive fields)
%%% and spike-triggered covariance.
%%% (1) Examples of constructing visual Gabor filters and
%%% filtering an image.
%%% (2) Examples of spike-triggered aproaches to find filters.

%%% Read the comments and copy each line of code into matlab.
%%% Type 'help <function_name>' in Matlab window for any
%%% function you would like more information on.
%%% Type 'which <function_name>' in Matlab window for the
%%% location of a file.

%%% Odelia Schwartz 2012
% This is a simplified version of:
% Spike-triggered neural characterization
% Schwartz, Pillow, Rust, Simoncelli
% Journal of Vision, 2006

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% (1) Visual filters and images
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% (1a) Gabor filters and images

%% Set parameters of sinusoid
sz = 20;
ctr = [round(sz/2+1), round(sz/2+1)];
period = 5;
direction = pi;
phase = 2*pi;
theSine = mkSine(sz, period, direction, 1, phase);
%%% Plot the sinusoid
subplot(2,2,1); showIm(theSine)

%% Make a 2 dimensional Gaussian and plot it
thesig = 2;
theGauss = mkGaussian(sz, thesig, ctr,1);
subplot(2,2,2); showIm(theGauss)

%% Make a Gabor filter, by multiplying a sinusoid with a Gaussian.
theFilt = theSine .* theGauss;
subplot(2,2,3); showIm(theFilt)

%%% 2. Load and filter an image

%% Load an image
im = pgmRead('einstein.pgm');
figure(2); subplot(2,2,1); showIm(im);

%% (1b) Convolve the image with the filter
response = rconv2(im, theFilt); 
figure(2); subplot(2,2,2); showIm(response)

%%% *To do*:
%%% 1. Try making different Gabor filters by varying parameters above
%%% (e.g., direction, priod, phase of the grating; and thesig of the Gaussian)

%%% (2) Spike-triggered approaches
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Close all figures
close all; 

%%% We have constructed in advance few model neurons. We will use spike-triggered
%%% approaches to figure out the receptive field properties of the neurons.

%% 2a. Generate random stimuli to "probe" the neuron with
nFrames = 500000;
xDim = 8;
kernelX = xDim;                 % spatial size of noise stimulus
kernelT = 6;                    % temporal size of noise stimulus
kernelSize = kernelX * kernelT;
allStim = randn(nFrames, kernelSize);

%% Show example frames of the white noise stimuli
figure(1); clf;
for theframe=1:16
subplot(4,4,theframe); showIm(reshape(allStim(theframe,:),6,8))
end

%% 2b. Generate spikes from a model neuron
%%% This can be toggled for different model neurons; choose from:
%%% spikeResp = ClassModel1(allStim);
%%% spikeResp = ClassModel2(allStim);
%%% spikeResp = ClassModel3(allStim);
spikeResp = ClassModel1(allStim);

%% Plot the spiking activity for the first 100 frames
clf; stem(spikeResp(1:100),'o')

%% Compute the spike-triggered average

%%% First find the frames for which the model neuron spiked
spikeInd=(find(spikeResp>0.5));

%% Then find the spike-triggered stimuli, i.e., the stimuli for which
%%% the neuron spiked
spikeStim = allStim(spikeInd,:);
numspikes = length(spikeInd);

%% Plot some example stimulus frames of the spike-triggered stimuli
%%% Can you tell by eye what in the stimulus is triggering a spike?
figure(1); clf;
for theframe=1:16
subplot(4,4,theframe); showIm(reshape(spikeStim(theframe,:),6,8))
end

%% We'll plot the spike-triggered average (STA)
%%% Is it a structured receptive field?
figure(2); clf;
sta = mean(spikeStim);
showIm(reshape(sta,6,8));
title('Spike-triggered average')

%%% Are there other receptive fields this model neuron is using to compute its
%%% spikes?

%% We'll do a simple version of a spike-triggered covariance
%%% This is a Principal Component Analysis, computing the eigenvalues
%%% (variances along each receptive field axes) and the eigenvectors
%%% (the receptive field axes).
thecov = (spikeStim'*spikeStim)./(numspikes-1);

[eigvec,eigval] = sortedEig(thecov);

%%% Plot the (sorted) eigenvalues 
%%% How many appear significant?
figure(3); clf;
plot(diag(eigval), 'o');

%% Plot the corresponding eigenvector (here set to the first).
%%% Note: If the STA was structured, the last eigenvector could just be the 
%%% first eigenvector, possibly negated 
thenum1 = 1;
figure(4); clf;
subplot(2,2,1);
showIm(reshape(eigvec(:,thenum1),6,8))
title(sprintf('Eigenvector %d',thenum1));

%% Plot another eigenvector
%%% Here set to the second, but change as needed...
%%% Is it structured? Do we expect it to be based on the eigenvalues?
thenum2 = 48;
figure(4);
subplot(2,2,2);
showIm(reshape(eigvec(:,thenum2),6,8))
title(sprintf('Eigenvector %d',thenum2));

%% Look at scatter plots onto two eigenvectors or receptive fields.
%%% We will compare the responses to the spike-triggered stimuli with
%%% those to the full stimulus set. We will match the number of stimuli
%%% for readability of the plots.

%%% The two receptive fields
%%% Toggle appropriate choices
basis1 = eigvec(:,thenum1); 
basis2 = eigvec(:,thenum2); 
%%% Responses of the two receptive fields to all stimuli
allProj = [allStim*basis1, allStim*basis2];
%%% And to the spike-triggered stimuli
spikeProj = [spikeStim*basis1, spikeStim*basis2];

figure(5); clf;
thenum = min(20000, numspikes);
plot(allProj(1:thenum,1), allProj(1:thenum,2), 'bo')
hold on;
plot(spikeProj(1:thenum,1), spikeProj(1:thenum,2), 'ro');
axis square;
xlabel (sprintf('Receptive field %d',thenum1));
ylabel(sprintf('Receptive field %d',thenum2));

%% Plot ellipse signifying the variances found by the Principal Component Analysis
d=diag(eigval);
angles=2*pi*[0:30]'/30;
%%% Variance along the 2 receptive fields
ellipse = 3*[sqrt(d(thenum1))*cos(angles), sqrt(d(thenum2))*sin(angles)];
%%% Variance along 2 other directions that are not structured
ellipse_other = 3*[sqrt(d(24))*cos(angles), sqrt(d(25))*sin(angles)];

%% Plot the ellipses
plot(ellipse_other(:,1),ellipse_other(:,2),'c','lineWidth',2);
plot(ellipse(:,1),ellipse(:,2),'m','lineWidth',2); 
axis square
hold off;
axis([-5 5 -5 5])
legend('All stim', 'Spike stim', 'All stim', 'Spike stim')

%%% What are the properties of this neuron?
%%% Go through this code with the remaining model neurons above.


