diff spp/src/peaks.cpp @ 6:ce08b0efa3fd draft

Uploaded

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/spp/src/peaks.cpp	Tue Nov 27 16:11:40 2012 -0500
@@ -0,0 +1,804 @@
+#include <vector>
+#include <string.h>
+#include <iostream>
+#include <string>
+#include <set>
+
+extern "C" {
+#include "R.h"
+#include "Rmath.h"
+#include "Rinternals.h"
+#include "Rdefines.h"
+}
+
+using namespace std;
+using namespace __gnu_cxx; 
+
+/**
+ * Calculate all local peaks
+ */
+
+//#define DEBUG 1
+
+extern "C" {
+  SEXP find_peaks(SEXP x_R,SEXP thr_R,SEXP max_span_R) {
+
+#ifdef DEBUG  
+    Rprintf("start\n");
+#endif
+    double* x=REAL(x_R);
+    int nx=LENGTH(x_R);
+    int max_span=*INTEGER(max_span_R);
+    double thr=REAL(thr_R)[0];
+#ifdef DEBUG  
+    Rprintf("n=%d; thr=%f; max_span=%d\n",nx,thr,max_span);
+#endif
+
+    vector<int> pos;
+  
+    double pv=x[0];
+    double ppv=0; // previous peak value
+    int ppp=-max_span-1; // previous peak position
+    
+    for(int i=1;i<(nx-1);i++) {
+      if(x[i]>pv && x[i]>=thr && x[i]>x[i+1]) {
+	if(max_span>2) {
+	  //Rprintf("i=%d; ppp=%d\n",i,ppp);
+	  if(i-ppp > max_span) {
+	    if(ppp>=0) {
+	      pos.push_back(ppp);
+	    }
+	    //Rprintf("recorded %d; now %d\n",ppp,i);
+	    ppp=i; ppv=x[i];
+	  } else {
+	    if(x[i]>ppv) {
+	      //Rprintf("reset from %d to %d\n",ppp,i);
+	      ppp=i; ppv=x[i];
+	    }
+	  }
+	} else {
+	  pos.push_back(i);
+	}
+      }
+      if(x[i]!=x[i+1]) { pv=x[i]; }
+    }
+
+    // add remaining peak
+    if(max_span>2 && ppp>=0) {
+      pos.push_back(ppp);
+    }
+
+    SEXP nv;
+    PROTECT(nv=allocVector(INTSXP,pos.size())); 
+    int* i_nv=INTEGER(nv);
+    int i=0;
+    for(vector<int> ::const_iterator pi=pos.begin();pi!=pos.end();++pi) {
+      i_nv[i++]=1+(*pi);
+    }
+  
+    UNPROTECT(1);
+    return(nv);
+  }
+
+
+
+
+  /************************************************************************/
+  // given a data vector d (positive values) and a set of signed center coordinates pos,
+  // returns coordinates of data points relative to the centers
+  // size is the size of the region around the centers
+  // return: vector of relative coordinates (x) and indecies of centers relative the coordinate
+  // was calculated (i).
+  SEXP get_relative_coordinates(SEXP d_R,
+				SEXP pos_R,
+				SEXP size_R)
+  {
+    int *d, *pos; 
+    int npos,nd,size;
+  
+    d = INTEGER(d_R); pos = INTEGER(pos_R);
+    npos=LENGTH(pos_R);  nd=LENGTH(d_R);
+    size = INTEGER(size_R)[0];
+#ifdef DEBUG  
+    Rprintf("|d|=%d, |c|=%d, size=%d\n",nd,npos,size);
+#endif
+
+    vector<int> x; vector<int> xi;
+    int k=0; // current pos index
+    
+    for(int i=0;i<nd;i++) {
+      // increment k until pos[k]+size>=d[i]
+      while((abs(pos[k])+size) < d[i]) { k++; if(k==npos) { break; };
+#ifdef DEBUG  
+	Rprintf("advancing k to %d\n",k);
+#endif
+      }
+      if(k==npos) { break; };
+      // increment i until d[i]>=pos[k]-size
+      while((abs(pos[k])-size) > d[i]) { i++; if(i==nd) { break; }
+#ifdef DEBUG  
+	Rprintf("advancing i to %d\n",i);
+#endif
+      }
+      if(i==nd) { break; }
+
+
+      int l=k;
+      while((l<npos) && ((abs(pos[l])-size) <= d[i])) { l++; 
+#ifdef DEBUG  
+	Rprintf("advancing l to %d\n",l);
+#endif
+      }
+      for(int j=k;j<l;j++) {
+	int pd=d[i]-abs(pos[j]);
+	if(abs(pd)<=size) {
+	  // record
+	  if(pos[j]>0) {
+	    x.push_back(pd);
+	  } else {
+	    x.push_back(-1*pd);
+	  }
+	  xi.push_back(j);
+#ifdef DEBUG  	
+	  Rprintf("recorded i=%d, j=%d\n",i,j);
+#endif
+	} else {
+	  break;
+	}
+      }
+    }
+    
+    SEXP xv_R,xiv_R;
+    PROTECT(xv_R=allocVector(INTSXP,x.size())); 
+    PROTECT(xiv_R=allocVector(INTSXP,x.size())); 
+    int* xv=INTEGER(xv_R);
+    int* xiv=INTEGER(xiv_R);
+
+    int i=0;
+    for(vector<int> ::const_iterator pi=x.begin();pi!=x.end();++pi) {
+      xv[i++]=*pi;
+    }
+    i=0;
+    for(vector<int> ::const_iterator pi=xi.begin();pi!=xi.end();++pi) {
+      xiv[i++]=1+(*pi);
+    }
+    
+    SEXP ans_R, names_R;
+    PROTECT(names_R = allocVector(STRSXP, 2));
+    SET_STRING_ELT(names_R, 0, mkChar("x"));
+    SET_STRING_ELT(names_R, 1, mkChar("i"));
+    
+    PROTECT(ans_R = allocVector(VECSXP, 2));
+    SET_VECTOR_ELT(ans_R, 0, xv_R);
+    SET_VECTOR_ELT(ans_R, 1, xiv_R);
+    setAttrib(ans_R, R_NamesSymbol, names_R);
+  
+    UNPROTECT(4);
+    return(ans_R);
+  }
+
+
+  // determines a set of points within a set of fragments
+  // note: all vectors sorted in ascending order
+  // note: all vectors are integers
+  // x_R - vector of point positions 
+  // se_R - vector of start and end positions 
+  // fi_R - vector of signed fragment indecies
+  // return_list_R - whether a list of fragments should be returned for each point
+  // return_unique_R - whether points in multiple fragments should be omitted
+  SEXP points_within(SEXP x_R,SEXP se_R,SEXP fi_R,SEXP return_list_R,SEXP return_unique_R,SEXP return_point_counts_R) {
+#ifdef DEBUG  
+    Rprintf("start\n");
+#endif
+    int* x=INTEGER(x_R);
+    int nx=LENGTH(x_R);
+    int* se=INTEGER(se_R);
+    int* fi=INTEGER(fi_R);
+    int nf=LENGTH(se_R);
+
+    int return_list=*(INTEGER(return_list_R));
+    int return_unique=*(INTEGER(return_unique_R));
+    int return_point_counts=*(INTEGER(return_point_counts_R));
+
+#ifdef DEBUG  
+    Rprintf("nf=%d; nx=%d, return_list=%d, return_unique=%d, return_point_counts=%d\n",nf/2,nx,return_list,return_unique,return_point_counts);
+#endif
+    set<int> fset;
+
+
+    SEXP nv; int *i_nv;
+    int np=0;
+    if(return_point_counts) {
+      PROTECT(nv = allocVector(INTSXP, nf/2)); np++;      
+      i_nv=INTEGER(nv);
+      for(int i=0;i<nf/2;i++) { i_nv[i]=0; }
+    } else if(return_list) {
+      PROTECT(nv = allocVector(VECSXP, nx)); np++;
+    } else {
+      PROTECT(nv=allocVector(INTSXP,nx));  np++;
+      i_nv=INTEGER(nv);
+    }
+
+    int j=0;
+
+    for(int i=0;i<nx;i++) {
+      // advance j
+      while(j<nf && se[j]<x[i]) {
+	int frag=fi[j];
+	if(frag>0) { // insert
+	  fset.insert(frag);
+#ifdef DEBUG  
+	  Rprintf("inserted frag %d, size=%d\n",frag,fset.size());
+#endif
+	} else { // remove
+	  fset.erase(-frag);
+#ifdef DEBUG  
+	  Rprintf("removed frag %d, size=%d\n",-frag,fset.size());
+#endif
+	}
+	j++;
+      }
+#ifdef DEBUG  
+      Rprintf("i=%d j=%d\n",i,j);
+#endif
+      if(return_list) {
+	if(fset.empty() || (return_unique && fset.size()>1)) {
+	  // assign null list?
+	} else {
+	  SEXP fil_R;
+	  PROTECT(fil_R=allocVector(INTSXP,fset.size()));  np++;
+	  int* fil=INTEGER(fil_R);
+	  int k=0;
+	  for(set<int>::const_iterator ki=fset.begin();ki!=fset.end();++ki) {
+	    fil[k]=*ki; k++;
+	  }
+	  SET_VECTOR_ELT(nv, i, fil_R);
+	  UNPROTECT(1); np--;
+	}
+      } else {
+	if(return_point_counts) {
+	  for(set<int>::const_iterator ki=fset.begin();ki!=fset.end();++ki) {
+	    i_nv[*ki-1]++;
+	  }
+	} else {
+	  if(fset.empty() || (return_unique && fset.size()>1)) {
+	    i_nv[i]=-1;
+	  } else {
+	    i_nv[i]=*fset.begin();
+	  }
+	}
+      }
+    }
+
+    UNPROTECT(np);
+    return(nv);
+  }
+
+
+  SEXP expuni_lr(SEXP x_R,      // positions and their number (assumed sorted in ascending order)
+		 SEXP mdist_R,  // max distance at which points should be considered
+		 SEXP lambda_R,  // lambda value
+		 SEXP spos_R,  // starting position
+		 SEXP epos_R,  // ending position
+		 SEXP step_R,  // step size
+		 SEXP return_peaks_R, // whether peak positions should be returned, or entire score vector
+		 SEXP min_peak_lr_R // min peak height (lr)
+		 ) 
+  {
+
+#ifdef DEBUG  
+    Rprintf("start\n");
+#endif
+    int* x=INTEGER(x_R);
+    int nx=LENGTH(x_R);
+    int mdist=INTEGER(mdist_R)[0];
+    double lambda=*(REAL(lambda_R));
+
+    int return_peaks=*(INTEGER(return_peaks_R));
+    double min_peak=*(REAL(min_peak_lr_R));
+
+    int spos=*(INTEGER(spos_R));
+    int epos=*(INTEGER(epos_R));
+    int step=*(INTEGER(step_R));
+
+    int nsteps=(int) (epos-spos)/step;
+
+
+#ifdef DEBUG  
+    Rprintf("n=%d; lambda=%f; mdist=%d; spos=%d; epos=%d; step=%d; nsteps=%d\n",nx,lambda,mdist,spos,epos,step,nsteps);
+#endif
+
+    
+    SEXP nv;
+    double *d_nv;
+    if(!return_peaks) {
+      PROTECT(nv=allocVector(REALSXP,nsteps+1)); 
+      d_nv=REAL(nv);
+    }
+
+
+    int i=0; // current index of the first point being used in the calculations
+    int j=0; // current index of the last point being used in the calculations
+    int sx=0; // current sum of all positions
+    int n=0;
+    
+    for(int k=0; k<=nsteps; k++) {
+      int cpos=spos+k*step;
+      // increase i until x[i]>=cpos-mdist; remove x from sx; decrement n;
+      while(i<nx && x[i]<(cpos-mdist)) {
+	n--; sx-=x[i]; i++;
+	//Rprintf("incremented i: i=%d; n=%d; sx=%d; cpos-mdist=%d; x[i-1]=%d\n",i,n,sx,cpos-mdist,x[i-1]);
+      }
+      //Rprintf("stable i: i=%d; n=%d; sx=%d; cpos-mdist=%d; x[i-1]=%d\n",i,n,sx,cpos-mdist,x[i-1]);
+
+      //if(i>j) { j=i; }
+
+      // increase j until x[j]>cpos
+      while(j<nx && x[j]<=cpos) {
+	n++; sx+=x[j]; j++;
+	//Rprintf("incremented j: j=%d; n=%d; sx=%d; cpos=%d; x[j-1]=%d\n",j,n,sx,cpos,x[j-1]);
+      }
+      //Rprintf("stable j: j=%d; n=%d; sx=%d; cpos=%d; x[j-1]=%d\n",j,n,sx,cpos,x[j]);
+      
+      // calculate lr
+      d_nv[k]=((double)(1-n))*log(lambda)-lambda*((double)(n*(cpos+1)-sx));
+      //Rprintf("recorded lr[%d]=%f\n",k-1,d_nv[k-1]);
+    }
+    UNPROTECT(1);
+    return(nv);
+  }
+
+
+  SEXP allpdist(SEXP x_R,SEXP max_dist_R) {
+
+#ifdef DEBUG  
+    Rprintf("start\n");
+#endif
+    double* x=REAL(x_R);
+    int nx=LENGTH(x_R);
+    double max_dist=*REAL(max_dist_R);
+#ifdef DEBUG  
+    Rprintf("n=%d; max_dist=%d\n",nx,max_dist);
+#endif
+
+    vector<double> dist;
+    
+    for(int i=0;i<nx;i++) {
+      for(int j=i+1;j<nx;j++) {
+
+	double d=x[j]-x[i];
+#ifdef DEBUG  
+	Rprintf("i=%d; j=%d; d=%f\n",i,j,d);
+#endif
+	if(d<=max_dist) {
+	  dist.push_back(d);
+	} else {
+	  break;
+	}
+      }
+    }
+
+    SEXP nv;
+    PROTECT(nv=allocVector(REALSXP,dist.size())); 
+    double* i_nv=REAL(nv);
+    int i=0;
+    for(vector<double> ::const_iterator pi=dist.begin();pi!=dist.end();++pi) {
+      i_nv[i++]=*pi;
+    }
+  
+    UNPROTECT(1);
+    return(nv);
+  }
+
+  // same as above, but for two different sets
+  SEXP allxpdist(SEXP x_R,SEXP y_R, SEXP max_dist_R) {
+
+#ifdef DEBUG  
+    Rprintf("start\n");
+#endif
+    double* x=REAL(x_R);
+    double* y=REAL(y_R);
+    int nx=LENGTH(x_R);
+    int ny=LENGTH(y_R);
+    double max_dist=*REAL(max_dist_R);
+#ifdef DEBUG  
+    Rprintf("nx=%d; ny=%d; max_dist=%d\n",nx,ny,max_dist);
+#endif
+
+    vector<double> dist;
+    int yi=0; // latest y start index
+
+    for(int i=0;i<nx;i++) {
+      // adjust yi so that yi>=x[i]-max_dist_R
+      while(y[yi]<(x[i]-max_dist) && yi<ny) { yi++; }
+      if(yi==ny) { break; }
+
+      for(int j=yi;j<ny;j++) {
+        double d=y[j]-x[i];
+#ifdef DEBUG  
+        Rprintf("i=%d; j=%d; d=%f\n",i,j,d);
+#endif
+        if(d<=max_dist) {
+          dist.push_back(d);
+        } else {
+          break;
+        }
+      }
+    }
+
+    SEXP nv;
+    PROTECT(nv=allocVector(REALSXP,dist.size()));
+    double* i_nv=REAL(nv);
+    int i=0;
+    for(vector<double> ::const_iterator pi=dist.begin();pi!=dist.end();++pi) {
+      i_nv[i++]=*pi;
+    }
+
+    UNPROTECT(1);
+    return(nv);
+  }
+
+  // returns a vector giving for each point,
+  // number of points within a given max_dist
+  SEXP nwithindist(SEXP x_R,SEXP max_dist_R) {
+
+#ifdef DEBUG  
+    Rprintf("start\n");
+#endif
+    double* x=REAL(x_R);
+    int nx=LENGTH(x_R);
+    double max_dist=*REAL(max_dist_R);
+
+    SEXP nv;
+    PROTECT(nv=allocVector(REALSXP,nx)); 
+    double* i_nv=REAL(nv);
+    for(int i=0;i<nx;i++) { i_nv[i]=0; }
+
+#ifdef DEBUG  
+    Rprintf("n=%d; max_dist=%d\n",nx,max_dist);
+#endif
+
+    for(int i=0;i<nx;i++) {
+      for(int j=i+1;j<nx;j++) {
+
+	double d=x[j]-x[i];
+#ifdef DEBUG  
+	Rprintf("i=%d; j=%d; d=%f\n",i,j,d);
+#endif
+	if(d<=max_dist) {
+	  i_nv[i]++;
+	  i_nv[j]++;
+	} else {
+	  break;
+	}
+      }
+    }
+  
+    UNPROTECT(1);
+    return(nv);
+  }
+
+
+
+
+  // given a list of sorted chromosome signal and background vectors (unscaled), determine 
+  // cluster contigs exceeding thr poisson P value, based on a whs window size,
+  // and satisfying mcs cluster size
+  SEXP find_poisson_enrichment_clusters(SEXP pos_R,SEXP flag_R,SEXP wsize_R,SEXP thr_R,SEXP mcs_R,SEXP bgm_R,SEXP mintag_R,SEXP either_R) {
+
+#ifdef DEBUG  
+    Rprintf("start\n");
+#endif
+    double* pos=REAL(pos_R);
+    int* flag=INTEGER(flag_R);
+    int nt=LENGTH(pos_R);
+    
+    int mcs=*INTEGER(mcs_R);
+    int wsize=*INTEGER(wsize_R);
+    int either=*INTEGER(either_R);
+    double thr=REAL(thr_R)[0];
+    double bgm=REAL(bgm_R)[0];
+    double mintag=REAL(mintag_R)[0];
+
+#ifdef DEBUG  
+    Rprintf("nt=%d; wsize=%d; thr=%f; mcs=%d; min.tag=%f; bgm=%f\n",nt,wsize,thr,mcs,mintag,bgm);
+#endif
+    
+    
+    vector< pair<double,double> > contigs;
+    
+    // running indecies (start and end)
+    int si=0; 
+    int ei=0;
+    
+    // current window coordinate
+    double ws=pos[0];
+    
+    // current window tag counts
+    int cc[2]={0,0};
+
+
+    if(nt>0) {
+      cc[flag[si]]++;    
+      // increment window end
+      while(ei<(nt-1) && (pos[ei+1]-ws) <= wsize) {
+	ei++;
+	cc[flag[ei]]++;
+      }
+
+
+      // cluster start,end positions
+      double cs,ce;
+      int inclust=0;
+
+      while(si<nt-1) {
+      
+	if((pos[si+1]-ws) > (pos[ei+1] - ws - wsize) && ei!=(nt-1)) {
+	  // move end boudnary
+	  ei++;
+	  ws=pos[ei]-wsize;
+	  cc[flag[ei]]++;
+	  while(ei<(nt-1) && pos[ei+1]==ws+wsize) {
+	    ei++;
+	    cc[flag[ei]]++;
+	  }
+	
+	  // increment window start
+	  while(si<(nt-1) && pos[si] < ws) {
+	    cc[flag[si]]--;
+	    si++;
+	  }
+
+	} else {
+	  // move up start boundary
+	  ws=pos[si+1];
+	  cc[flag[si]]--;
+	  si++;
+	  while(si<(nt-1) && pos[si+1]==ws) {
+	    cc[flag[si]]--;
+	    si++;
+	  }
+	
+	  // increment window end
+	  while(ei<(nt-1) && (pos[ei+1] - ws) <= wsize) {
+	    ei++;
+	    cc[flag[ei]]++;
+	  }
+
+	}
+
+	// calculate z score
+	double dc0=((double)cc[0])+0.5;
+	double dc1=((double)cc[1])+0.5;
+	double rte=dc0+dc1-0.25*thr*thr;
+	double lb;
+	if(rte<=0) { 
+	  lb=0; 
+	} else {
+	  lb=(sqrt(dc1*dc0) - 0.5*thr*sqrt(rte))/(dc0 - 0.25*thr*thr);
+	  if(lb<0) { lb=0; }
+	  lb*=lb;
+	}
+
+	//Rprintf("%f=f(%f,%f,%f); %f=f(%f,%f,%f)\n",lb,1.0-thr,2.0*dc1,2.0*dc0,ub,thr,2.0*dc1,2.0*dc0);
+      
+#ifdef DEBUG  
+	//double ub=gsl_cdf_fdist_Qinv(thr,2.0*dc1,2.0*dc0)*dc1/dc0;
+	double ub=(sqrt(dc1*dc0) + 0.5*thr*sqrt(rte))/(dc0 - 0.25*thr*thr);
+	ub*=ub;
+	Rprintf("s=%d (%f); e=%d (%f); window: %f-%f; cc=[%d,%d]; lb=%f; ub=%f\n",si,pos[si],ei,pos[ei],ws,ws+wsize,cc[0],cc[1],lb,ub);
+#endif
+      
+	int bc=lb>=bgm && cc[1]>=mintag;
+	if(either) {
+	  bc=lb>=bgm || cc[1]>=mintag;
+	}
+	if(bc) {
+	  if(inclust) {
+	    double nce=ws+wsize/2.0;
+	    if(nce-ce > wsize/2.0) {
+	      // next point is too far removed, end cluster
+	      if(ce-cs >= mcs) {
+		contigs.push_back(pair<double,double>(cs,ce));
+#ifdef DEBUG  
+		Rprintf("recorded cluster %f-%f\n",cs,ce);
+#endif
+	      }
+	      inclust=0;
+	    } else {
+	      ce=nce;
+	    }
+	  } else {
+	    inclust=1;
+	    cs=ws+wsize/2.0;
+	    ce=cs;
+	  }
+	} else {
+	  if(inclust) {
+	    if(ce-cs >= mcs) {
+	      contigs.push_back(pair<double,double>(cs,ce));
+#ifdef DEBUG  
+	      Rprintf("recorded cluster %f-%f\n",cs,ce);
+#endif
+	    }
+	    inclust=0;
+	  }
+	}
+    
+      }
+
+      if(inclust) {
+	if(ce-cs >= mcs) {
+	  contigs.push_back(pair<double,double>(cs,ce));
+#ifdef DEBUG  
+	  Rprintf("recorded cluster %f-%f\n",cs,ce);
+#endif
+	}
+	inclust=0;
+      }
+    }
+    
+    SEXP cs_R,ce_R;
+    PROTECT(cs_R=allocVector(REALSXP,contigs.size())); 
+    PROTECT(ce_R=allocVector(REALSXP,contigs.size())); 
+    double* csa=REAL(cs_R);
+    double* cea=REAL(ce_R);
+
+    int i=0;
+    for(vector< pair<double,double> >::const_iterator ci=contigs.begin(); ci!=contigs.end();++ci) {
+      csa[i]=ci->first;
+      cea[i]=ci->second;
+      i++;
+    }
+    
+    SEXP ans_R, names_R;
+    PROTECT(names_R = allocVector(STRSXP, 2));
+    SET_STRING_ELT(names_R, 0, mkChar("s"));
+    SET_STRING_ELT(names_R, 1, mkChar("e"));
+    
+    PROTECT(ans_R = allocVector(VECSXP, 2));
+    SET_VECTOR_ELT(ans_R, 0, cs_R);
+    SET_VECTOR_ELT(ans_R, 1, ce_R);
+    setAttrib(ans_R, R_NamesSymbol, names_R);
+  
+    UNPROTECT(4);
+    return(ans_R);
+
+  }
+
+
+  // finds intersection between a list of regions
+  // the flag has +n/-n value, corresponding to the start/end of a segment in n-th regionset
+  // max_val: 1 - report max overlapping value, -1: report min, 0 - don't look at values
+  // returns: $s, $e, ($v) lists
+  SEXP region_intersection(SEXP n_R,SEXP pos_R,SEXP flags_R,SEXP vals_R,SEXP max_val_R,SEXP union_R) {
+    const int max_val=*INTEGER(max_val_R);
+    const int unionr=*INTEGER(union_R);
+    const int n=*INTEGER(n_R);
+    double* pos=REAL(pos_R);
+    int* flags=INTEGER(flags_R);
+    double* val=REAL(vals_R);
+    
+#ifdef DEBUG  
+    Rprintf("n=%d; npos=%d; max_val=%d\n",n,LENGTH(pos_R),max_val);
+#endif
+
+    int s[n]; // flag status for each set
+    double mv[n]; // max/min value of current clusters
+
+    for(int i=0;i<n;i++) { s[i]=0; }
+    
+    vector<double> starts;
+    vector<double> ends;
+    vector<double> values;
+
+    int start=-1;
+    double mval=0;
+    for(int i=0;i<LENGTH(pos_R);i++) {
+      // update flags
+      int f=flags[i];
+      if(f>0) {
+	s[abs(f)-1]++;
+      } else {
+	s[abs(f)-1]--;
+      }
+      
+      if(max_val!=0 && val[i]*max_val > mval*max_val) { mval=val[i]; }
+
+      // joined status
+      int all;
+      if(unionr) {
+	all=0;
+	for(int j=0;j<n;j++) { if(s[j]>0) { all=1; break;} }
+      } else {
+	all=1;
+	for(int j=0;j<n;j++) { all=all & (s[j]>0); }
+      }
+      
+      
+      //Rprintf("i=%d; s=[",i);
+      //for(int j=0;j<n;j++) { Rprintf("%d",s[j]); }
+      //Rprintf("]; all=%d; start=%d\n",all,start);
+
+      if(start>=0) {
+	// in fragment
+	if(!all) { 
+	  // end fragment
+	  starts.push_back(pos[start]);
+	  ends.push_back(pos[i]);
+	  start=-1;
+	  if(max_val!=0) { values.push_back(mval); }
+
+#ifdef DEBUG  
+	  Rprintf("recorded new fragment (s=%f,e=%f,v=%f);\n",pos[start],pos[i],mval);
+#endif
+	}
+      } else {
+	// should a fragment be started?
+	if(all) {
+	  start=i;
+	  if(max_val!=0) { mval=val[i]; }
+#ifdef DEBUG  
+	  Rprintf("starting new fragment (s=%f,i=%d);\n",pos[start],i);
+#endif
+	}
+      }
+    }
+    SEXP cs_R,ce_R,cv_R;
+    PROTECT(cs_R=allocVector(REALSXP,starts.size())); 
+    PROTECT(ce_R=allocVector(REALSXP,ends.size())); 
+    
+    double* csa=REAL(cs_R);
+    int i=0;
+    for(vector<double>::const_iterator ci=starts.begin(); ci!=starts.end(); ++ci) {
+      csa[i]=*ci; i++;
+    }
+
+    csa=REAL(ce_R);
+    i=0;
+    for(vector<double>::const_iterator ci=ends.begin(); ci!=ends.end(); ++ci) {
+      csa[i]=*ci; i++;
+    }
+    
+    if(max_val!=0) {
+      PROTECT(cv_R=allocVector(REALSXP,values.size())); 
+      csa=REAL(cv_R);
+      i=0;
+      for(vector<double>::const_iterator ci=values.begin(); ci!=values.end(); ++ci) {
+	csa[i]=*ci; i++;
+      }
+    }
+
+    SEXP ans_R, names_R;
+    if(max_val!=0) {
+      PROTECT(names_R = allocVector(STRSXP, 3));
+      SET_STRING_ELT(names_R, 0, mkChar("s"));
+      SET_STRING_ELT(names_R, 1, mkChar("e"));
+      SET_STRING_ELT(names_R, 2, mkChar("v"));
+
+      PROTECT(ans_R = allocVector(VECSXP, 3));
+      SET_VECTOR_ELT(ans_R, 0, cs_R);
+      SET_VECTOR_ELT(ans_R, 1, ce_R);
+      SET_VECTOR_ELT(ans_R, 2, cv_R);
+    } else {
+      PROTECT(names_R = allocVector(STRSXP, 2));
+      SET_STRING_ELT(names_R, 0, mkChar("s"));
+      SET_STRING_ELT(names_R, 1, mkChar("e"));
+
+      PROTECT(ans_R = allocVector(VECSXP, 2));
+      SET_VECTOR_ELT(ans_R, 0, cs_R);
+      SET_VECTOR_ELT(ans_R, 1, ce_R);
+    }
+    
+    setAttrib(ans_R, R_NamesSymbol, names_R);
+    
+    if(max_val!=0) {
+      UNPROTECT(5);
+    } else {
+      UNPROTECT(4);
+    }
+    return(ans_R);
+  }
+
+}
+