/* closefuns.h Function definitions to be #included in closepair.c several times with different values of macros. Macros used: CLOSEFUN name of function for 'closepairs' CROSSFUN name of function for 'crosspairs' DIST if defined, also return d COORDS if defined, also return xi, yi, xj, yj, dx, dy THRESH if defined, also return 1(d < s) ZCOORD if defined, coordinates are 3-dimensional SINGLE if defined, capture only i < j $Revision: 1.9 $ $Date: 2015/12/30 04:01:51 $ */ #ifdef ZCOORD #define SPACEDIM 3 #else #define SPACEDIM 2 #endif SEXP CLOSEFUN(SEXP xx, SEXP yy, #ifdef ZCOORD SEXP zz, #endif SEXP rr, #ifdef THRESH SEXP ss, #endif SEXP nguess) { double *x, *y; double xi, yi, rmax, r2max, rmaxplus, dx, dy, d2; #ifdef ZCOORD double *z; double zi, dz; #endif int n, k, kmax, kmaxold, maxchunk, i, j, m; /* local storage */ int *iout, *jout; /* R objects in return value */ SEXP Out, iOut, jOut; /* external storage pointers */ int *iOutP, *jOutP; #ifdef COORDS double *xiout, *yiout, *xjout, *yjout, *dxout, *dyout; SEXP xiOut, yiOut, xjOut, yjOut, dxOut, dyOut; double *xiOutP, *yiOutP, *xjOutP, *yjOutP, *dxOutP, *dyOutP; #ifdef ZCOORD double *ziout, *zjout, *dzout; SEXP ziOut, zjOut, dzOut; double *ziOutP, *zjOutP, *dzOutP; #endif #endif #ifdef DIST double *dout; SEXP dOut; double *dOutP; #endif #ifdef THRESH double s, s2; int *tout; SEXP tOut; int *tOutP; #endif /* protect R objects from garbage collector */ PROTECT(xx = AS_NUMERIC(xx)); PROTECT(yy = AS_NUMERIC(yy)); #ifdef ZCOORD PROTECT(zz = AS_NUMERIC(zz)); #endif PROTECT(rr = AS_NUMERIC(rr)); PROTECT(nguess = AS_INTEGER(nguess)); #ifdef THRESH PROTECT(ss = AS_NUMERIC(ss)); #define NINPUTS (3+SPACEDIM) #else #define NINPUTS (2+SPACEDIM) #endif /* Translate arguments from R to C */ x = NUMERIC_POINTER(xx); y = NUMERIC_POINTER(yy); #ifdef ZCOORD z = NUMERIC_POINTER(zz); #endif n = LENGTH(xx); rmax = *(NUMERIC_POINTER(rr)); kmax = *(INTEGER_POINTER(nguess)); r2max = rmax * rmax; rmaxplus = rmax + rmax/16.0; #ifdef THRESH s = *(NUMERIC_POINTER(ss)); s2 = s * s; #endif k = 0; /* k is the next available storage location and also the current length of the list */ if(n > 0 && kmax > 0) { /* allocate space */ iout = (int *) R_alloc(kmax, sizeof(int)); jout = (int *) R_alloc(kmax, sizeof(int)); #ifdef COORDS xiout = (double *) R_alloc(kmax, sizeof(double)); yiout = (double *) R_alloc(kmax, sizeof(double)); xjout = (double *) R_alloc(kmax, sizeof(double)); yjout = (double *) R_alloc(kmax, sizeof(double)); dxout = (double *) R_alloc(kmax, sizeof(double)); dyout = (double *) R_alloc(kmax, sizeof(double)); #ifdef ZCOORD ziout = (double *) R_alloc(kmax, sizeof(double)); zjout = (double *) R_alloc(kmax, sizeof(double)); dzout = (double *) R_alloc(kmax, sizeof(double)); #endif #endif #ifdef DIST dout = (double *) R_alloc(kmax, sizeof(double)); #endif #ifdef THRESH tout = (int *) R_alloc(kmax, sizeof(int)); #endif /* loop in chunks of 2^16 */ i = 0; maxchunk = 0; while(i < n) { R_CheckUserInterrupt(); maxchunk += 65536; if(maxchunk > n) maxchunk = n; for(; i < maxchunk; i++) { xi = x[i]; yi = y[i]; #ifdef ZCOORD zi = z[i]; #endif #ifndef SINGLE if(i > 0) { /* scan backward */ for(j = i - 1; j >= 0; j--) { dx = x[j] - xi; if(dx < -rmaxplus) break; dy = y[j] - yi; d2 = dx * dx + dy * dy; #ifdef ZCOORD if(d2 <= r2max) { dz = z[j] - zi; d2 = d2 + dz * dz; #endif if(d2 <= r2max) { /* add this (i, j) pair to output */ if(k >= kmax) { /* overflow; allocate more space */ kmaxold = kmax; kmax = 2 * kmax; iout = intRealloc(iout, kmaxold, kmax); jout = intRealloc(jout, kmaxold, kmax); #ifdef COORDS xiout = dblRealloc(xiout, kmaxold, kmax); yiout = dblRealloc(yiout, kmaxold, kmax); xjout = dblRealloc(xjout, kmaxold, kmax); yjout = dblRealloc(yjout, kmaxold, kmax); dxout = dblRealloc(dxout, kmaxold, kmax); dyout = dblRealloc(dyout, kmaxold, kmax); #ifdef ZCOORD ziout = dblRealloc(ziout, kmaxold, kmax); zjout = dblRealloc(zjout, kmaxold, kmax); dzout = dblRealloc(dzout, kmaxold, kmax); #endif #endif #ifdef DIST dout = dblRealloc(dout, kmaxold, kmax); #endif #ifdef THRESH tout = intRealloc(tout, kmaxold, kmax); #endif } jout[k] = j + 1; /* R indexing */ iout[k] = i + 1; #ifdef COORDS xiout[k] = xi; yiout[k] = yi; xjout[k] = x[j]; yjout[k] = y[j]; dxout[k] = dx; dyout[k] = dy; #ifdef ZCOORD ziout[k] = zi; zjout[k] = z[j]; dzout[k] = dz; #endif #endif #ifdef DIST dout[k] = sqrt(d2); #endif #ifdef THRESH tout[k] = (d2 <= s2) ? 1 : 0; #endif ++k; } #ifdef ZCOORD } #endif } } #endif if(i + 1 < n) { /* scan forward */ for(j = i + 1; j < n; j++) { dx = x[j] - xi; if(dx > rmaxplus) break; dy = y[j] - yi; d2 = dx * dx + dy * dy; #ifdef ZCOORD if(d2 <= r2max) { dz = z[j] - zi; d2 = d2 + dz * dz; #endif if(d2 <= r2max) { /* add this (i, j) pair to output */ if(k >= kmax) { /* overflow; allocate more space */ kmaxold = kmax; kmax = 2 * kmax; iout = intRealloc(iout, kmaxold, kmax); jout = intRealloc(jout, kmaxold, kmax); #ifdef COORDS xiout = dblRealloc(xiout, kmaxold, kmax); yiout = dblRealloc(yiout, kmaxold, kmax); xjout = dblRealloc(xjout, kmaxold, kmax); yjout = dblRealloc(yjout, kmaxold, kmax); dxout = dblRealloc(dxout, kmaxold, kmax); dyout = dblRealloc(dyout, kmaxold, kmax); #ifdef ZCOORD ziout = dblRealloc(ziout, kmaxold, kmax); zjout = dblRealloc(zjout, kmaxold, kmax); dzout = dblRealloc(dzout, kmaxold, kmax); #endif #endif #ifdef DIST dout = dblRealloc(dout, kmaxold, kmax); #endif #ifdef THRESH tout = intRealloc(tout, kmaxold, kmax); #endif } jout[k] = j + 1; /* R indexing */ iout[k] = i + 1; #ifdef COORDS xiout[k] = xi; yiout[k] = yi; xjout[k] = x[j]; yjout[k] = y[j]; dxout[k] = dx; dyout[k] = dy; #ifdef ZCOORD ziout[k] = zi; zjout[k] = z[j]; dzout[k] = dz; #endif #endif #ifdef DIST dout[k] = sqrt(d2); #endif #ifdef THRESH tout[k] = (d2 <= s2) ? 1 : 0; #endif ++k; } #ifdef ZCOORD } #endif } } /* end of i loop */ } } } /* return a list of vectors */ PROTECT(iOut = NEW_INTEGER(k)); PROTECT(jOut = NEW_INTEGER(k)); #ifdef COORDS PROTECT(xiOut = NEW_NUMERIC(k)); PROTECT(yiOut = NEW_NUMERIC(k)); PROTECT(xjOut = NEW_NUMERIC(k)); PROTECT(yjOut = NEW_NUMERIC(k)); PROTECT(dxOut = NEW_NUMERIC(k)); PROTECT(dyOut = NEW_NUMERIC(k)); #ifdef ZCOORD PROTECT(ziOut = NEW_NUMERIC(k)); PROTECT(zjOut = NEW_NUMERIC(k)); PROTECT(dzOut = NEW_NUMERIC(k)); #endif #endif #ifdef DIST PROTECT(dOut = NEW_NUMERIC(k)); #endif #ifdef THRESH PROTECT(tOut = NEW_INTEGER(k)); #endif if(k > 0) { iOutP = INTEGER_POINTER(iOut); jOutP = INTEGER_POINTER(jOut); #ifdef COORDS xiOutP = NUMERIC_POINTER(xiOut); yiOutP = NUMERIC_POINTER(yiOut); xjOutP = NUMERIC_POINTER(xjOut); yjOutP = NUMERIC_POINTER(yjOut); dxOutP = NUMERIC_POINTER(dxOut); dyOutP = NUMERIC_POINTER(dyOut); #ifdef ZCOORD ziOutP = NUMERIC_POINTER(ziOut); zjOutP = NUMERIC_POINTER(zjOut); dzOutP = NUMERIC_POINTER(dzOut); #endif #endif #ifdef DIST dOutP = NUMERIC_POINTER(dOut); #endif #ifdef THRESH tOutP = INTEGER_POINTER(tOut); #endif for(m = 0; m < k; m++) { iOutP[m] = iout[m]; jOutP[m] = jout[m]; #ifdef COORDS xiOutP[m] = xiout[m]; yiOutP[m] = yiout[m]; xjOutP[m] = xjout[m]; yjOutP[m] = yjout[m]; dxOutP[m] = dxout[m]; dyOutP[m] = dyout[m]; #ifdef ZCOORD ziOutP[m] = ziout[m]; zjOutP[m] = zjout[m]; dzOutP[m] = dzout[m]; #endif #endif #ifdef DIST dOutP[m] = dout[m]; #endif #ifdef THRESH tOutP[m] = tout[m]; #endif } } #define HEAD 2 #ifdef THRESH #define NECK 1 #else #define NECK 0 #endif #ifdef COORDS #define MIDDLE (3*SPACEDIM) #else #define MIDDLE 0 #endif #ifdef DIST #define TAIL 1 #else #define TAIL 0 #endif PROTECT(Out = NEW_LIST(HEAD+NECK+MIDDLE+TAIL)); SET_VECTOR_ELT(Out, 0, iOut); SET_VECTOR_ELT(Out, 1, jOut); #ifdef THRESH SET_VECTOR_ELT(Out, HEAD, tOut); #endif #ifdef COORDS #ifdef ZCOORD SET_VECTOR_ELT(Out, HEAD+NECK, xiOut); SET_VECTOR_ELT(Out, HEAD+NECK+1, yiOut); SET_VECTOR_ELT(Out, HEAD+NECK+2, ziOut); SET_VECTOR_ELT(Out, HEAD+NECK+3, xjOut); SET_VECTOR_ELT(Out, HEAD+NECK+4, yjOut); SET_VECTOR_ELT(Out, HEAD+NECK+5, zjOut); SET_VECTOR_ELT(Out, HEAD+NECK+6, dxOut); SET_VECTOR_ELT(Out, HEAD+NECK+7, dyOut); SET_VECTOR_ELT(Out, HEAD+NECK+8, dzOut); #else SET_VECTOR_ELT(Out, HEAD+NECK, xiOut); SET_VECTOR_ELT(Out, HEAD+NECK+1, yiOut); SET_VECTOR_ELT(Out, HEAD+NECK+2, xjOut); SET_VECTOR_ELT(Out, HEAD+NECK+3, yjOut); SET_VECTOR_ELT(Out, HEAD+NECK+4, dxOut); SET_VECTOR_ELT(Out, HEAD+NECK+5, dyOut); #endif #endif #ifdef DIST SET_VECTOR_ELT(Out, HEAD+NECK+MIDDLE, dOut); #endif UNPROTECT(NINPUTS+1+HEAD+NECK+MIDDLE+TAIL); /* 1 is for 'Out' itself */ return(Out); } #undef NINPUTS #undef HEAD #undef NECK #undef MIDDLE #undef TAIL /* ........................................................ */ SEXP CROSSFUN(SEXP xx1, SEXP yy1, #ifdef ZCOORD SEXP zz1, #endif SEXP xx2, SEXP yy2, #ifdef ZCOORD SEXP zz2, #endif SEXP rr, #ifdef THRESH SEXP ss, #endif SEXP nguess) { /* input vectors */ double *x1, *y1, *x2, *y2; #ifdef ZCOORD double *z1, *z2; #endif /* lengths */ int n1, n2, nout, noutmax, noutmaxold, maxchunk; /* distance parameter */ double rmax, r2max, rmaxplus; /* indices */ int i, j, jleft, m; /* temporary values */ double x1i, y1i, xleft, dx, dy, dx2, d2; #ifdef ZCOORD double z1i, dz; #endif /* local storage */ int *iout, *jout; /* R objects in return value */ SEXP Out, iOut, jOut; /* external storage pointers */ int *iOutP, *jOutP; #ifdef COORDS SEXP xiOut, yiOut, xjOut, yjOut, dxOut, dyOut; double *xiOutP, *yiOutP, *xjOutP, *yjOutP, *dxOutP, *dyOutP; double *xiout, *yiout, *xjout, *yjout, *dxout, *dyout; #ifdef ZCOORD SEXP ziOut, zjOut, dzOut; double *ziOutP, *zjOutP, *dzOutP; double *ziout, *zjout, *dzout; #endif #endif #ifdef DIST SEXP dOut; double *dOutP; double *dout; #endif #ifdef THRESH double s, s2; int *tout; SEXP tOut; int *tOutP; #endif /* protect R objects from garbage collector */ PROTECT(xx1 = AS_NUMERIC(xx1)); PROTECT(yy1 = AS_NUMERIC(yy1)); PROTECT(xx2 = AS_NUMERIC(xx2)); PROTECT(yy2 = AS_NUMERIC(yy2)); #ifdef ZCOORD PROTECT(zz1 = AS_NUMERIC(zz1)); PROTECT(zz2 = AS_NUMERIC(zz2)); #endif PROTECT(rr = AS_NUMERIC(rr)); PROTECT(nguess = AS_INTEGER(nguess)); #ifdef THRESH PROTECT(ss = AS_NUMERIC(ss)); #define NINPUTS (2*SPACEDIM + 3) #else #define NINPUTS (2*SPACEDIM + 2) #endif /* Translate arguments from R to C */ x1 = NUMERIC_POINTER(xx1); y1 = NUMERIC_POINTER(yy1); x2 = NUMERIC_POINTER(xx2); y2 = NUMERIC_POINTER(yy2); #ifdef ZCOORD z1 = NUMERIC_POINTER(zz1); z2 = NUMERIC_POINTER(zz2); #endif n1 = LENGTH(xx1); n2 = LENGTH(xx2); rmax = *(NUMERIC_POINTER(rr)); noutmax = *(INTEGER_POINTER(nguess)); r2max = rmax * rmax; rmaxplus = rmax + rmax/16.0; #ifdef THRESH s = *(NUMERIC_POINTER(ss)); s2 = s * s; #endif nout = 0; /* nout is the next available storage location and also the current length of the list */ if(n1 > 0 && n2 > 0 && noutmax > 0) { /* allocate space */ iout = (int *) R_alloc(noutmax, sizeof(int)); jout = (int *) R_alloc(noutmax, sizeof(int)); #ifdef COORDS xiout = (double *) R_alloc(noutmax, sizeof(double)); yiout = (double *) R_alloc(noutmax, sizeof(double)); xjout = (double *) R_alloc(noutmax, sizeof(double)); yjout = (double *) R_alloc(noutmax, sizeof(double)); dxout = (double *) R_alloc(noutmax, sizeof(double)); dyout = (double *) R_alloc(noutmax, sizeof(double)); #ifdef ZCOORD ziout = (double *) R_alloc(noutmax, sizeof(double)); zjout = (double *) R_alloc(noutmax, sizeof(double)); dzout = (double *) R_alloc(noutmax, sizeof(double)); #endif #endif #ifdef DIST dout = (double *) R_alloc(noutmax, sizeof(double)); #endif #ifdef THRESH tout = (int *) R_alloc(noutmax, sizeof(int)); #endif jleft = 0; i = 0; maxchunk = 0; while(i < n1) { R_CheckUserInterrupt(); maxchunk += 65536; if(maxchunk > n1) maxchunk = n1; for( ; i < maxchunk; i++) { x1i = x1[i]; y1i = y1[i]; #ifdef ZCOORD z1i = z1[i]; #endif /* adjust starting point jleft */ xleft = x1i - rmaxplus; while((x2[jleft] < xleft) && (jleft+1 < n2)) ++jleft; /* process from j = jleft until dx > rmax + epsilon */ for(j=jleft; j < n2; j++) { /* squared interpoint distance */ dx = x2[j] - x1i; if(dx > rmaxplus) break; dx2 = dx * dx; dy = y2[j] - y1i; d2 = dx2 + dy * dy; #ifdef ZCOORD if(d2 <= r2max) { dz = z2[j] - z1i; d2 = d2 + dz * dz; #endif if(d2 <= r2max) { /* add this (i, j) pair to output */ if(nout >= noutmax) { /* overflow; allocate more space */ noutmaxold = noutmax; noutmax = 2 * noutmax; iout = intRealloc(iout, noutmaxold, noutmax); jout = intRealloc(jout, noutmaxold, noutmax); #ifdef COORDS xiout = dblRealloc(xiout, noutmaxold, noutmax); yiout = dblRealloc(yiout, noutmaxold, noutmax); xjout = dblRealloc(xjout, noutmaxold, noutmax); yjout = dblRealloc(yjout, noutmaxold, noutmax); dxout = dblRealloc(dxout, noutmaxold, noutmax); dyout = dblRealloc(dyout, noutmaxold, noutmax); #ifdef ZCOORD ziout = dblRealloc(ziout, noutmaxold, noutmax); zjout = dblRealloc(zjout, noutmaxold, noutmax); dzout = dblRealloc(dzout, noutmaxold, noutmax); #endif #endif #ifdef DIST dout = dblRealloc(dout, noutmaxold, noutmax); #endif #ifdef THRESH tout = intRealloc(tout, noutmaxold, noutmax); #endif } iout[nout] = i + 1; /* R indexing */ jout[nout] = j + 1; #ifdef COORDS xiout[nout] = x1i; yiout[nout] = y1i; xjout[nout] = x2[j]; yjout[nout] = y2[j]; dxout[nout] = dx; dyout[nout] = dy; #ifdef ZCOORD ziout[nout] = z1i; zjout[nout] = z2[j]; dzout[nout] = dz; #endif #endif #ifdef DIST dout[nout] = sqrt(d2); #endif #ifdef THRESH tout[nout] = (d2 <= s2) ? 1 : 0; #endif ++nout; } #ifdef ZCOORD } #endif } } } } /* return a list of vectors */ PROTECT(iOut = NEW_INTEGER(nout)); PROTECT(jOut = NEW_INTEGER(nout)); #ifdef COORDS PROTECT(xiOut = NEW_NUMERIC(nout)); PROTECT(yiOut = NEW_NUMERIC(nout)); PROTECT(xjOut = NEW_NUMERIC(nout)); PROTECT(yjOut = NEW_NUMERIC(nout)); PROTECT(dxOut = NEW_NUMERIC(nout)); PROTECT(dyOut = NEW_NUMERIC(nout)); #ifdef ZCOORD PROTECT(ziOut = NEW_NUMERIC(nout)); PROTECT(zjOut = NEW_NUMERIC(nout)); PROTECT(dzOut = NEW_NUMERIC(nout)); #endif #endif #ifdef DIST PROTECT(dOut = NEW_NUMERIC(nout)); #endif #ifdef THRESH PROTECT(tOut = NEW_INTEGER(nout)); #endif if(nout > 0) { iOutP = INTEGER_POINTER(iOut); jOutP = INTEGER_POINTER(jOut); #ifdef COORDS xiOutP = NUMERIC_POINTER(xiOut); yiOutP = NUMERIC_POINTER(yiOut); xjOutP = NUMERIC_POINTER(xjOut); yjOutP = NUMERIC_POINTER(yjOut); dxOutP = NUMERIC_POINTER(dxOut); dyOutP = NUMERIC_POINTER(dyOut); #ifdef ZCOORD ziOutP = NUMERIC_POINTER(ziOut); zjOutP = NUMERIC_POINTER(zjOut); dzOutP = NUMERIC_POINTER(dzOut); #endif #endif #ifdef DIST dOutP = NUMERIC_POINTER(dOut); #endif #ifdef THRESH tOutP = INTEGER_POINTER(tOut); #endif for(m = 0; m < nout; m++) { iOutP[m] = iout[m]; jOutP[m] = jout[m]; #ifdef COORDS xiOutP[m] = xiout[m]; yiOutP[m] = yiout[m]; xjOutP[m] = xjout[m]; yjOutP[m] = yjout[m]; dxOutP[m] = dxout[m]; dyOutP[m] = dyout[m]; #ifdef ZCOORD ziOutP[m] = ziout[m]; zjOutP[m] = zjout[m]; dzOutP[m] = dzout[m]; #endif #endif #ifdef DIST dOutP[m] = dout[m]; #endif #ifdef THRESH tOutP[m] = tout[m]; #endif } } #define HEAD 2 #ifdef THRESH #define NECK 1 #else #define NECK 0 #endif #ifdef COORDS #define MIDDLE (3*SPACEDIM) #else #define MIDDLE 0 #endif #ifdef DIST #define TAIL 1 #else #define TAIL 0 #endif PROTECT(Out = NEW_LIST(HEAD+NECK+MIDDLE+TAIL)); SET_VECTOR_ELT(Out, 0, iOut); SET_VECTOR_ELT(Out, 1, jOut); #ifdef THRESH SET_VECTOR_ELT(Out, HEAD, tOut); #endif #ifdef COORDS #ifdef ZCOORD SET_VECTOR_ELT(Out, HEAD+NECK, xiOut); SET_VECTOR_ELT(Out, HEAD+NECK+1, yiOut); SET_VECTOR_ELT(Out, HEAD+NECK+2, ziOut); SET_VECTOR_ELT(Out, HEAD+NECK+3, xjOut); SET_VECTOR_ELT(Out, HEAD+NECK+4, yjOut); SET_VECTOR_ELT(Out, HEAD+NECK+5, zjOut); SET_VECTOR_ELT(Out, HEAD+NECK+6, dxOut); SET_VECTOR_ELT(Out, HEAD+NECK+7, dyOut); SET_VECTOR_ELT(Out, HEAD+NECK+8, dzOut); #else SET_VECTOR_ELT(Out, HEAD+NECK, xiOut); SET_VECTOR_ELT(Out, HEAD+NECK+1, yiOut); SET_VECTOR_ELT(Out, HEAD+NECK+2, xjOut); SET_VECTOR_ELT(Out, HEAD+NECK+3, yjOut); SET_VECTOR_ELT(Out, HEAD+NECK+4, dxOut); SET_VECTOR_ELT(Out, HEAD+NECK+5, dyOut); #endif #endif #ifdef DIST SET_VECTOR_ELT(Out, HEAD+NECK+MIDDLE, dOut); #endif UNPROTECT(NINPUTS+1+HEAD+NECK+MIDDLE+TAIL); /* 1 is for 'Out' itself */ return(Out); } #undef NINPUTS #undef HEAD #undef NECK #undef MIDDLE #undef TAIL /* ........................................................ */ /* Alternative code for CLOSEFUN, based on algorithm in CROSSFUN */ #define ALT_ALGO(NAME) ALT_PREFIX(NAME) #define ALT_PREFIX(NAME) alt ## NAME SEXP ALT_ALGO(CLOSEFUN)(SEXP xx, SEXP yy, #ifdef ZCOORD SEXP zz, #endif SEXP rr, #ifdef THRESH SEXP ss, #endif SEXP nguess) { /* input vectors */ double *x, *y; #ifdef ZCOORD double *z; #endif /* lengths */ int n, nout, noutmax, noutmaxold, maxchunk; /* distance parameter */ double rmax, r2max, rmaxplus; /* indices */ int i, j, jleft, m; /* temporary values */ double xi, yi, xleft, dx, dy, dx2, d2; #ifdef ZCOORD double zi, dz; #endif /* local storage */ int *iout, *jout; /* R objects in return value */ SEXP Out, iOut, jOut; /* external storage pointers */ int *iOutP, *jOutP; #ifdef COORDS SEXP xiOut, yiOut, xjOut, yjOut, dxOut, dyOut; double *xiOutP, *yiOutP, *xjOutP, *yjOutP, *dxOutP, *dyOutP; double *xiout, *yiout, *xjout, *yjout, *dxout, *dyout; #ifdef ZCOORD SEXP ziOut, zjOut, dzOut; double *ziOutP, *zjOutP, *dzOutP; double *ziout, *zjout, *dzout; #endif #endif #ifdef DIST SEXP dOut; double *dOutP; double *dout; #endif #ifdef THRESH double s, s2; int *tout; SEXP tOut; int *tOutP; #endif /* protect R objects from garbage collector */ PROTECT(xx = AS_NUMERIC(xx)); PROTECT(yy = AS_NUMERIC(yy)); #ifdef ZCOORD PROTECT(zz = AS_NUMERIC(zz)); #endif PROTECT(rr = AS_NUMERIC(rr)); PROTECT(nguess = AS_INTEGER(nguess)); #ifdef THRESH PROTECT(ss = AS_NUMERIC(ss)); #define NINPUTS (SPACEDIM + 3) #else #define NINPUTS (SPACEDIM + 2) #endif /* Translate arguments from R to C */ x = NUMERIC_POINTER(xx); y = NUMERIC_POINTER(yy); #ifdef ZCOORD z = NUMERIC_POINTER(zz); #endif n = LENGTH(xx); rmax = *(NUMERIC_POINTER(rr)); noutmax = *(INTEGER_POINTER(nguess)); r2max = rmax * rmax; rmaxplus = rmax + rmax/16.0; #ifdef THRESH s = *(NUMERIC_POINTER(ss)); s2 = s * s; #endif nout = 0; /* nout is the next available storage location and also the current length of the list */ if(n > 0 && noutmax > 0) { /* allocate space */ iout = (int *) R_alloc(noutmax, sizeof(int)); jout = (int *) R_alloc(noutmax, sizeof(int)); #ifdef COORDS xiout = (double *) R_alloc(noutmax, sizeof(double)); yiout = (double *) R_alloc(noutmax, sizeof(double)); xjout = (double *) R_alloc(noutmax, sizeof(double)); yjout = (double *) R_alloc(noutmax, sizeof(double)); dxout = (double *) R_alloc(noutmax, sizeof(double)); dyout = (double *) R_alloc(noutmax, sizeof(double)); #ifdef ZCOORD ziout = (double *) R_alloc(noutmax, sizeof(double)); zjout = (double *) R_alloc(noutmax, sizeof(double)); dzout = (double *) R_alloc(noutmax, sizeof(double)); #endif #endif #ifdef DIST dout = (double *) R_alloc(noutmax, sizeof(double)); #endif #ifdef THRESH tout = (int *) R_alloc(noutmax, sizeof(int)); #endif jleft = 0; i = 0; maxchunk = 0; while(i < n) { R_CheckUserInterrupt(); maxchunk += 65536; if(maxchunk > n) maxchunk = n; for( ; i < maxchunk; i++) { xi = x[i]; yi = y[i]; #ifdef ZCOORD zi = z[i]; #endif /* adjust starting point jleft */ xleft = xi - rmaxplus; while((x[jleft] < xleft) && (jleft+1 < n)) ++jleft; /* process from j = jleft until dx > rmax + epsilon */ for(j=jleft; j < n; j++) { /* squared interpoint distance */ dx = x[j] - xi; if(dx > rmaxplus) break; dx2 = dx * dx; dy = y[j] - yi; d2 = dx2 + dy * dy; #ifdef ZCOORD if(d2 <= r2max) { dz = z[j] - zi; d2 = d2 + dz * dz; #endif if(d2 <= r2max) { /* add this (i, j) pair to output */ if(nout >= noutmax) { /* overflow; allocate more space */ noutmaxold = noutmax; noutmax = 2 * noutmax; iout = intRealloc(iout, noutmaxold, noutmax); jout = intRealloc(jout, noutmaxold, noutmax); #ifdef COORDS xiout = dblRealloc(xiout, noutmaxold, noutmax); yiout = dblRealloc(yiout, noutmaxold, noutmax); xjout = dblRealloc(xjout, noutmaxold, noutmax); yjout = dblRealloc(yjout, noutmaxold, noutmax); dxout = dblRealloc(dxout, noutmaxold, noutmax); dyout = dblRealloc(dyout, noutmaxold, noutmax); #ifdef ZCOORD ziout = dblRealloc(ziout, noutmaxold, noutmax); zjout = dblRealloc(zjout, noutmaxold, noutmax); dzout = dblRealloc(dzout, noutmaxold, noutmax); #endif #endif #ifdef DIST dout = dblRealloc(dout, noutmaxold, noutmax); #endif #ifdef THRESH tout = intRealloc(tout, noutmaxold, noutmax); #endif } iout[nout] = i + 1; /* R indexing */ jout[nout] = j + 1; #ifdef COORDS xiout[nout] = xi; yiout[nout] = yi; xjout[nout] = x[j]; yjout[nout] = y[j]; dxout[nout] = dx; dyout[nout] = dy; #ifdef ZCOORD ziout[nout] = zi; zjout[nout] = z[j]; dzout[nout] = dz; #endif #endif #ifdef DIST dout[nout] = sqrt(d2); #endif #ifdef THRESH tout[nout] = (d2 <= s2) ? 1 : 0; #endif ++nout; } #ifdef ZCOORD } #endif } } } } /* return a list of vectors */ PROTECT(iOut = NEW_INTEGER(nout)); PROTECT(jOut = NEW_INTEGER(nout)); #ifdef COORDS PROTECT(xiOut = NEW_NUMERIC(nout)); PROTECT(yiOut = NEW_NUMERIC(nout)); PROTECT(xjOut = NEW_NUMERIC(nout)); PROTECT(yjOut = NEW_NUMERIC(nout)); PROTECT(dxOut = NEW_NUMERIC(nout)); PROTECT(dyOut = NEW_NUMERIC(nout)); #ifdef ZCOORD PROTECT(ziOut = NEW_NUMERIC(nout)); PROTECT(zjOut = NEW_NUMERIC(nout)); PROTECT(dzOut = NEW_NUMERIC(nout)); #endif #endif #ifdef DIST PROTECT(dOut = NEW_NUMERIC(nout)); #endif #ifdef THRESH PROTECT(tOut = NEW_INTEGER(nout)); #endif if(nout > 0) { iOutP = INTEGER_POINTER(iOut); jOutP = INTEGER_POINTER(jOut); #ifdef COORDS xiOutP = NUMERIC_POINTER(xiOut); yiOutP = NUMERIC_POINTER(yiOut); xjOutP = NUMERIC_POINTER(xjOut); yjOutP = NUMERIC_POINTER(yjOut); dxOutP = NUMERIC_POINTER(dxOut); dyOutP = NUMERIC_POINTER(dyOut); #ifdef ZCOORD ziOutP = NUMERIC_POINTER(ziOut); zjOutP = NUMERIC_POINTER(zjOut); dzOutP = NUMERIC_POINTER(dzOut); #endif #endif #ifdef DIST dOutP = NUMERIC_POINTER(dOut); #endif #ifdef THRESH tOutP = INTEGER_POINTER(tOut); #endif for(m = 0; m < nout; m++) { iOutP[m] = iout[m]; jOutP[m] = jout[m]; #ifdef COORDS xiOutP[m] = xiout[m]; yiOutP[m] = yiout[m]; xjOutP[m] = xjout[m]; yjOutP[m] = yjout[m]; dxOutP[m] = dxout[m]; dyOutP[m] = dyout[m]; #ifdef ZCOORD ziOutP[m] = ziout[m]; zjOutP[m] = zjout[m]; dzOutP[m] = dzout[m]; #endif #endif #ifdef DIST dOutP[m] = dout[m]; #endif #ifdef THRESH tOutP[m] = tout[m]; #endif } } #define HEAD 2 #ifdef THRESH #define NECK 1 #else #define NECK 0 #endif #ifdef COORDS #define MIDDLE (3*SPACEDIM) #else #define MIDDLE 0 #endif #ifdef DIST #define TAIL 1 #else #define TAIL 0 #endif PROTECT(Out = NEW_LIST(HEAD+NECK+MIDDLE+TAIL)); SET_VECTOR_ELT(Out, 0, iOut); SET_VECTOR_ELT(Out, 1, jOut); #ifdef THRESH SET_VECTOR_ELT(Out, HEAD, tOut); #endif #ifdef COORDS #ifdef ZCOORD SET_VECTOR_ELT(Out, HEAD+NECK, xiOut); SET_VECTOR_ELT(Out, HEAD+NECK+1, yiOut); SET_VECTOR_ELT(Out, HEAD+NECK+2, ziOut); SET_VECTOR_ELT(Out, HEAD+NECK+3, xjOut); SET_VECTOR_ELT(Out, HEAD+NECK+4, yjOut); SET_VECTOR_ELT(Out, HEAD+NECK+5, zjOut); SET_VECTOR_ELT(Out, HEAD+NECK+6, dxOut); SET_VECTOR_ELT(Out, HEAD+NECK+7, dyOut); SET_VECTOR_ELT(Out, HEAD+NECK+8, dzOut); #else SET_VECTOR_ELT(Out, HEAD+NECK, xiOut); SET_VECTOR_ELT(Out, HEAD+NECK+1, yiOut); SET_VECTOR_ELT(Out, HEAD+NECK+2, xjOut); SET_VECTOR_ELT(Out, HEAD+NECK+3, yjOut); SET_VECTOR_ELT(Out, HEAD+NECK+4, dxOut); SET_VECTOR_ELT(Out, HEAD+NECK+5, dyOut); #endif #endif #ifdef DIST SET_VECTOR_ELT(Out, HEAD+NECK+MIDDLE, dOut); #endif UNPROTECT(NINPUTS+1+HEAD+NECK+MIDDLE+TAIL); /* 1 is for 'Out' itself */ return(Out); } #undef NINPUTS #undef HEAD #undef NECK #undef MIDDLE #undef TAIL #undef ALT_ALGO #undef ALT_PREFIX