# Matlab Code for Spike Time Distances, Parallel in q

#### Jonathan Victor: jdvicto@med.cornell.edu

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

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

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

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