% generate spikes:
n = 1000;
misi = 0.01;
isis = exprnd( misi, n, 1 );
isis = isis + 0.01;
spikes = cumsum( isis );
p = rand( size( spikes ) );

% spike waveforms:
dt = 0.0001;
x = -0.01:dt:0.01;
y1 = 2.0*exp( -(x-0.0003).^2/2.0/0.0005^2 ) - 1.4*exp( -(x-0.0005).^2/2.0/0.0005^2 );
y2 = exp( -x.^2/2.0/0.002^2 ).*cos(2.0*pi*x/0.005);
p12 = 0.5;

% voltage trace:
noise = 0.2;
t = 0:dt:spikes(end)+3.0*x(end);
v = noise*randn( 1, length( t ) );
for k = 1:length( spikes )
    inx = ceil( spikes(k)/dt );
    if p(k) < p12
        v(inx:inx+length(x)-1) = v(inx:inx+length(x)-1) + y1; 
    else
        v(inx:inx+length(x)-1) = v(inx:inx+length(x)-1) + y2; 
    end
end

figure( 1 );
clf;
plot( t(t<1.0), v(t<1.0) );
hold on;

% find peaks:
thresh = 0.7;
inx = find( v(2:end-1) > thresh & v(1:end-2)<v(2:end-1) & v(2:end-1) > v(3:end) ) + 1;
spiket = t(inx);
spikev = v(inx);
tinx = inx;
for k=1:2
    inx = find( ( spiket(2:end-1)-spiket(1:end-2)>0.005 | spikev(2:end-1) > spikev(1:end-2) ) & ( spiket(3:end)-spiket(2:end-1)>0.005 | spikev(2:end-1) > spikev(3:end) ) )+1;
    spiket = spiket(inx);
    spikev = spikev(inx);
    tinx = tinx(inx);
end
scatter( spiket(spiket<1.0), spikev(spiket<1.0), 100.0, 'r', 'filled' );
%scatter( t(tinx), v(tinx), 100.0, 'r', 'filled' );
hold off;

% spike waveform snippets:
w = ceil( 0.005/dt );
vs = [];
for k=1:length(tinx)
    vs = [ vs; v(tinx(k)-w:tinx(k)+w-1) ];
end
ts = t(1:size(vs,2));
ts = ts - ts(floor(length(ts)/2));
figure( 2 );
clf;
hold on
for k=1:size(vs,1)
    plot( ts, vs(k,:), '-b' );
end
hold off

% pca:
cv = cov( vs );
[ ev , ed ] = eig( cv );
[d,dinx] = sort( diag(ed), 'descend' );
% spectrum of eigenvalues:
figure( 3 );
clf;
subplot( 4, 1, 1 );
scatter( 1:length(d), d, 'b', 'filled' );
% project onto eigenvectors:
nx = vs * ev(:,dinx(1));
ny = vs * ev(:,dinx(2));
%scatter( nx, ny, 'b', 'filled', 'MarkerEdgeColor', 'white' );
%hold on;

% features:
subplot( 4, 2, 5 );
plot( 1000.0*ts, ev(:,dinx(1)), 'g', 'LineWidth', 2 );
subplot( 4, 2, 6 );
plot( 1000.0*ts, ev(:,dinx(2)), 'r', 'LineWidth', 2 );

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

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

% 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),:), '-g' );
end
plot( 1000.0*x, y1, '-k', 'LineWidth', 2 );
xlim( [ 1000.0*ts(1) 1000.0*ts(end) ] )
hold off
subplot( 4, 2, 8 );
hold on
kinx2 = find(kx==2);
for k=1:length(kinx2)
    plot( 1000.0*ts, vs(kinx2(k),:), '-r' );
end
plot( 1000.0*x, y2, '-k', 'LineWidth', 2 );
xlim( [ 1000.0*ts(1) 1000.0*ts(end) ] )
hold off

% spike trains:
figure( 1 );
hold on;
six1 = find( spiket(kinx1)<1.0 );
scatter( spiket(kinx1(six1)), spikev(kinx1(six1)), 100.0, 'g', 'filled' );
six2 = find( spiket(kinx2)<1.0 );
scatter( spiket(kinx2(six2)), spikev(kinx2(six2)), 100.0, 'r', 'filled' );
hold off;