|
6
|
1 #include <math.h>
|
|
|
2 #include "R.h"
|
|
|
3 #include "Rmath.h"
|
|
|
4 #include "Rinternals.h"
|
|
|
5
|
|
|
6
|
|
|
7 #undef DEBUG 1
|
|
|
8
|
|
|
9 // dout is npos-length output array.
|
|
|
10 // n - number of positions in pos (and length of tc count array)
|
|
|
11 // spos - starting position
|
|
|
12 void cdensum(int *n, double *pos, double *tc, double *spos, int *bw,int *dw, int *npos, int *step,double *dout)
|
|
|
13 {
|
|
|
14 int i,j;
|
|
|
15
|
|
|
16 double epos= *spos + ((double) *npos);
|
|
|
17 double dbw=(double) *bw;
|
|
|
18 for(i = 0; i< *n; i++) {
|
|
|
19 // size of the window to which the contributions should be added
|
|
|
20 int in=(int) (pos[i]- *spos);
|
|
|
21 int ic=tc[i];
|
|
|
22 int whs=(*dw)*(*bw)*ic;
|
|
|
23 int ws=(int) floor((in-whs)/(*step));
|
|
|
24 int we=(int) ceil((in+whs)/(*step));
|
|
|
25 if(ws<0) { ws=0; }
|
|
|
26 if(we>= *npos) { we= *npos -1; }
|
|
|
27
|
|
|
28 for(j=ws;j<we;j++) {
|
|
|
29 double beta=((double)(j*(*step)-in))/dbw;
|
|
|
30 dout[j]+=((double)ic)*exp(-0.5*beta*beta);
|
|
|
31 }
|
|
|
32 }
|
|
|
33 }
|
|
|
34
|
|
|
35
|
|
|
36 // window tag counts
|
|
|
37 // dout is npos-length output array that will contain window tag counts
|
|
|
38 // windows are of a specified size, moved at a specified step
|
|
|
39 // n - number of positions in sorted tag array (positive only)
|
|
|
40 // spos - starting position
|
|
|
41 void window_n_tags(int *n, double *pos, double *spos, int *window_size, int *window_step, int *npos, int *dout)
|
|
|
42 {
|
|
|
43 int i;
|
|
|
44 int cs=0; int ce=0; // current array start/end indecies
|
|
|
45 int ctc=0; // current tag count
|
|
|
46 double wpos=*spos-(*window_size)/2; // left-edge position
|
|
|
47 //Rprintf("n=%d; window_size=%d, window_step=%d, npos=%d, spos=%f\n",*n,*window_size,*window_step,*npos,*spos);
|
|
|
48 for(i=0;i<*npos;i++) {
|
|
|
49 // advance end if needed
|
|
|
50 double ep=wpos+(*window_size);
|
|
|
51 while(ce<(*n) && pos[ce]<=ep) {
|
|
|
52 ctc++; ce++;
|
|
|
53 }
|
|
|
54 // advance start
|
|
|
55 while(cs<*n && pos[cs]<wpos) {
|
|
|
56 ctc--; cs++;
|
|
|
57 }
|
|
|
58 dout[i]=ctc;
|
|
|
59 // advance window position
|
|
|
60 wpos+=*window_step;
|
|
|
61 }
|
|
|
62 }
|
|
|
63
|
|
|
64 // window tag counts
|
|
|
65 // windows are of a specified size, moved at a specified step
|
|
|
66 // pos - tag positions (positive, pre-shifted)y
|
|
|
67 // spos - starting position
|
|
|
68 // returns nsteps-length output array that will contain window tag counts
|
|
|
69 SEXP cwindow_n_tags(SEXP pos_R, SEXP spos_R, SEXP window_size_R, SEXP window_step_R, SEXP nsteps_R) {
|
|
|
70 double* pos=REAL(pos_R);
|
|
|
71 int n=LENGTH(pos_R);
|
|
|
72 int window_size=*INTEGER(window_size_R);
|
|
|
73 int window_step=*INTEGER(window_step_R);
|
|
|
74 int nsteps=*INTEGER(nsteps_R);
|
|
|
75 double spos=*REAL(spos_R);
|
|
|
76
|
|
|
77 // allocate return array
|
|
|
78 SEXP tc_R;
|
|
|
79 PROTECT(tc_R=allocVector(INTSXP,nsteps));
|
|
|
80 int* dout=INTEGER(tc_R);
|
|
|
81
|
|
|
82 int i;
|
|
|
83 int cs=0; int ce=0; // current array start/end indecies
|
|
|
84 int ctc=0; // current tag count
|
|
|
85 double wpos=spos-window_size/2; // left-edge position
|
|
|
86 //Rprintf("n=%d; window_size=%d, window_step=%d, npos=%d, spos=%f\n",n,window_size,window_step,nsteps,spos);
|
|
|
87 for(i=0;i<nsteps;i++) {
|
|
|
88 // advance end if needed
|
|
|
89 double ep=wpos+window_size;
|
|
|
90 while(ce<n && pos[ce]<=ep) {
|
|
|
91 ctc++; ce++;
|
|
|
92 }
|
|
|
93 // advance start
|
|
|
94 while(cs<n && pos[cs]<wpos) {
|
|
|
95 ctc--; cs++;
|
|
|
96 }
|
|
|
97 dout[i]=ctc;
|
|
|
98 // advance window position
|
|
|
99 wpos+=window_step;
|
|
|
100 }
|
|
|
101 UNPROTECT(1);
|
|
|
102 return(tc_R);
|
|
|
103 }
|
|
|
104
|
|
|
105 // tag counts in windows around specified positions
|
|
|
106 // pos - tag positions
|
|
|
107 // ntags - number of tags in each position
|
|
|
108 // wpos - window positions
|
|
|
109 // returns a pos-length vector giving number of tags that fall within window_half_size from the provided positions
|
|
|
110 SEXP cwindow_n_tags_around(SEXP pos_R, SEXP ntags_R, SEXP wpos_R, SEXP window_half_size_R) {
|
|
|
111 double* pos=REAL(pos_R);
|
|
|
112 int* ntags=INTEGER(ntags_R);
|
|
|
113 int n=LENGTH(pos_R);
|
|
|
114 double* wpos=REAL(wpos_R);
|
|
|
115 int nw=LENGTH(wpos_R); // number of windows
|
|
|
116 double whs=(double) *INTEGER(window_half_size_R);
|
|
|
117
|
|
|
118 // allocate return array
|
|
|
119 SEXP tc_R;
|
|
|
120 PROTECT(tc_R=allocVector(INTSXP,nw));
|
|
|
121 int* dout=INTEGER(tc_R);
|
|
|
122
|
|
|
123 int i;
|
|
|
124 int cs=0; int ce=0; // current array start/end indecies
|
|
|
125 int ctc=0; // current tag count
|
|
|
126 for(i=0;i<nw;i++) {
|
|
|
127 //if(i>(nw-2)) { Rprintf("-i=%d; cs=%d, ce=%d; ctc=%d\n",i,cs,ce,ctc); }
|
|
|
128 // advance end if needed
|
|
|
129 double ep=wpos[i]+whs;
|
|
|
130 while(ce<n && pos[ce]<=ep) {
|
|
|
131 ctc+=ntags[ce]; ce++;
|
|
|
132 }
|
|
|
133 // advance start
|
|
|
134 double sp=wpos[i]-whs;
|
|
|
135 while(cs<n && pos[cs]<sp) {
|
|
|
136 ctc-=ntags[cs]; cs++;
|
|
|
137 }
|
|
|
138 dout[i]=ctc;
|
|
|
139 // if(i>(nw-2)) { Rprintf("+i=%d; cs=%d, ce=%d; ctc=%d\n",i,cs,ce,ctc); }
|
|
|
140 }
|
|
|
141 UNPROTECT(1);
|
|
|
142 return(tc_R);
|
|
|
143 }
|
|
|
144
|
