/*===========================================================================
  C-equivalent of R-function lagderiv
  =========================================================================== */
SEXP getLagDeriv(SEXP T, SEXP nr)
{
    SEXP value;
    int i, ilen, interval;
    double t;

    ilen = LENGTH(nr);
    if (initialisehist == 0)
        error("pastgradient can only be called from 'func' or 'res' when triggered by appropriate integrator.");
    if (!isNumeric(T)) error("'t' should be numeric");

    t = *NUMERIC_POINTER(T);
    interval = findHistInt (t);

    if ((ilen ==1) && (INTEGER(nr)[0] == 0)) {
        PROTECT(value=NEW_NUMERIC(n_eq));
        for(i=0; i<n_eq; i++) {
            NUMERIC_POINTER(value)[i] = past(i, interval, t, 2);
        }
    } else {
        PROTECT(value=NEW_NUMERIC(ilen));
        for(i=0; i<ilen; i++) {
            NUMERIC_POINTER(value)[i] = past(INTEGER(nr)[i]-1, interval, t, 2);
        }
    }
    UNPROTECT(1);
    return(value);
}

STGM::CBoolSphereSystem * InitSphereSystem(SEXP R_param, SEXP R_cond) {
  SEXP R_box;
  PROTECT( R_box  = getListElement( R_cond, "box"));
  double *boxX = NUMERIC_POINTER( getListElement( R_box, "xrange"));
  double *boxY = NUMERIC_POINTER( getListElement( R_box, "yrange"));
  double *boxZ = NUMERIC_POINTER( getListElement( R_box, "zrange"));

  double lam   = asReal(AS_NUMERIC( getListElement( R_param, "lam")));

  /*  print level */
  PL = asInteger(getListElement( R_cond,"pl"));

  /* simulation box */
  STGM::CBox3 box(boxX,boxY,boxZ);

  /* set up sphere system */
  STGM::CBoolSphereSystem *sp = (STGM::CBoolSphereSystem*)Calloc(1,STGM::CBoolSphereSystem);

  try {
      new(sp)STGM::CBoolSphereSystem(box,lam);
  } catch(...) {
      error(_("InitSpheroidSystem(): Memory allocation error for sphere system."));
  }

  UNPROTECT(1);
  return sp;
}

//////////////////////////////////////////////////
// addXAxis - add X axis information
unsigned int addXAxis(SEXP data, SEXP dataNames, unsigned int j, TH1* hist)
{
  int n = hist->GetNbinsX();
  TAxis* axis = hist->GetXaxis();
  
  // Determine breaks--
  // Add to list
  SEXP breaks = addNumericVector(data, dataNames, j++, n+1, "breaks");

  // Get information
  for ( unsigned int i=0; i<n; ++i ) {
    NUMERIC_POINTER(breaks)[i] = axis->GetBinLowEdge(i+1);
  }

  // Add the high edge
  NUMERIC_POINTER(breaks)[n] = axis->GetBinUpEdge(n);

  // Determine mids--
  SEXP mids = addNumericVector(data, dataNames, j++, n, "mids");

  // Get information
  for ( unsigned int i=0; i<n; ++i ) {
    NUMERIC_POINTER(mids)[i] = axis->GetBinCenter(i+1);
  }

  // Get name of axis
  SEXP xname = addCharVector(data, dataNames, j++, 1, "xname");
  SET_STRING_ELT( xname, 0, mkChar( axis->GetTitle() ) );
  
  // Done
  return j;
}

SEXP R_RngStreams_GetData (SEXP R_stream)
     /*----------------------------------------------------------------------*/
     /* Get data structure of Stream object.                                 */
     /* (For the name of the Stream object use R_RngStreams_GetName() ).     */
     /*                                                                      */
     /* parameters:                                                          */
     /*   R_stream ... (pointer) ... pointer the Stream object               */
     /*                                                                      */
     /* return:                                                              */
     /*   data (double[20])                                                  */
     /*----------------------------------------------------------------------*/
{
  RngStream stream;
  SEXP R_stream_data;

  /* check pointer */
  CHECK_STREAM_PTR(R_stream);

  /* Extract pointer to generator */
  stream = R_ExternalPtrAddr(R_stream);
  CHECK_NULL(stream);

  PROTECT(R_stream_data = NEW_NUMERIC(20));
  memcpy(NUMERIC_POINTER(R_stream_data)   , stream->Cg, 6*sizeof(double));
  memcpy(NUMERIC_POINTER(R_stream_data)+ 6, stream->Bg, 6*sizeof(double));
  memcpy(NUMERIC_POINTER(R_stream_data)+12, stream->Ig, 6*sizeof(double));
  NUMERIC_POINTER(R_stream_data)[18] = (double) stream->Anti;
  NUMERIC_POINTER(R_stream_data)[19] = (double) stream->IncPrec;

  UNPROTECT(1);
  
  /* return data to R */
  return R_stream_data;

} /* end of R_RngStreams_GetData() */

SEXP point_in_polygon(SEXP px, SEXP py, SEXP polx, SEXP poly) {
	int i;
	PLOT_POINT p;
	POLYGON pol;
	SEXP ret;

	S_EVALUATOR
	pol.lines = LENGTH(polx); /* check later that first == last */
	pol.p = (PLOT_POINT *) Calloc(pol.lines, PLOT_POINT); /* Calloc does error handling */
	for (i = 0; i < LENGTH(polx); i++) {
		pol.p[i].x = NUMERIC_POINTER(polx)[i];
		pol.p[i].y = NUMERIC_POINTER(poly)[i];
	}
    pol.close = (pol.p[0].x == pol.p[pol.lines - 1].x && 
			pol.p[0].y == pol.p[pol.lines - 1].y);
	setup_poly_minmax(&pol);

	PROTECT(ret = NEW_INTEGER(LENGTH(px)));
	for (i = 0; i < LENGTH(px); i++) {
		p.x = NUMERIC_POINTER(px)[i];
		p.y = NUMERIC_POINTER(py)[i];
		if ((p.x > pol.mbr.min.x) & (p.x <= pol.mbr.max.y) & (p.y > pol.mbr.min.y) & (p.y <= pol.mbr.max.y)) {
			INTEGER_POINTER(ret)[i] = InPoly(p, &pol);
		} else {
			INTEGER_POINTER(ret)[i] = 0;
		}
	}
	Free(pol.p);
	UNPROTECT(1);
	return(ret);
}

SEXP seqlib_tm_santa_lucia(SEXP sequences,SEXP ct)
{
	int vlen,i,wg=0;
	sequence_tp*ms;
	double Ct;
	SEXP res;
	if(!isString(sequences))
	   error("sequence must have character type");
	if (!isReal(ct) || length(ct) != 1)
	   error("ct value must be single real");
	vlen = length(sequences);
	Ct = REAL(ct)[0];
	PROTECT(res = NEW_NUMERIC(vlen));
	for (i=0; i< vlen; i++)
	{
		ms = sequence_from_string(CHAR(STRING_ELT(sequences,i)));
		if (sequence_conv_to_acgt_only(ms))
			NUMERIC_POINTER(res)[i] = sequence_melt_nn_SantaLucia(ms,Ct);
		else 
		{
			NUMERIC_POINTER(res)[i] = NA_REAL;
			if (!wg) 
			{
				warning("Non-determined nucleotides in sequences");
				wg = 1;
			}
		}
		free(ms);
	}
	UNPROTECT(1);
	return res;
}

// version for .Call, faster because nothing is duplicated
SEXP jarowinklerCALL(SEXP str1EXP, SEXP str2EXP, SEXP W_1EXP, SEXP W_2EXP, 
  SEXP W_tEXP, SEXP rEXP)
{
  const char *str_1, *str_2;
  double *W_1, *W_2, *W_t, *r, *ans;
  int length_1, length_2, maxlen;
  SEXP ret;
  W_1 = NUMERIC_POINTER(W_1EXP);
  W_2 = NUMERIC_POINTER(W_2EXP);
  W_t = NUMERIC_POINTER(W_tEXP);
  r = NUMERIC_POINTER(rEXP);
  length_1 = LENGTH(str1EXP);
  length_2 = LENGTH(str2EXP);
  maxlen = length_1 > length_2 ? length_1 : length_2; 
  PROTECT(ret = NEW_NUMERIC(maxlen));
  ans = NUMERIC_POINTER(ret);
  for (int str_ind=0; str_ind < maxlen; str_ind++)
  {
    str_1=CHAR(STRING_ELT(str1EXP, str_ind % length_1));
    str_2=CHAR(STRING_ELT(str2EXP, str_ind % length_2));
    ans[str_ind]=jarowinkler_core(str_1, str_2, *W_1, *W_2, *W_t, *r);
  }
  UNPROTECT(1);
  return(ret);
}

// Return closest point to given distance within geometry.
// 'spgeom' must be a LineString
SEXP rgeos_interpolate(SEXP env, SEXP spgeom, SEXP d, SEXP normalized) {

    GEOSContextHandle_t GEOShandle = getContextHandle(env);
    GEOSGeom geom = rgeos_convert_R2geos(env, spgeom);

    GEOSGeom res_geos;
    double dist;

    int nlines = length(GET_SLOT(spgeom, install("lines")));
    if (nlines < 1) {
        error("rgeos_project: invalid number of lines");
    }

    int n = LENGTH(d);
    if (n < 1) {
        error("rgeos_interpolate: invalid number of requested points");
    }

    int pc = 0;
    SEXP crd;
    PROTECT(crd = NEW_NUMERIC(n*2)); pc++;

    double x;
    double y;
    SEXP ans;

    // select interpolation function (normalized/unnormalized)
    GEOSGeometry GEOS_DLL *(*interp_fun)(GEOSContextHandle_t,
                                         const GEOSGeometry*,
                                         double);

    if (LOGICAL_POINTER(normalized)[0]) {
        interp_fun = &GEOSInterpolateNormalized_r;
    } else {
        interp_fun = &GEOSInterpolate_r;
    }

    // interpolate points and store result in coord matrix
    for (int i = 0; i < n; i++) {

        dist = NUMERIC_POINTER(d)[i];

        res_geos = (*interp_fun)(GEOShandle, geom, dist);

        rgeos_Pt2xy(env, res_geos, &x, &y);

        NUMERIC_POINTER(crd)[i] = x;
        NUMERIC_POINTER(crd)[n+i] = y;
    }

    GEOSGeom_destroy_r(GEOShandle, geom);
    GEOSGeom_destroy_r(GEOShandle, res_geos);

    // return coordinates as matrix
    PROTECT(ans = rgeos_formatcrdMat(crd, n)); pc++;

    UNPROTECT(pc);
    return(ans);
}

/**
* @brief Summarizes a list of vectors into a list of binned vectors of equal length. Each vector bin summarizes an approximately equal amount of values.
*
* @param method Charater array defining the method to be used for binning. Can be 'mean' 'median' or 'max'
* @param score_list List with numeric vectors
* @param window_size Window width of the vectors that will be returned
* @return List with updated vectors
* @details Walks through the vectors and calls shrink or expand to set vectors to equal widths
* @note Nothing
* @todo Nothing
*/
SEXP approx_window(SEXP window_count, SEXP score_list, SEXP method) {
    const char *methodn = STRING_VALUE(method);
    const int wsize=INTEGER_VALUE(window_count);

    SEXP lnames = getAttrib(score_list, R_NamesSymbol);
    SEXP ori_vec,new_vec,out_names,out_list;
    int elcount=0,elements=LENGTH(lnames),upc=0,olen;
    signal(SIGINT,SIG_DFL);
    PROTECT(lnames = AS_CHARACTER(lnames));
    upc++;
    PROTECT(out_list = allocVector(VECSXP, elements));
    upc++;
    PROTECT(out_names = allocVector(STRSXP,elements));
    upc++;

    //Select proper call back
    double (*summarizep)(int *,int,double *);
    if(!strcmp(methodn,"mean")) {
        summarizep=mean_dble;
    } else if(!strcmp(methodn,"median")) {
        summarizep=median_dble;
    } else if(!strcmp(methodn,"max")) {
        summarizep=vect_max_dble;
    } else {
        error("%s not known",methodn);
        goto FINALIZE;
    }


    for(; elcount<elements; ++elcount) {
        PROTECT(ori_vec=AS_NUMERIC(VECTOR_ELT(score_list, elcount)));
        PROTECT(new_vec = NEW_NUMERIC(wsize));
        olen=LENGTH(ori_vec);
        double *ori_vecp= NUMERIC_POINTER(ori_vec);
        double *new_vecp= NUMERIC_POINTER(new_vec);
        SET_STRING_ELT(out_names,elcount,mkChar(CHAR(STRING_ELT(lnames, elcount))));
        if(olen>wsize) {
            shrink_dble(ori_vecp,new_vecp,olen,wsize,summarizep);
            SET_VECTOR_ELT(out_list, elcount, new_vec);
        } else if(olen<wsize) {
            expand_dble(ori_vecp,new_vecp,olen,wsize);
            SET_VECTOR_ELT(out_list, elcount, new_vec);
        } else {
            SET_VECTOR_ELT(out_list, elcount, ori_vec);
        }

        UNPROTECT(2);
    }
    setAttrib(out_list, R_NamesSymbol, out_names);

FINALIZE:
    UNPROTECT(upc);
    return(out_list);
}

  SEXP fastcluster_correlation_distances(SEXP matrix_, SEXP const nrow_, SEXP const ncol_, SEXP const type_) {
    SEXP distance = NULL; // return value

    try{
			PROTECT(type_);
			const int type = *INTEGER_POINTER(type_);
			UNPROTECT(1); // type_

      PROTECT(nrow_);
      if (!IS_INTEGER(nrow_) || LENGTH(nrow_)!=1) Rf_error("'nrow' must be a single integer.");
      const int nrow = *INTEGER_POINTER(nrow_);
      UNPROTECT(1); // nrow_
      if (nrow<2) Rf_error("nrow must be at least 2.");

      PROTECT(ncol_);
      if (!IS_INTEGER(ncol_) || LENGTH(ncol_)!=1) Rf_error("'ncol' must be a single integer.");
      const int ncol = *INTEGER_POINTER(ncol_);
      UNPROTECT(1); // ncol_
      if (ncol<2) Rf_error("ncol must be at least 2.");

      const std::ptrdiff_t N = static_cast<std::ptrdiff_t>(nrow*ncol);

      PROTECT(matrix_ = AS_NUMERIC(matrix_));
      if (LENGTH(matrix_)!=N) Rf_error("Improperly specified matrix dimensions.");
      const double * const matrix = NUMERIC_POINTER(matrix_);

			// R defaults to by-column comparisons
			const std::ptrdiff_t dsize = static_cast<std::ptrdiff_t>((ncol)*(ncol-1)/2);

      PROTECT(distance = NEW_NUMERIC(dsize));
      double * const d = NUMERIC_POINTER(distance);

			if(type==2) pearson_distances_pairwise_complete_obs_variant(d, matrix, nrow, ncol);
			else if(type==3) spearman_distances_pairwise_complete_obs(d, matrix, nrow, ncol);
			else pearson_distances_pairwise_complete_obs(d, matrix, nrow, ncol);

      UNPROTECT(2); // matrix_ and distance

    } // try
    catch (const std::bad_alloc&) {
      Rf_error( "Memory overflow.");
    }
    catch(const std::exception& e){
      Rf_error( e.what() );
    }
    catch(const nan_error&){
      Rf_error("NaN dissimilarity value.");
    }
    catch(...){
      Rf_error( "C++ exception (unknown reason)." );
    }

    return distance;
  }

SEXP MRF_Stat(SEXP _crf, SEXP _instances)
{
	CRF crf(_crf);

	int nInstances = INTEGER_POINTER(GET_DIM(_instances))[0];
	int nPar = INTEGER_POINTER(AS_INTEGER(GetVar(_crf, "n.par")))[0];

	PROTECT(_instances = AS_NUMERIC(_instances));
	double *instances = NUMERIC_POINTER(_instances);

	SEXP _nodePar;
	PROTECT(_nodePar = AS_INTEGER(GetVar(_crf, "node.par")));
	int *nodePar = INTEGER_POINTER(_nodePar);

	SEXP _edgePar = GetVar(_crf, "edge.par");
	int **edgePar = (int **) R_alloc(crf.nEdges, sizeof(int *));
	SEXP _edgeParI, _temp;
	PROTECT(_edgeParI = NEW_LIST(crf.nEdges));
	for (int i = 0; i < crf.nEdges; i++)
	{
		SET_VECTOR_ELT(_edgeParI, i, _temp = AS_INTEGER(GetListElement(_edgePar, i)));
		edgePar[i] = INTEGER_POINTER(_temp);
	}

	SEXP _stat;
	PROTECT(_stat = NEW_NUMERIC(nPar));
	double *stat = NUMERIC_POINTER(_stat);
	SetValues(_stat, stat, 0.0);

	int *y = (int *) R_allocVector<int>(crf.nNodes);

	for (int n = 0; n < nInstances; n++)
	{
		for (int i = 0; i < crf.nNodes; i++)
		{
			y[i] = instances[n + nInstances * i] - 1;
			int p = nodePar[i + crf.nNodes * y[i]] - 1;
			if (p >= 0 && p < nPar)
				stat[p]++;
		}

		for (int i = 0; i < crf.nEdges; i++)
		{
			int p = edgePar[i][y[crf.EdgesBegin(i)] + crf.nStates[crf.EdgesBegin(i)] * y[crf.EdgesEnd(i)]] - 1;
			if (p >= 0 && p < nPar)
				stat[p]++;
		}
	}

	UNPROTECT(4);

	return(_stat);
}

// Return distance of points 'spppoints' projected on 'spgeom' from origin
// of 'spgeom'. Geometry 'spgeom' must be a lineal geometry
SEXP rgeos_project(SEXP env, SEXP spgeom, SEXP sppoint, SEXP normalized) {

    GEOSContextHandle_t GEOShandle = getContextHandle(env);
    GEOSGeom geom = rgeos_convert_R2geos(env, spgeom);

    SEXP crds = GET_SLOT(sppoint, install("coords"));
    SEXP dim = getAttrib(crds, install("dim"));

    int nlines = length(GET_SLOT(spgeom, install("lines")));
    if (nlines < 1) {
        error("rgeos_project: invalid number of lines");
    }

    int n = INTEGER_POINTER(dim)[0];
    if (n < 1) {
        error("rgeos_project: invalid number of points");
    }

    int pc = 0;
    SEXP ans;
    PROTECT(ans = NEW_NUMERIC(n)); pc++;

    GEOSGeom p;

    // select projection function (normalized/unnormalized)
    double GEOS_DLL (*proj_fun)(GEOSContextHandle_t,
                                const GEOSGeometry*,
                                const GEOSGeometry*);

    if (LOGICAL_POINTER(normalized)[0]) {
        proj_fun = &GEOSProjectNormalized_r;
    } else {
        proj_fun = &GEOSProject_r;
    }

    // project points to line geometry
    for (int i = 0; i < n; i++) {

        p = rgeos_xy2Pt(env,
                        NUMERIC_POINTER(crds)[i],
                        NUMERIC_POINTER(crds)[i+n]);

        NUMERIC_POINTER(ans)[i] = (*proj_fun)(GEOShandle, geom, p);
    }

    GEOSGeom_destroy_r(GEOShandle, geom);
    GEOSGeom_destroy_r(GEOShandle, p);

    UNPROTECT(pc);

    return(ans);
}

SEXP lmin22(SEXP nb, SEXP y, SEXP cy, SEXP card, SEXP beta) {
    int i, j, k, nswitch=0, n=length(card), pc=0;
    SEXP ans;
    double t1, t2, ytemp, yhat;
    double *Y, *CY, *B;

    Y = (double *) R_alloc((size_t) n, sizeof(double));
    CY = (double *) R_alloc((size_t) n, sizeof(double));
    B = (double *) R_alloc((size_t) length(beta), sizeof(double));

    for (i=0; i<n; i++) {
        Y[i] = NUMERIC_POINTER(y)[i];
        CY[i] = NUMERIC_POINTER(cy)[i];
    }
    for (i=0; i<length(beta); i++) {
        B[i] = NUMERIC_POINTER(beta)[i];
    }

    PROTECT(ans = NEW_LIST(2)); pc++;
    SET_VECTOR_ELT(ans, 0, NEW_NUMERIC(n));
    SET_VECTOR_ELT(ans, 1, NEW_INTEGER(1));

    for (i=0; i<n; i++) {
      if (INTEGER_POINTER(card)[i] > 0) {
        t1 = fabs(Y[i] - CY[i]);
        yhat = B[0] + B[1]*CY[i];
        t2 = fabs(yhat - CY[i]);
        for (j=0; j<INTEGER_POINTER(card)[i]; j++) {
            k = INTEGER_POINTER(VECTOR_ELT(nb, i))[j]-ROFFSET;
            t1 = t1 + fabs(Y[k] - CY[k]);
            t2 = t2 + fabs(Y[k] - (CY[k] - Y[i] + B[0] + B[1]*CY[i]));
        }
        if (t1 <= t2) {
            nswitch++;
            ytemp = Y[i];
            Y[i] = yhat;
            for (j=0; j<INTEGER_POINTER(card)[i]; j++) {
                k = INTEGER_POINTER(VECTOR_ELT(nb, i))[j]-ROFFSET;
                CY[k] = CY[k] - ytemp + Y[i];
            }
        }
      }
    }

    for (i=0; i<n; i++) {
        NUMERIC_POINTER(VECTOR_ELT(ans, 0))[i] = Y[i];
    }

    INTEGER_POINTER(VECTOR_ELT(ans, 1))[0] = nswitch;
    UNPROTECT(pc); /* ans */
    return(ans);
}

SEXP 
R_cv_svd_wold (SEXP xx, SEXP kk, SEXP maxrankmaxrank, SEXP toltol, 
               SEXP maxitermaxiter, SEXP setssets)
{
    bcv_error_t err = 0;
    bcv_index_t m, n, i, k, maxiter, maxrank;
    bcv_svd_wold_t *wold = NULL;
    double tol, *msep;
    SEXP msepmsep, dimdim;

    m       = INTEGER (getAttrib (xx, R_DimSymbol))[0];
    n       = INTEGER (getAttrib (xx, R_DimSymbol))[1];
    k       = asInteger (kk);
    maxrank = asInteger (maxrankmaxrank);
    tol     = asReal (toltol);
    maxiter = asInteger (maxitermaxiter);

    PROTECT (msepmsep = allocVector (REALSXP, (maxrank + 1) * k));
    PROTECT (dimdim   = allocVector (INTSXP, 2));
    INTEGER (dimdim) [0] = maxrank + 1;
    INTEGER (dimdim) [1] = k;
    setAttrib (msepmsep, R_DimSymbol, dimdim);
    msep = NUMERIC_POINTER (msepmsep);

    bcv_matrix_t x       = { m, n, NUMERIC_POINTER (xx), BCV_MAX (m,1) };
    bcv_partition_t part = { m*n, k, INTEGER_POINTER (setssets) };

    wold = bcv_svd_wold_alloc (m*n, m, n);
    if (!wold)
        error ("could not allocate enough memory for Wold "
               " cross-validation of a %d-by-%d matrix", m, n);

    bcv_svd_wold_init (wold, &x, &part);
    
    for (i = 0; i < k; i++)
    {
        R_CheckUserInterrupt ();
        err = bcv_svd_wold_get_msep (wold, i, tol, maxiter, msep, maxrank);
            
        if (err)
            error ("the SVD algorithm did not converge for the (%d)"
                   " holdout", i);
        
        msep += maxrank + 1;
    }
    
    bcv_svd_wold_free (wold);
    
    UNPROTECT (2);
    return msepmsep;
}

SEXP spOverlap(SEXP bbbi, SEXP bbbj) {

	int pc=0,overlap=1;
	double bbi[4], bbj[4];
	SEXP ans;

	PROTECT(ans = NEW_INTEGER(1)); pc++;
	bbi[0] = NUMERIC_POINTER(bbbi)[0];
	bbi[1] = NUMERIC_POINTER(bbbi)[1];
	bbi[2] = NUMERIC_POINTER(bbbi)[2];
	bbi[3] = NUMERIC_POINTER(bbbi)[3];
	bbj[0] = NUMERIC_POINTER(bbbj)[0];
	bbj[1] = NUMERIC_POINTER(bbbj)[1];
	bbj[2] = NUMERIC_POINTER(bbbj)[2];
	bbj[3] = NUMERIC_POINTER(bbbj)[3];

        if ((bbi[0]>bbj[2]) | (bbi[1]>bbj[3]) | 
		(bbi[2]<bbj[0]) | (bbi[3]<bbj[1]) ) {
		overlap=0;
	}

	INTEGER_POINTER(ans)[0] = overlap;		
	UNPROTECT(pc); /* ans */
	return(ans);
}

//////////////////////////////////////////////////
// addXContents
unsigned int addXContents(SEXP data, SEXP dataNames, unsigned int j, TH1* hist)
{
  int n = hist->GetNbinsX();
  
  // Determine counts and uncert--
  // Add to list
  SEXP counts =  addNumericVector(data, dataNames, j++, n, "counts");
  SEXP uncert =  addNumericVector(data, dataNames, j++, n, "uncert");
  
  for ( unsigned int i=0; i<n; ++i ) {
    NUMERIC_POINTER(counts)[i] = hist->GetBinContent(i+1);
    NUMERIC_POINTER(uncert)[i] = hist->GetBinError(i+1);
  }

  // Set under and over flows
  SEXP uof = addNumericVector(data, dataNames, j++, 2, "underOverFlows");
  NUMERIC_POINTER(uof)[0] = hist->GetBinContent(0);
  NUMERIC_POINTER(uof)[1] = hist->GetBinContent(n+1);

  // Get mean
  SEXP mean = addNumericVector(data, dataNames, j++, 2, "mean");
  NUMERIC_POINTER(mean)[0] = hist->GetMean(1);
  NUMERIC_POINTER(mean)[1] = hist->GetMean(11);

  // Get rms
  SEXP rms = addNumericVector(data, dataNames, j++, 2, "rms");
  NUMERIC_POINTER(rms)[0] = hist->GetRMS(1);
  NUMERIC_POINTER(rms)[1] = hist->GetRMS(11);
  
  return j;
}

SEXP mom_calc_int2(SEXP is, SEXP m, SEXP nb, SEXP weights, SEXP card) {
    SEXP Omega;
    int hm = INTEGER_POINTER(m)[0];
    int n = length(card);
    double *eta, *zeta, *omega, sum, res;
    int i, ii, j, k1, k2, k3;
    int iis = length(is);

    omega = (double *) R_alloc((size_t) hm, sizeof(double));
    eta = (double *) R_alloc((size_t) n, sizeof(double));
    zeta = (double *) R_alloc((size_t) n, sizeof(double));
    for (j=0; j<hm; j++) omega[j] = 0.0;

    for (ii=0; ii<iis; ii++) {
        R_CheckUserInterrupt();
        i = INTEGER_POINTER(is)[ii]-ROFFSET;
        for (j=0; j<n; j++) eta[j] = 0.0;
        eta[i] = 1.0;
        for (j=1; j<hm; j=j+2) {
            for (k1=0; k1<n; k1++) {
                k3 = INTEGER_POINTER(card)[k1];
                if (k3 == 0) {
                    zeta[k1] = 0.0;
                } else {
                    sum = 0.0;
                    for (k2=0; k2<k3; k2++) {
                        sum += eta[INTEGER_POINTER(VECTOR_ELT(nb, k1))[k2]
                            - ROFFSET] * NUMERIC_POINTER(VECTOR_ELT(weights,
                            k1))[k2];
                    }
                    zeta[k1] = sum;
                }
            }
            res = F77_CALL(ddot)(&n, zeta, &c__1, eta, &c__1);
            if (R_FINITE(res)) omega[(j-1)] += res;
            else error("non-finite dot product %d, %d", i, j);
            res = F77_CALL(ddot)(&n, zeta, &c__1, zeta, &c__1);
            if (R_FINITE(res)) omega[j] += res;
            else error("non-finite dot product %d, %d", i, j);
            for (k1=0; k1<n; k1++) eta[k1] = zeta[k1];
        }
    }

    PROTECT(Omega = NEW_NUMERIC(hm));
    for (j=0; j<hm; j++) NUMERIC_POINTER(Omega)[j] = omega[j];

    UNPROTECT(1);
    return(Omega);
}

SEXP R_RngStreams_Sample (SEXP R_stream, SEXP R_size)
     /*----------------------------------------------------------------------*/
     /* Sample from Stream object.                                           */
     /*                                                                      */
     /* parameters:                                                          */
     /*   R_stream ... (pointer) ... pointer the Stream object               */
     /*   R_size   ... (int)     ... sample size                             */
     /*                                                                      */
     /* return:                                                              */
     /*   pointer to stream object                                           */
     /*----------------------------------------------------------------------*/
{
  RngStream stream;
  int n = INTEGER(R_size)[0];
  int i;
  SEXP R_sample;

  /* check pointer */
  CHECK_STREAM_PTR(R_stream);

  /* Extract pointer to generator */
  stream = R_ExternalPtrAddr(R_stream);
  CHECK_NULL(stream);

  /* generate random sample of size n */
  PROTECT(R_sample = NEW_NUMERIC(n));
  for (i=0; i<n; i++)
    NUMERIC_POINTER(R_sample)[i] = RngStream_RandU01(stream);
  UNPROTECT(1);

  /* return sample to R */
  return R_sample;
 
} /* end of R_RngStreams_sample() */

SEXP match3bytes(SEXP buf, SEXP m1, SEXP m2, SEXP m3)
{
  int i, j, n, n_match;
  double *resp;
  unsigned char *bufp, *m1p, *m2p, *m3p;
  SEXP res;
  PROTECT(buf = AS_RAW(buf));
  PROTECT(m1 = AS_RAW(m1));
  PROTECT(m2 = AS_RAW(m2));
  PROTECT(m3 = AS_RAW(m3));
  bufp = RAW_POINTER(buf);
  m1p = RAW_POINTER(m1);
  m2p = RAW_POINTER(m2);
  m3p = RAW_POINTER(m3);
  n = LENGTH(buf);
  n_match = 0;
  for (i = 0; i < n - 2; i++) {
    if (bufp[i] == *m1p && bufp[i + 1] == *m2p && bufp[i + 2] == *m3p) {
      n_match++;
      ++i;   /* skip */
      ++i;   /* skip */
    }
  }
  PROTECT(res = NEW_NUMERIC(n_match));
  resp = NUMERIC_POINTER(res);
  j = 0;
  for (i = 0; i < n - 2; i++) {
    if (j <= n_match && bufp[i] == *m1p && bufp[i + 1] == *m2p && bufp[i + 2] == *m3p) {
      resp[j++] = i + 1; /* the 1 is to offset from C to R */
    }
  }
  UNPROTECT(5);
  return(res);
}

SEXP seqlib_dg_selfhyb(SEXP sequences,SEXP t,SEXP naC,SEXP mgC,SEXP acidtype)
{
	int vlen,i,atype;
	double ct,nac,mgc,rv;
	sequence_tp *ms;
	SEXP res;
	if(!isString(sequences))
	   error("sequence must have character type");
	vlen = length(sequences);
	PROTECT(res = NEW_NUMERIC(vlen));

	if (!isReal(t) || length(t) != 1)
              error("t value must be single real");
        ct = REAL(t)[0];
	if (!isReal(naC) || length(naC) != 1)
              error("naC value must be single real");
        nac = REAL(naC)[0];
	if (!isReal(mgC) || length(mgC) != 1)
              error("mgC value must be single real");
        mgc = REAL(mgC)[0];
	if (!isInteger(acidtype) || length(acidtype) != 1)
              error("acidtype value must be single int");
        atype = INTEGER(acidtype)[0];
	for (i=0; i< vlen; i++)
	{
		ms = sequence_from_string(CHAR(STRING_ELT(sequences,i)));
		rv = sequence_hybrid_ss_min(ms,ct,nac,mgc,atype);
		NUMERIC_POINTER(res)[i] = rv;
		free(ms);
	}
	UNPROTECT(1);
	return res;
}