diff spp/src/cdensum.c @ 6:ce08b0efa3fd draft

Uploaded

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/spp/src/cdensum.c	Tue Nov 27 16:11:40 2012 -0500
@@ -0,0 +1,144 @@
+#include <math.h>
+#include "R.h"
+#include "Rmath.h"
+#include "Rinternals.h"
+
+
+#undef DEBUG 1
+
+// dout is npos-length output array.
+// n - number of positions in pos (and length of tc count array)
+// spos - starting position
+void cdensum(int *n, double *pos, double *tc, double *spos, int *bw,int *dw, int *npos, int *step,double *dout)
+{
+  int i,j;
+ 
+  double epos= *spos + ((double) *npos);
+  double dbw=(double) *bw;
+  for(i = 0; i< *n; i++) {
+    // size of the window to which the contributions should be added
+    int in=(int) (pos[i]- *spos);
+    int ic=tc[i];
+    int whs=(*dw)*(*bw)*ic;
+    int ws=(int) floor((in-whs)/(*step));
+    int we=(int) ceil((in+whs)/(*step));
+    if(ws<0) { ws=0; } 
+    if(we>= *npos) { we= *npos -1; }
+    
+    for(j=ws;j<we;j++) {
+      double beta=((double)(j*(*step)-in))/dbw;
+      dout[j]+=((double)ic)*exp(-0.5*beta*beta);
+    }
+  }
+}
+
+
+// window tag counts
+// dout is npos-length output array that will contain window tag counts
+// windows are of a specified size, moved at a specified step
+// n - number of positions in sorted tag array (positive only)
+// spos - starting position
+void window_n_tags(int *n, double *pos, double *spos, int *window_size, int *window_step, int *npos, int *dout)
+{
+  int i;
+  int cs=0; int ce=0; // current array start/end indecies
+  int ctc=0; // current tag count
+  double wpos=*spos-(*window_size)/2; // left-edge position
+  //Rprintf("n=%d; window_size=%d, window_step=%d, npos=%d, spos=%f\n",*n,*window_size,*window_step,*npos,*spos);
+  for(i=0;i<*npos;i++) {
+    // advance end if needed
+    double ep=wpos+(*window_size);
+    while(ce<(*n) && pos[ce]<=ep) {
+      ctc++; ce++;
+    }
+    // advance start
+    while(cs<*n && pos[cs]<wpos) {
+      ctc--; cs++;
+    }
+    dout[i]=ctc;
+    // advance window position
+    wpos+=*window_step;
+  }
+}
+
+// window tag counts
+// windows are of a specified size, moved at a specified step
+// pos - tag positions (positive, pre-shifted)y
+// spos - starting position
+// returns nsteps-length output array that will contain window tag counts
+SEXP cwindow_n_tags(SEXP pos_R, SEXP spos_R, SEXP window_size_R, SEXP window_step_R, SEXP nsteps_R) {
+  double* pos=REAL(pos_R);
+  int n=LENGTH(pos_R);
+  int window_size=*INTEGER(window_size_R);
+  int window_step=*INTEGER(window_step_R);
+  int nsteps=*INTEGER(nsteps_R);
+  double spos=*REAL(spos_R);
+  
+  // allocate return array
+  SEXP tc_R;
+  PROTECT(tc_R=allocVector(INTSXP,nsteps));
+  int* dout=INTEGER(tc_R);
+
+  int i;
+  int cs=0; int ce=0; // current array start/end indecies
+  int ctc=0; // current tag count
+  double wpos=spos-window_size/2; // left-edge position
+  //Rprintf("n=%d; window_size=%d, window_step=%d, npos=%d, spos=%f\n",n,window_size,window_step,nsteps,spos);
+  for(i=0;i<nsteps;i++) {
+    // advance end if needed
+    double ep=wpos+window_size;
+    while(ce<n && pos[ce]<=ep) {
+      ctc++; ce++;
+    }
+    // advance start
+    while(cs<n && pos[cs]<wpos) {
+      ctc--; cs++;
+    }
+    dout[i]=ctc;
+    // advance window position
+    wpos+=window_step;
+  }
+  UNPROTECT(1);
+  return(tc_R);
+}
+
+// tag counts in windows around specified positions
+// pos - tag positions 
+// ntags - number of tags in each position
+// wpos - window positions
+// returns a pos-length vector giving number of tags that fall within window_half_size from the provided positions
+SEXP cwindow_n_tags_around(SEXP pos_R, SEXP ntags_R, SEXP wpos_R, SEXP window_half_size_R) {
+  double* pos=REAL(pos_R);
+  int* ntags=INTEGER(ntags_R);
+  int n=LENGTH(pos_R);
+  double* wpos=REAL(wpos_R);
+  int nw=LENGTH(wpos_R); // number of windows
+  double whs=(double) *INTEGER(window_half_size_R);
+  
+  // allocate return array
+  SEXP tc_R;
+  PROTECT(tc_R=allocVector(INTSXP,nw));
+  int* dout=INTEGER(tc_R);
+
+  int i;
+  int cs=0; int ce=0; // current array start/end indecies
+  int ctc=0; // current tag count
+  for(i=0;i<nw;i++) {
+    //if(i>(nw-2)) {      Rprintf("-i=%d; cs=%d, ce=%d; ctc=%d\n",i,cs,ce,ctc);    }    
+    // advance end if needed
+    double ep=wpos[i]+whs;
+    while(ce<n && pos[ce]<=ep) {
+      ctc+=ntags[ce]; ce++;
+    }
+    // advance start
+    double sp=wpos[i]-whs;
+    while(cs<n && pos[cs]<sp) {
+      ctc-=ntags[cs]; cs++;
+    }
+    dout[i]=ctc;
+    // if(i>(nw-2)) {      Rprintf("+i=%d; cs=%d, ce=%d; ctc=%d\n",i,cs,ce,ctc);    }
+  }
+  UNPROTECT(1);
+  return(tc_R);
+}
+