% load data into time, voltage and spiketimes
load( 'extdata' );

% indices into voltage trace of spike times:
dt = time( 2) - time(1);
tinx = round(spiketimes/dt)+1;

% plot voltage trace with dettected spikes:
figure( 1 );
clf;
plot( time, voltage, '-b' )
hold on
scatter( time(tinx), voltage(tinx), 'r', 'filled' );
xlabel( 'time [ms]' );
ylabel( 'voltage' );
hold off

% spike waveform snippets:
w = ceil( 0.005/dt );
vs = [];
for k=1:length(tinx)
    vs = [ vs; voltage(tinx(k)-w:tinx(k)+w-1) ];
end
ts = time(1:size(vs,2));
ts = ts - ts(floor(length(ts)/2));

% % plot snippets:
figure( 2 );
clf;
hold on
for k=1:size(vs,1)
    plot( ts, vs(k,:), '-b' );
end
xlabel( 'time [ms]' );
ylabel( 'voltage' );
hold off

% pca:
cv = cov( vs );
[ ev , ed ] = eig( cv );
[d,dinx] = sort( diag(ed), 'descend' );

figure( 3 );
clf;
subplot( 4, 2, 1 );
imagesc( cv );
xlabel( 'time bin' );
ylabel( 'time bin' );
title( 'covariance matrix' );
caxis([-0.1 0.1])

% spectrum of eigenvalues:
subplot( 4, 2, 2 );
scatter( 1:length(d), d, 'b', 'filled' );
xlabel( 'index' );
ylabel( 'eigenvalue' );

% features:
subplot( 4, 2, 5 );
plot( 1000.0*ts, ev(:,dinx(1)), 'r', 'LineWidth', 2 );
xlabel( 'time [ms]' );
ylabel( 'eigenvector 1' );
subplot( 4, 2, 6 );
plot( 1000.0*ts, ev(:,dinx(2)), 'g', 'LineWidth', 2 );
xlabel( 'time [ms]' );
ylabel( 'eigenvector 2' );

% project onto eigenvectors:
nx = vs * ev(:,dinx(1));
ny = vs * ev(:,dinx(2));
%scatter( nx, ny, 'b', 'filled', 'MarkerEdgeColor', 'white' );

% clustering (two clusters):
%kx = kmeans( [ nx, ny ], 2 );
% nx smaller or greater a threshold:
kthresh = 1.6;
kx = ones( size( nx ) );
kx(nx<kthresh) = 2;

subplot( 4, 1, 2 );
scatter( nx(kx==1), ny(kx==1), 'r', 'filled', 'MarkerEdgeColor', 'white' );
hold on;
scatter( nx(kx==2), ny(kx==2), 'g', 'filled', 'MarkerEdgeColor', 'white' );
hold off;
xlabel( 'projection onto eigenvector 1' );
ylabel( 'projection onto eigenvector 2' );

% show sorted spike waveforms:
subplot( 4, 2, 7 );
hold on
kinx1 = find(kx==1);
for k=1:length(kinx1)
    plot( 1000.0*ts, vs(kinx1(k),:), '-r' );
end
plot( 1000.0*waveformt, waveform2, '-k', 'LineWidth', 2 );
xlim( [ 1000.0*ts(1) 1000.0*ts(end) ] )
xlabel( 'time [ms]' );
ylabel( 'waveform 1' );
hold off

subplot( 4, 2, 8 );
hold on
kinx2 = find(kx==2);
for k=1:length(kinx2)
    plot( 1000.0*ts, vs(kinx2(k),:), '-g' );
end
plot( 1000.0*waveformt, waveform1, '-k', 'LineWidth', 2 );
xlim( [ 1000.0*ts(1) 1000.0*ts(end) ] )
xlabel( 'time [ms]' );
ylabel( 'waveform 2' );
hold off

% spike trains:
figure( 1 );
hold on;
scatter( time(tinx(kinx1)), voltage(tinx(kinx1)), 100.0, 'r', 'filled' );
scatter( time(tinx(kinx2)), voltage(tinx(kinx2)), 100.0, 'g', 'filled' );
hold off;