%% Part 1, computing an STA

%remove figures, clear variables, clear command window

close all;

clear all;

clc;



%load .mat datafile

load staData;



%1) look at spikes, find spike at 3ms

%use figure('position', [100 100 300 600]) to change size/location

figure;

stem(1:100, staData.n1Spikes(1:100))

xlabel('time (miliseconds)');

title('Spike Raster');



%unpack random stimuli from a given seed

rng(1337);

staData.stim = (2.*rand(33,33,10000)-1);



%2) display stim at 3ms, can see nothing

figure;

imagesc(staData.X, staData.Y, staData.stim(:,:,3));

xlabel('Horizontal Space');

ylabel('Vertical Space');

colormap(gray);

colorbar;

axis xy;

title('Stimuli');



%3) generate STA

sta = zeros(length(staData.X), length(staData.Y));

for t = 1:10000 %for every time bin

    if staData.n3Spikes(t) == true %check if a spike occured

        sta = sta + staData.stim(:,:,t); % if so, add stim to STA

    end %otherwise, do nothing

end %loop until t = 10000 then stop



%normalize sta

sta = sta / sum(staData.n3Spikes);



%4) display sta

figure;

imagesc(staData.X, staData.Y, sta);

xlabel('Horizontal Space');

ylabel('Vertical Space');

colormap(gray);

colorbar;

axis xy;

title('STA');



%% Part 2, generating a 2D RF



%make 2D arrays of X and Y positions

[DX, DY] = meshgrid(staData.X, staData.Y);



%inline function definition

%you can compose multiple inline functions into simple expressions

g = @(D,mu,sigma,alpha) alpha*exp(-(D-mu).^2./(2*sigma^2));

g2 = @(D1, D2, m1, m2, s, a) a*g(D1, m1, s, 1).*g(D2, m2, s, 1);



%make 2d RF

rf = g2(DX, DY, 0, 5, 2, 2);



%display RF

imagesc(staData.X, staData.Y, rf);

colorbar;



%% Part 3, simulating a LNP neuron



%create raster variable

mySpikes = false(1,10000);



%for every time bin

for t = 1:10000

    

    %take the dot product of the stim and RF

    r = sum(sum(staData.stim(:,:,t).*rf));

    

    %pass through nonlinearity

    lambda = max(0,r);

    

    %coinflip for spike

    if rand < lambda

        mySpikes(t) = true;

    end

end



%plot raster

figure;



%subplot puts multiple plots in same figure window

subplot(2,1,1);

stem(1:100, mySpikes(1:100))

xlabel('time (miliseconds)');

title('Spike Raster');



%see how many spikes occured, divided by duration in seconds.

%can scale down RF gain to compensate

spkCount = sum(mySpikes)/10;

disp(['Rate: ' num2str(spkCount), ' Hz']);



% fancy vectorized sta code

subplot(2,1,2);

imagesc(staData.X, staData.Y, mean(staData.stim(:,:,mySpikes),3));