c/Mex Code for Spike Time Distances Between Spike
Trains
main cost-based metrics page
algorithm page for cost-based metrics
c/MEX Code for Spike Time Metric
<![CDATA[
function d=spkdl(stime,sstart,send,cost)
%SPKDL. Calculates the "spike time" distance matrix at fixed cost. (MEX file version.)
% D=SPKDL(STIME,SSTART,SEND,COST). Calculates the full distance matrix for the "spike time" metric
% (Victor & Purpura 1996) for a single cost. Compiled as a MEX file (SPKDL.DLL) on Windows NT, but
% compilation seems to work on Windows 98 and Windows 2000 platforms.
%
% STIME list of spike times (row vector)
% SSTART list of trial start indices (row vector, length: # stimulus XMPsentations)
% SEND list of trial end indices (row vector, length: same as SSTART); if a trial has no
% spikes in it, the SEND value should be one less than the corresponding SSTART value
% COST cost per unit time to move a spike
%
% D the distance , as a row vector (length: #XMPsentations^2)
%
% Copyright (c) 1999 by Daniel Reich and Jonathan Victor.
% Written as a MEX file by Daniel Reich from FORTRAN code by Jonathan Victor.
/*=================================================================
*
* SPKDL.C .MEX file corresponding to SPKD.M
* CALCULATES DISTANCE BETWEEN TWO SPIKE TRAINS
* IN THE SPIKE TIME METRIC BY A CONTINUUM
* MODIFICATION OF THE SELLERS ALGORITHM
*
* INPUT IS ACTUALLY A SERIES OF SPIKE TRAINS WITH INDICES
*
* INPUT VARIABLES:
* STIME: LIST OF SPIKE TIMES (1...n)
* SSTART: LIST OF TRIAL START INDICES
* SEND: LIST OF TRIAL END INDICES
* COST: THE COST PER UNIT TIME TO MOVE A SPIKE
*
* D: THE DISTANCE
*
* MAXCOST is set to 100000
*
* Copyright (c) 1999 by Daniel Reich and Jonathan Victor
* All rights reserved.
*=================================================================*/
/* $Revision: 1.0 $ */
#include "mex.h"
/* Input Arguments */
#define STIME prhs[0]
#define SSTART prhs[1]
#define SEND prhs[2]
#define COST prhs[3]
/* Output Arguments */
#define D plhs[0]
#if !defined(MAX)
#define MAX(A, B) ((A) > (B) ? (A) : (B))
#endif
#if !defined(MIN)
#define MIN(A, B) ((A) < (B) ? (A) : (B))
#endif
#define PI 3.14159265
#define MAXCOST 10000
#define EPS 0.000001
#define NR_END 1
#define FREE_ARG char*
void nrerror(char error_text[])
/* Numerical Recipes standard error handler */
{
fprintf(stderr,"Numerical Recipes run-time error...\n");
fprintf(stderr,"%s\n",error_text);
fprintf(stderr,"...now exiting to system...\n");
exit(1);
}
double *dvector(long nl, long nh)
/* allocate a double vector with subscript range v[nl..nh] */
{
double *v;
v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double)));
if (!v) nrerror("allocation failure in dvector()");
return v-nl+NR_END;
}
double **dmatrix(long nrl, long nrh, long ncl, long nch)
/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
{
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
double **m;
/* allocate pointers to rows */
m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));
if (!m) nrerror("allocation failure 1 in matrix()");
m += NR_END;
m -= nrl;
/* allocate rows and set pointers to them */
m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;
for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
void free_dvector(double *v, long nl, long nh)
/* free a double vector allocated with dvector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch)
/* free a double matrix allocated by dmatrix() */
{
free((FREE_ARG) (m[nrl]+ncl-NR_END));
free((FREE_ARG) (m+nrl-NR_END));
}
double dabs(double x)
{
if (x<0) return -x;
else return x;
}
static void getdist(
double *d,
unsigned int nspi,
double tli[],
unsigned int nspj,
double tlj[],
double cost
)
{
double **scr;
unsigned int i,j;
if (cost<EPS) {
*d=dabs((double)nspi-(double)nspj);
/*printf("\n%d %d %f %f",nspi,nspj,dabs((double)nspi-(double)nspj),*d);*/
}
else if (cost>=MAXCOST)
*d=nspi+nspj;
else {
/* INITIALIZE MARGINS WITH COST OF ADDING A SPIKE */
scr=dmatrix(0,nspi,0,nspj);
for (i=0; i<=nspi; i++)
scr[i][0]=i;
for (j=0; j<=nspj; j++)
scr[0][j]=j;
if (nspi!=0 && nspj!=0) {
/* THE HEART OF THE ALGORITHM */
for (i=1; i<=nspi; i++)
for (j=1; j<=nspj; j++)
scr[i][j]=MIN(MIN(scr[i-1][j]+1,scr[i][j-1]+1),scr[i-1][j-1]+cost*dabs(tli[i-1]-tlj[j-1]));
}
*d=scr[nspi][nspj];
free_dmatrix(scr,0,nspi,0,nspj);
}
return;
}
static void getdistl(
double *d,
unsigned int nstart,
double stime[],
double sstart[],
double send[],
double cost
)
{
double *tli, *tlj;
double dist;
unsigned int i,j,itli,itlj;
unsigned int nspi,nspj;
double junk;
for (i=0; i<nstart; i++) {
nspi=(int)(send[i]-sstart[i]+1);
if (nspi>0) {
tli=dvector(0,nspi-1);
for (itli=0; itli<nspi; itli++) {
tli[itli]=stime[itli+(int)sstart[i]-1];
/*printf("\n%d %f",itli,tli[itli]);*/
}
}
for (j=i+1; j<nstart; j++) {
nspj=(int)(send[j]-sstart[j]+1);
/*printf("\nnspi=%d\tnspj=%d",nspi,nspj);*/
if (nspi>0 && nspj>0) {
tlj=dvector(0,nspj-1);
for (itlj=0; itlj<nspj; itlj++) {
tlj[itlj]=stime[itlj+(int)sstart[j]-1];
/*printf("\n%d %f",itlj,tlj[itlj]);*/
}
getdist(&dist,nspi,tli,nspj,tlj,cost);
/*printf("\t%f",dist);*/
free_dvector(tlj,0,nspj-1);
}
if (nspi==0 && nspj>0) dist=nspj;
else if (nspj==0 && nspi>0) dist=nspi;
else if (nspj==0 && nspi==0) dist=0;
/*printf("\n%d %d %d %f",i,j,i*nstart+j,dist);*/
d[i*nstart+j]=d[j*nstart+i]=dist;
}
if (nspi>0) free_dvector(tli,0,nspi-1);
}
return;
}
void mexFunction( int nlhs, mxArray *plhs[],
int nrhs, const mxArray*prhs[] )
{
double *stime, cost;
double *sstart, *send;
double *d;
unsigned int nspikes,nstart,nend,junk;
/* Check for proper number of arguments */
if (nrhs != 4) {
mexErrMsgTxt("Four input arguments required.");
} else if (nlhs > 1) {
mexErrMsgTxt("Too many output arguments.");
}
/* Check the dimensions of STIME. STIME can be n X 1 or 1 X n. */
nspikes = mxGetM(STIME);
junk = mxGetN(STIME);
if (!mxIsDouble(STIME) || (MIN(nspikes,junk) > 1)) {
mexErrMsgTxt("SPKDL requires that STIME be a row or column vector.");
}
nspikes=MAX(nspikes,junk);
nstart = mxGetM(SSTART);
junk = mxGetN(SSTART);
nstart = MAX(nstart,junk);
nend = mxGetM(SEND);
junk = mxGetN(SEND);
nend= MAX(nend,junk);
if (nstart != nend) {
mexErrMsgTxt("SPKDL requires that SSTART and SEND have the same dimensions.");
}
cost = mxGetM(COST);
junk = mxGetN(COST);
if (!mxIsDouble(COST) || (MAX(cost,junk) != 1) || (MIN(cost,junk) != 1)) {
mexErrMsgTxt("SPKD requires that COST be a scalar.");
}
/* Create a matrix for the return argument */
D = mxCreateDoubleMatrix(nstart*nstart, 1, mxREAL);
/* Assign pointers to the various parameters */
d = mxGetPr(D);
stime = mxGetPr(STIME);
sstart = mxGetPr(SSTART);
send = mxGetPr(SEND);
cost = mxGetScalar(COST);
/* Do the actual computations in a subroutine */
getdistl(d,nstart,stime,sstart,send,cost);
return;
}
]]>