Matlab Code for Spike Time Distances, Parallel in q
main cost-based metrics page
algorithm page for cost-based metrics
Matlab code for Spike Time Metric, Parallel in q
There are three modules: spkdallq_recur, spkdallqk_final, spkdallqk_dist,
with the following relationships:
- spkdallq_recur
carries out the core dynamic programming algorithm and produces
two output arrays, clast and c. clast applies to the entire spike trains; c can be
used to calculate distances from truncated spike trains.
- spkdallq_final
carries out the final step of the algorithm to calculate distances
from the array clast produced by spkdallq_recur.
- spkdallq_dist
calls spkdallq_recur and then spkdallq_final.
spkdallq_recur
<![CDATA[
function [clast,rmax,c]=spkdallq_recur(sa,sb)
% function [clast,rmax,c]=spkdallqk_recur(sa,sb) does the recursion for the DP
% algorithm for spike time distance. It uses the method of spkdallq, but
% reformats and produces arguments to be parallel to spkdallqk_recur.
%
% sa: vector of spike times for first spike train
% sb: vector of spike times for second spike train
% c: array of critical lengthts of linkages. size(c)=[1+length(sa) 1+length(sb) 1+rmax]
% rmax: maximum number of linkages (=min(length(sa),length(sb))
% clast: vector of critical lengths for full spike trains, necessary for calculating all
% distances, which is done in spkdallq_dist, via spkdallq_final.
%
% Copyright (c) by Jonathan Victor.
%
% See also SPKD, SPKDALLQ, SPKDALLQK_RECUR.
%
tli=sa; %reformat the input arguments
tlj=sb;
nspi=length(tli);
nspj=length(tlj);
nij=min(nspi,nspj);
lc=repmat(NaN,[nspi+1 nspj+1 nij]); % assume no length of movement
if (nij>0)
%
% best linkage length is a recursion based on linkages for shorter trains
%
for i=2:nspi+1
for j=2:nspj+1
td=abs(tli(i-1)-tlj(j-1));
li=squeeze(lc(i-1,j,:));
lj=squeeze(lc(i,j-1,:));
lij=td+[0;squeeze(lc(i-1,j-1,1:nij-1))];
lc(i,j,:)=min([li,lj,lij],[],2);
end
end
end
rmax=nij;
c=cat(3,zeros(nspi+1,nspj+1),lc);
clast=squeeze(c(end,end,:));
return
]]>
spkdallq_final
<![CDATA[
function dists=spkdallq_final(qlist,sa,sb,clast,rmax)
% dists=spkdallq_final(qlist,sa,sb,clast,rmax) does the
% final portion of parallel DP spike tima algorithm, extended to simultaneous
% calculation for all values of q.
%
% dists: a column vector of the distances
% qlist: a column vector of values of q. qklist(:,1): q-values
% sa, sb: spike times on the two spike trains
% clast, rmax: calculated by spkdallq_recur
%
% See spkdallqk.doc.
%
% Copyright (c) by Jonathan Victor.
%
% See also SPKDALLQ_RECUR, SPKDALLQ_DIST, SPKD, SPKDALLQ.
%
%
na=length(sa);
nb=length(sb);
nq=length(qlist);
qlist=reshape(qlist,[nq 1]);
%
%make column vector of na+nb-2r, indicating costs of deletions and insertions
%
nanbrs=(na+nb)*ones(1+rmax,1)-2*[0:rmax]';
%
% find the best strategy (all r (matched links) and all s (mismatched links)
%
clast=reshape(clast,[1+rmax 1]);
for iq=1:nq
posscosts=qlist(iq,1)*clast+nanbrs;
dists(iq,1)=min(min(posscosts));
end
return
]]>
spkdallqk_dist
<![CDATA[
function [dists,clast,rmax]=spkdallq_dist(qlist,sa,sb)
% function [dists,clast,rmax]=spkdallq_dist(qlist,sa,sb) does the
% recursion and final portion of the parallel DP single unit algorithm, extended to simultaneous
% calculation for all values of q.
%
% qlist: a column vector of values of q. qklist(:,1): q-values.
% sa, sb: spike times on the two spike trains
% clast, rmax: calculated by spkdallq_recur.
%
% See spkdallqk.doc.
%
% Copyright (c) by Jonathan Victor.
%
% See also SPKDALLQ_RECUR, SPKDALLQ_DIST.
%
%do the recursion
[clast,rmax,c]=spkdallq_recur(sa,sb);
%do the final stage
dists=spkdallq_final(qlist,sa,sb,clast,rmax);
return
]]>