:SPUD feature extraction algorithm implementation as appearing in:
:Data-mining of time-domain features from neural extracellular field data (2007 - in press)
:S Neymotin, DJ Uhlrich, KA Manning, WW Lytton
NEURON {
SUFFIX nothing
: BVBASE is bit vector base number (typically 0 or -1)
GLOBAL SPUD_INSTALLED, SHM_SPUD, NOV_SPUD, DEBUG_SPUD
}
PARAMETER {
SPUD_INSTALLED=0
DEBUG_SPUD=0
SHM_SPUD=4 : used in spud() for measuring sharpness
NOV_SPUD=1 : used in spud() to eliminate overlap of spikes
CREEP_SPUD=0 : used in spud() to allow left/right "creep" to local minima
}
VERBATIM
#include <stdlib.h>
#ifndef NRN_VERSION_GTEQ_8_2_0
#include <math.h>
#include <sys/time.h>
extern double* hoc_pgetarg();
extern double hoc_call_func(Symbol*, int narg);
extern FILE* hoc_obj_file_arg(int narg);
extern Object** hoc_objgetarg();
extern int vector_instance_px();
extern void* vector_arg();
extern double* vector_vec();
extern double hoc_epsilon;
extern double chkarg();
extern void set_seed();
extern int ivoc_list_count(Object*);
extern Object* ivoc_list_item(Object*, int);
extern int hoc_is_double_arg(int narg);
extern char* hoc_object_name(Object*);
char ** hoc_pgargstr();
#endif
int list_vector_px3(Object*, int, double**, void**);
int ismono1(double*,int,int);
static void hxf(void *ptr) { free(ptr); hoc_execerror("",0); }
ENDVERBATIM
:* src.spud(thresh,dlist,nqslist)
: dest.spud(src) -- default thresh=0; returns indices
: look for multiple threshold crossings to define peaks
: creates multiple parallel vectors for an NQS db
: counts peaks pointing upward -- should be all pos
: see eg decnqs.hoc:fudup() for usage
VERBATIM
//** declarations
#define UDSL 50
#define UDNQ 11
// nq=new NQS("LOC","PEAK","WIDTH","BASE","HEIGHT","START","SLICES","SHARP","INDEX","FILE","NESTED")
#define LOC nq[0] // loc of peak of spike
#define PEAK nq[1] // value at peak (absolute height)
#define WIDTH nq[2] // rt flank - lt flanks (? isn't it rt flank - LOC ?)
#define BASE nq[3] // height at base
#define HEIGHT nq[4] // peak - base
#define START nq[5] // left flank of spike?
#define SLICES nq[6] // how many slices found this spike
#define SHARP nq[7] // 2nd deriv at peak
#define INDEX nq[8] // consecutive numbering of spikes
// nq[9] // will use to fill in trace's file name at hoc level
#define NESTED nq[10] // how many bumps are nested within this one
//** procedure spud()
static double spud (void* vv) {
int i, k, m, n, nqsz, nsrc, jj[UDSL], f[UDSL], lc, dsz[UDSL], nqmax, thsz, lc2, done, dbn;
double *src, *tvec, *th, *dest[UDSL], *nq[UDNQ], *tmp, *dbx, lt, thdist;
Object *ob, *ob2;
void *vvd[UDSL], *vvth, *vnq[UDNQ];
//** read in vectors and verify sizes, etc
nsrc = vector_instance_px(vv, &src); // trace to analyze
thsz = vector_arg_px(1, &th); // vector of thresholds to check
ob = *hoc_objgetarg(2); // storage for values for each threshold
ob2 = *hoc_objgetarg(3); // list of NQS vectors for returning values
tmp = (double *)ecalloc(nsrc, sizeof(double)); // tmp is size of trace
lc = ivoc_list_count(ob);
lc2 = ivoc_list_count(ob2);
if (lc>UDSL) {printf("spud ERRF mismatch: max slice list:%d %d\n",UDSL,lc); hxf(tmp);}
if (lc2!=UDNQ){printf("spud ERRB mismatch: NQS sz is %d (%d in list)\n",UDNQ,lc2);hxf(tmp);}
if (nsrc<lc) {printf("spud ERRC mismatch: %d %d\n",lc,nsrc); hxf(tmp);} // ??
if (lc!=thsz) {printf("spud ERRA mismatch: %d %d\n",lc,thsz); hxf(tmp);}
if (!ismono1(th,thsz,-1)) {printf("spud ERRD: not mono dec %g %d\n",th[0],thsz); hxf(tmp);}
// thdist=(th[thsz-2]-th[thsz-1])/2; // NOT BEING USED: the smallest spike we will accept
for (k=0;k <lc;k++) dsz[k] =list_vector_px3(ob , k, &dest[k], &vvd[k]);
for (k=0;k<lc2;k++) {
i=list_vector_px3(ob2, k, &nq[k], &vnq[k]);
if (k==0) nqmax=i; else if (i!=nqmax) { // all NQ vecs same size
printf("spud ERRE mismatch: %d %d %d\n",k,i,nqmax); hxf(tmp); }
}
//** store crossing points and midpoints in dest[k]
// dest vectors dest[k] will store crossing points and midpoints at each th[k] slice location
// as triplets: up/max/down
for (k=0; k<lc; k++) { // iterate thru thresholds
jj[k]=f[k]=0; // jj[k] is ind into dest[k]; f[k] is flag for threshold crossings
for (i=0;i<nsrc && src[i]>th[k];i++) {} // start somewhere below this thresh th[k]
for (; i<nsrc; i++) { // iterate through trace
if (src[i]>th[k]) {
if (f[k]==0) { // ? passing thresh
if (jj[k]>=dsz[k]){printf("(%d,%d,%d) :: ",k,jj[k],dsz[k]);
hoc_execerror("Dest vec too small in spud ", 0); }
dest[k][jj[k]++] = (i-1) + (th[k]-src[i-1])/(src[i]-src[i-1]); // interpolate
f[k]=1;
tmp[k]=-1e9; dest[k][jj[k]]=-1.; // flag in tmp says that a thresh found here
}
if (f[k]==1 && src[i]>tmp[k]) { // use tmp[] even more temporarily
tmp[k]=src[i]; // pick out max
dest[k][jj[k]] = (double)i; // location of this peak
}
} else { // below thresh
if (f[k]==1) { // just passed going down
jj[k]++; // triplet will be indices of cross-up/peak/cross-down
dest[k][jj[k]++] = (i-1) + (src[i-1]-th[k])/(src[i-1]-src[i]);
f[k]=0;
}
}
}
}
//** truncate dest vectors to multiples of 3:
for (k=0;k<lc;k++) vector_resize((IvocVect*)vvd[k],(int)(floor((double)jj[k]/3.)*3.));
for (i=0; i<nsrc; i++) tmp[i]=0.; // clear temp space
//** go through all the slices to find identical peaks and save widths and locations
// tmp[] uses triplets centered around a location corresponding to a max loc in the
// original vector; the widest flanks for each are then on either side of this loc
for (k=0;k<lc;k++) { // need to go from top to bottom to widen flanks
for (i=1;i<jj[k];i+=3) { // through centers (peaks)
m=(int)dest[k][i]; // hash: place center at location
if (tmp[m-2]<0 || tmp[m-1]<0 || tmp[m+1]<0 || tmp[m+2]<0) continue; // ignore; too crowded
tmp[m]--; // count how many slices have found this peak (use negative)
tmp[m-1]=dest[k][i-1]; tmp[m+1]=dest[k][i+1]; // flanks
}
}
//** 1st (of 2) loops through tmp[] -- pick up flanks
// step through tmp[] looking for negatives which indicate the slice count and pick up
// flanks from these
// nq=new NQS("LOC","PEAK","WIDTH","BASE","HEIGHT","START","SLICES","SHARP","INDEX","FILE")
for (i=0,k=0; i<nsrc; i++) if (tmp[i]<0.) { // tmp holds neg of count of slices
if (k>=nqmax) { printf("spud ERRG OOR in NQ db: %d %d\n",k,nqmax); hxf(tmp); }
LOC[k]=(double)i; // approx location of the peak of the spike
WIDTH[k]=tmp[i+1]; // location of right side -- temp storage
START[k]=tmp[i-1]; // start of spike (left side)
SLICES[k]=-tmp[i]; // # of slices
k++;
}
nqsz=k; // k ends up as size of NQS db
if (DEBUG_SPUD && ifarg(4)) { dbn=vector_arg_px(4, &dbx); // DEBUG -- save tmp vector
if (dbn<nsrc) printf("spud ERRH: Insufficient room in debug vec (%d<%d)\n",dbn,nsrc);
else for (i=0;i<nsrc;i++) dbx[i]=tmp[i];
}
//** adjust flanks to handle nested bumps
// 3 ways to handle spike nested in a spike or elongated base:
// NB always using same slice for both L and R flanks; NOV_SPUD flag: (no-overlap)
// 0. nested spike(s) share flanks determined by shared base
// 1. nested spike(s) have individual bases, 1st and last use flanks from base
// 2. nested spike(s) have individual bases, base flanks listed separately w/out peak
// here use
// search nq vecs to compare flanks to neighboring centers
// if flanks overlap the centers on LT or RT side,
// correct them by going back to original slice loc info (in dest[])
//*** look at left side -- is this flank to left of center of another bump?
if (NOV_SPUD) for (i=0;i<nqsz;i++) { // iterate through NQS db
if ((i-1)>0 && START[i] < LOC[i-1]) { // flank is to left of prior center
if (DEBUG_SPUD) printf("LT problem %d %g %g<%g\n",i,LOC[i],START[i],LOC[i-1]);
for (m=lc-1,done=0;m>=0 && !done;m--) { // m:go from bottom (widest) to top
for (n=1;n<jj[m] && !done;n+=3) { // n:through centers
// pick out lowest slice with this peak LOC whose flank is to RT of prior peak
if (floor(dest[m][n])==LOC[i] && dest[m][n-1]>LOC[i-1]) {
// ??[i]=START[i]; // temp storage for L end of this overlap
// replace both left and right flanks at this level -- #1 above
START[i]=dest[m][n-1]; WIDTH[i]=dest[m][n+1]; done=1;
}
}
}
}
//*** now look at RT side
if ((i+1)<nqsz && WIDTH[i]>LOC[i+1]) {
if (DEBUG_SPUD) printf("RT problem %d %g %g>%g\n",i,LOC[i],WIDTH[i],LOC[i+1]);
for (m=lc-1,done=0;m>=0 && !done;m--) { // m: go from bottom to top
for (n=1;n<jj[m] && !done;n+=3) { // n: through centers
// pick out lowest slice with this peak LOC whose flank is to LT of next peak
if (floor(dest[m][n])==LOC[i] && dest[m][n+1]<LOC[i+1]) {
// ??[i]=WIDTH[i]; // end of overlap
START[i]=dest[m][n-1]; WIDTH[i]=dest[m][n+1]; done=1;
}
}
}
}
}
//make sure left and right sides of bump occur at local minima
//shouldn't creeping be before NOV_SPUD=1 overlap check???
//creeping can result only in equal borders btwn two bumps
//on one side, so it should be ok here...
if(CREEP_SPUD) for(i=0,k=0;i<nsrc;i++) if(tmp[i]<0.){
//move left side to local minima
int idx = (int)START[k];
while(idx >= 1 && src[idx] >= src[idx-1]) idx--;
START[k] = idx;
//move right side to local minima
idx = (int)WIDTH[k];
while(idx < nsrc-1 && src[idx] >= src[idx+1]) idx++;
WIDTH[k] = idx;
k++;
}
//** 2nd loop through tmp[] used to fill in the rest of NQS
// needed to split into 2 loops so that could check for overlaps and correct those
// before filling in the rest of nq
for (i=0,k=0; i<nsrc; i++) if (tmp[i]<0.) { // tmp holds neg of count of slices
// calculate a base voltage lt as interpolated value on left side
lt=src[(int)floor(START[k])]+(START[k]-floor(START[k]))*\
(src[(int)floor(START[k]+1.)]-src[(int)floor(START[k])]);
BASE[k]=lt; // base voltage
PEAK[k]=src[i]; // peak voltage
WIDTH[k] = WIDTH[k] - START[k]; // width = RT_flank-LT_flank
HEIGHT[k]=PEAK[k]-BASE[k]; // redund measure -- can eliminate
// measure of sharpness diff of 1st derivs btwn peak and SHM_SPUD dist from peak
// to get 2nd deriv would be normalized by 2*SHM_SPUD*tstep
// ??could take an ave. or max first deriv for certain distance on either side
SHARP[k]=(src[i]-src[i-(int)SHM_SPUD])-(src[i+(int)SHM_SPUD]-src[i]);
INDEX[k]=(double)k;
k++;
}
int iNumBumps = k;
//count # of other bumps nested within each bump
if(!NOV_SPUD){
for(i=0; i<iNumBumps; i++){
NESTED[i] = 0;
int j = 0;
for(;j<iNumBumps;j++){
if(i!=j && LOC[j] >= START[i] && LOC[j] <= START[i]+WIDTH[i]){
NESTED[i]+=1.0;
}
}
}
} else for(i=0;i<iNumBumps;i++) NESTED[i]=0.0;
//** finish up
for (i=0;i<lc2;i++) vector_resize((IvocVect*)vnq[i], nqsz);
if (k!=nqsz) { printf("spud ERRI INT ERR: %d %d\n",k,nqsz); hxf(tmp); }
free(tmp);
return jj[0];
}
ENDVERBATIM
:* PROCEDURE install_spud()
PROCEDURE install_spud () {
if (SPUD_INSTALLED==1) {
printf("spud.mod installed\n")
} else {
SPUD_INSTALLED=1
VERBATIM {
install_vector_method("spud", spud);
}
ENDVERBATIM
}
}