/* * Compute global cost matrix - companion * to the dtw R package * (c) Toni Giorgino 2007-2012 * Distributed under GPL-2 with NO WARRANTY. * * \$Id: computeCM.c 343 2013-12-11 11:04:41Z tonig \$ * */ #include #include #include #ifdef TEST_UNIT // Define R-like functions - a bad idea #include #define R_NaInt INT_MIN #define R_alloc(n,size) alloca((n)*(size)) #define error(...) { fprintf (stderr, __VA_ARGS__); exit(-1); } #else #include #endif #ifndef NAN #error "This code requires native IEEE NAN support. Possible solutions: 1) verify you are using gcc with -std=gnu99; 2) use the fallback interpreted DTW version (should happen automatically); 3) ask the author" #endif /* undo R indexing */ #define EP(ii,jj) ((jj)*nsteps+(ii)) #define EM(ii,jj) ((jj)*n+(ii)) #define CLEARCLIST { \ for(int z=0; z=nsteps) { error("Error on pattern row %d, pattern number %d out of bounds\n", i,pn[i]+1); } } /* assuming pattern ids are in ascending order */ int npats=pn[nsteps-1]+1; /* prepare a cost list per pattern */ double *clist=(double*) R_alloc(npats,sizeof(double)); /* we do not initialize the seed - the caller is supposed to do so cm[0]=lm[0]; */ /* clear the direction matrix */ for(int i=0; i=0 && jj>=0) { /* address ok? C convention */ double cc=sc[s]; if(cc==-1.0) { clist[p]=cm[EM(ii,jj)]; } else { /* we rely on NAN to propagate */ clist[p] += cc*lm[EM(ii,jj)]; } } } int minc=argmin(clist,npats); if(minc>-1) { cm[EM(i,j)]=clist[minc]; sm[EM(i,j)]=minc+1; /* convert to 1-based */ } } } /* Memory alloc'd by R_alloc is automatically freed */ } /* -------------------------------------------------- * * Wrapper for .Call, avoids several copies. Returns a list with names * "costMatrix" and "directionMatrix" */ #ifndef TEST_UNIT #include #include SEXP computeCM_Call(SEXP wm, /* logical */ SEXP lm, /* double */ SEXP cm, /* double */ SEXP dir) { /* double */ /* Get problem size */ SEXP lm_dim; PROTECT(lm_dim = GET_DIM(lm)); /* ---- 1 */ int *p_lm_dim = INTEGER_POINTER(lm_dim); /* Get pattern size */ SEXP dir_dim; PROTECT(dir_dim = GET_DIM(dir)); /* ---- 2 */ int nsteps=INTEGER_POINTER(dir_dim)[0]; /* Cost matrix (input+output 1). */ SEXP cmo; PROTECT(cmo=duplicate(cm)); /* ---- 3 */ /* Output 2: smo, INTEGER */ SEXP smo; PROTECT(smo=allocMatrix(INTSXP,p_lm_dim[0],p_lm_dim[1])); /* ---- 4 */ /* Dispatch to C */ computeCM(p_lm_dim, LOGICAL_POINTER(wm), NUMERIC_POINTER(lm), &nsteps, NUMERIC_POINTER(dir), NUMERIC_POINTER(cmo), INTEGER_POINTER(smo)); /* cmo and smo are now set. Put them in a list. From S. Blay, http://www.sfu.ca/~sblay/R-C-interface.ppt */ SEXP list_names; PROTECT(list_names = allocVector(STRSXP,2)); /* ---- 5 */ SET_STRING_ELT(list_names,0,mkChar("costMatrix")); SET_STRING_ELT(list_names,1,mkChar("directionMatrix")); // Creating a list with 2 vector elements: SEXP list; PROTECT(list = allocVector(VECSXP, 2)); /* ---- 6 */ SET_VECTOR_ELT(list, 0, cmo); SET_VECTOR_ELT(list, 1, smo); // and attaching the vector names: setAttrib(list, R_NamesSymbol, list_names); UNPROTECT(6); return list; } #endif /* Test as follows: R CMD SHLIB -d computeCM.c dyn.load("computeCM.so") lm <- matrix(nrow = 6, ncol = 6, byrow = TRUE, c( 1, 1, 2, 2, 3, 3, 1, 1, 1, 2, 2, 2, 3, 1, 2, 2, 3, 3, 3, 1, 2, 1, 1, 2, 3, 2, 1, 2, 1, 2, 3, 3, 3, 2, 1, 2 )) step.matrix <- as.matrix(structure(c(1, 1, 2, 2, 3, 3, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 2, 0, -1, 1, -1, 1, -1, 1), .Dim = c(6L, 4L), class = "stepPattern", npat = 3, norm = "N")) nsteps<-dim(step.matrix)[1] wm <- matrix(TRUE,6,6) cm <- matrix(NA,6,6) cm[1,1] <- lm[1,1]; sm <- matrix(NA,6,6) out<-.C("computeCM",NAOK=TRUE, as.integer(dim(cm)), as.logical(wm), as.double(lm), as.integer(nsteps), as.double(step.matrix), cmo=as.double(cm), smo=as.integer(sm)) cmoo<-matrix(out\$cmo,6,6) smoo<-matrix(out\$smo,6,6) storage.mode(wm) <- "logical" storage.mode(lm) <- "double" storage.mode(cm) <- "double" storage.mode(step.matrix) <- "double" out2<-.Call("computeCM_Call", wm, lm, cm, step.matrix) */ #ifdef TEST_UNIT /* -------------------------------------------------- * Unit test - for debugging */ void tm_print(int *s, double *mm, double *r); /* test main equivalent to the following mylm<-outer(1:10,1:10) globalCostNative(mylm)->myg2 */ #define TS 5000 #define TSS (TS*TS) int main(int argc,char **argv) { int ts[]={TS,TS}; int *twm; double *tlm; int tnstepsp[]={6}; double tdir[]={1, 1, 2, 2, 3, 3, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0,-1, 1,-1, 1,-1, 1}; double *tcm; int *tsm; int i,j; twm=malloc(TSS*sizeof(int)); for( i=0;i