%% Exercise 1: Matrix Operations for a Feedforward Network
close all
clear all
clc

M = 10;
N = 50;

W = rand( M, N );
W = 2.*( W - 0.5 );

x = randn( 1, 50 )';

y = W * x;

figure
subplot( 3, 1, 1 )
plot( x )
box off
xlabel( 'sample number' )
ylabel( 'input values (a.u.)' )

subplot( 3, 1, 2 )
pcolor( W )
title( 'Weights matrix' )

colorbar
subplot( 3, 1, 3 )
plot( y )
box off
xlabel( 'sample number' )
ylabel( 'output values (a.u.)' )

%% Exercise 2: Logical operations, for-loops, and plotting (random walk)
clear all
close all
clc

T = 10000; % measured in ms
Vreset = -70; % measured in mV
Vthresh = -45; % measured in mV
V0 = -65; % measured in mV

% find p such that the firing rate is approximately 10Hz
firingRate = 0;
p = 0;
while 1
if abs( firingRate - 10 ) < 0.1
break

end

V = generatevoltage( p, T, Vreset, Vthresh, V0 );
firingRate = sum( V == Vthresh ) / ( T / 1000 );
p = p + 0.0001;

if p > 1
disp( 'did not find a suitable p' )
break

end

end

ts = 1 : T; % time vector in ms

figure
plot( ts, V )
box off
xlabel( 'time (ms)' )
ylabel( 'voltage (mV)' )
xlim( [ 0 1000 ] )
ylim( [ -80 -40 ] )
% title( [ 'probability that dV = 1 is ', num2str( p ) ] )

%% Exercise 3: Convolution to estimate voltage response to a spike train
clear all
close all
clc

% generating a 20 Hz Poisson Spike Train
targetRate = 20; % Measured in Hz.
N = 10000; % Measured in ms.
firingRate = 0;
p = 0;

while 1
if abs( firingRate - targetRate ) < 0.1
break
end

spikeTrain = rand( 1, N ) > ( 1 - p );
firingRate = sum( spikeTrain == 1 ) / ( N / 1000 );
p = p + 0.00001;

if p > 1
disp('did not find a suitable p')
break

end

end

% construct kernel
mu = 5;
k = exppdf( -50 : 50, mu );

% Find spike times first
ts = 1 : N; % time vector in ms
spkTimes = ts( spikeTrain );

% convolution
spikeTrain = double( spikeTrain );
estVoltage = conv( spikeTrain, k, 'same' );

% plot
figure
subplot( 2, 1, 1 )
plot( -50 : 50, k )
box off
title( 'kernel' )
xlabel( 'tau (ms)' )

subplot( 2, 1, 2 )
plot( ts, estVoltage )
hold on
yLevel = 0.3;
plot( spkTimes, yLevel, 'k', 'Marker', '*', 'MarkerSize', 4 )
box off
ylabel('voltage (a.u.)')
xlabel( 'time (ms)' )
xlim( [ 0 1000 ] )
ylim( [ 0 0.4 ] )
legend( 'voltage', 'spike times' )

%% Exercise 4: Convolution to detect edges in images
clear all
close all
clc

% Load image.
load( 'octopus.mat' )

figure
imagesc( octopus );
colormap gray

% Make the kernel.
k = [...
0, 0, 0;...
0, 1.125, 0 ;...
0, 0, 0 ] - 0.125 * ones( 3, 3 );
% k = double( k );
octopusEdges = conv2( octopus, k, 'same' );

figure
imagesc( octopusEdges );
colormap gray

edgeImage = abs( octopusEdges );

figure
imagesc( edgeImage );
colormap gray
caxis( [ -20 100 ] )

%% Exercise 5: Correlation to analyze premotor neural data
clear all
close all
clc

load( 'MackeviciusData.mat' )

T = length( song ) / fs;
dt = 1 / fs;
k = length( song );
ts = dt : dt : T;

figure
subplot( 2, 1, 1 )
plot( ts, song )
title( 'Bird song.' )
ylabel( 'Amplitude (mV)')

subplot( 2, 1, 2 )
plot( ts, units )
title( 'Units firing.' )
ylabel( 'Amplitude (\muV)')

% Play the song if you wish.
% sound( song, fs )

% Smooth song and spike trains.
gWin = gausswin( ceil( fs * 0.026 ) );
smoothUnits = 10 * log10( conv( units .^ 2, gWin, 'same') );
smoothSong = 10 * log10( conv( song.^2, gWin, 'same' ) );

figure
subplot( 2, 1, 1 )
plot( ts, smoothSong )
ylim( [ -30 10 ] )
title( 'Smoothed bird song power.' )
ylabel( 'Power (dB)')

subplot( 2, 1, 2 )
plot( ts, smoothUnits )
ylim( [ 0 7 ] )
title( 'Units firing power.' )
ylabel( 'Power (dB)')
xlabel( 'time (s)' )

% Compute the cross-correlation between the two signals. Limit the lags to
% [-0.5, 0.5] seconds.
maxLag = 0.5;
lagSamp = round( maxLag * fs );
[ songUnits, tauSamp ] = xcorr( smoothUnits, smoothSong, lagSamp );

% Generate tau lags for plotting.
% tau = ( 1 : length( songUnits ) ) .* dt - 0.5 * length( songUnits ) * dt;
tau = tauSamp ./ fs;

figure
plot( tau, songUnits )
hold on
xlim( [ -0.5 0.5 ] )
xlabel( '\tau (s)' )
ylabel( 'Corr. Value (a.u.)' )

% Find time at which units fire in relation to song production.
[ pks, locs ] = findpeaks( songUnits );
[ ~, maxIdx ] = max( pks );
maxTau = tau( locs( maxIdx ) );
fprintf( 'Unit firing precedes song production by %5.2g s.', -maxTau );

% Mark max tau.
plot( maxTau, pks( maxIdx ), 'ro',...
'MarkerFaceColor', 'r',...
'MarkerSize', 6 )

%% Exercise 6: Singular Value Decomposition (SVD) on images of faces
clear all;
close all
load( 'jaffe.mat' )

% i=1;
% figure()
% imagesc(reshape(IMS(i,:),137,86));
% colormap gray;

meanImages = mean( IMS, 2 );
IMS = IMS - repmat( meanImages, 1, size( IMS, 2 ) );
for i = 1 : size( IMS, 1 )
IMS( i, : ) = IMS( i, : ) ./ var( IMS( i, : ) );

end
[ U, S, V ] = svd( IMS );

% confirm IMS = U*S*V'

figure
plot( diag( S ) )
xlabel( 'ranking singular value label' )
ylabel( 'singular value' )

figure
for i = 1 : 9
subplot( 3, 3, i )
imagesc( reshape( V( :, i ), 137, 86 ) );
colormap gray
title( [ 'singular value ', num2str( i ) ] )
axis equal

end

n = 9;
VRedDim = V;
for i = ( n + 1 ): size( V, 2 )
VRedDim( :, i ) = 0;

end
I = U * S * VRedDim';

numImages = 5;
figure
count = 1;
for j = 1 : numImages
subplot( numImages, 2, count )
count = count + 1;
imagesc( reshape( IMS( j, : ), 137, 86 ) );
axis equal
colormap gray
subplot( numImages, 2, count )
count = count + 1;
imagesc( reshape( I( j, : ), 137, 86 ) );
colormap gray
axis equal

end

%% Exercise 7: Estimation of spectro-temporal receptive fields with spike-triggered averaging
clear all
close all
clc

stimulus = generatestimulus;

freq = logspace( 2, 4, 50 );% 50 log-spaced freqs 100 - 10,000 Hz
dt = 0.01; % sampling interval in seconds

% Listen to the first 3 seconds of the stimulus.
% playstim( stimulus( :, 1 : ( 3 / dt ) ), freq, dt);

kernel = generatekernel;
% playstim( Kernel, freq, 0.001 )

estResponse = zeros( 1, size( stimulus, 2 ) );
for i = size( kernel, 2 ) + 1 : size( stimulus, 2 )
response = kernel .* stimulus( :, i - size( kernel, 2 ) : i - 1 );
estResponse( i ) = sum( response( : ) );

end

sortedEstResp = flip( sort( estResponse ) );
thresh = sortedEstResp( 1e4 );
spikes = estResponse >= thresh;
spikes = double( spikes );

% Plot variables.
figure
subplot( 3, 1, 1 )
imagesc( stimulus( :, 1 : ( 5 / dt ) ) )

subplot( 3, 1, 2 )
plot( estResponse( 1 : ( 5 / dt ) ) )

subplot( 3, 1, 3)
plot( spikes( 1 : ( 5 / dt ) ) )

% Find STA
sta = zeros( 50, 150 );
for i = find( spikes == 1 )
sta = sta + ( 1 / sum( spikes ) ) .* stimulus( :, i - 150 : i - 1 );

end

figure
surf( 1 : 150, logspace( 2, 4, 50 ), sta, 'edgecolor', 'none' );
axis tight;
view( 0, 90 );
shg
set( gca, 'yscale', 'log' )
xlabel( 'time (ms)' )
ylabel( 'Frequency (Hz)' )
title( 'STA' )
colorbar;

figure
imagesc( kernel )
title( 'Kernel (STRF)' )

% playstim( sta, freq, 0.001 )