SEXP call_gambim(SEXP y, SEXP times, SEXP derivfunc, SEXP parms, SEXP rtol,
		SEXP atol, SEXP rho, SEXP Tcrit, SEXP jacfunc, SEXP initfunc,
    SEXP verbose, SEXP LRW, SEXP rWork, SEXP iWork, SEXP jT,
    SEXP nOut, SEXP Nrmas, SEXP masfunc, SEXP ML, SEXP MU, SEXP Hini,
    SEXP Rpar, SEXP Ipar, SEXP flist, SEXP Type)

{
/******************************************************************************/
/******                   DECLARATION SECTION                            ******/
/******************************************************************************/

  int  j, nt, latol, lrtol, imas, mlmas, mumas, type;
  int  isForcing, runOK;
  double *Atol, *Rtol, hini;
  int itol, ijac, ml, mu, iout, idid, liw, lrw, sum;

  /* pointers to functions passed to FORTRAN */
  C_jac_func_type_gb    *jac_func_gb = NULL;
  C_solout_type         *solout = NULL;
  C_solout_type_bim     *solout_bim = NULL;
  C_mas_type            *mas_func = NULL;

/******************************************************************************/
/******                         STATEMENTS                               ******/
/******************************************************************************/

/*                      #### initialisation ####                              */
  init_N_Protect();

  type  = INTEGER(Type)[0];     /* jacobian type */
  ijac  = INTEGER(jT)[0];       /* jacobian type */
  n_eq = LENGTH(y);             /* number of equations */
  nt   = LENGTH(times);         /* number of output times */
  maxt = nt;

  tt = (double *) R_alloc(nt, sizeof(double));
  for (j = 0; j < nt; j++) tt[j] = REAL(times)[j];

//  mflag = INTEGER(verbose)[0];

  /* is function a dll ?*/
  isDll = inherits(derivfunc, "NativeSymbol");

  /* initialise output ... */
  initOutC(isDll, n_eq, nOut, Rpar, Ipar);

  /* copies of variables that will be changed in the FORTRAN subroutine */
  xytmp = (double *) R_alloc(n_eq, sizeof(double));
  for (j = 0; j < n_eq; j++) xytmp[j] = REAL(y)[j];

  ytmp = (double *) R_alloc(n_eq, sizeof(double));
  for (j = 0; j < n_eq; j++) ytmp[j] = 0.;

  latol = LENGTH(atol);
  Atol = (double *) R_alloc((int) latol, sizeof(double));
  for (j = 0; j < latol; j++) Atol[j] = REAL(atol)[j];

  lrtol = LENGTH(rtol);
  Rtol = (double *) R_alloc((int) lrtol, sizeof(double));
  for (j = 0; j < lrtol; j++) Rtol[j] = REAL(rtol)[j];

  /* tolerance specifications */
  if (latol == 1 ) itol = 0;
  else             itol = 1;

  /* mass matrix */
  imas  = INTEGER(Nrmas)[0];
  mlmas = INTEGER(Nrmas)[1];
  mumas = INTEGER(Nrmas)[2];

  /* work vectors */
  if (type == 1) {
    liw = 27;
    lrw = 21;
  } else  {       //  if (type == 2)
    liw = n_eq + 40;
    lrw = INTEGER(LRW)[0];
 }

  iwork = (int *) R_alloc(liw, sizeof(int));
  for (j=0; j<LENGTH(iWork); j++) iwork[j] = INTEGER(iWork)[j];
  for (j=LENGTH(iWork); j<liw; j++) iwork[j] = 0;

  rwork = (double *) R_alloc(lrw, sizeof(double));
  for (j=0; j<length(rWork); j++) rwork[j] = REAL(rWork)[j];
  for (j=length(rWork); j<lrw; j++) rwork[j] = 0.;

  ml = INTEGER(ML)[0];
  mu = INTEGER(MU)[0];
  hini = REAL(Hini)[0];

  /* initialise global R-variables...  */
  initglobals (nt);

  /* Initialization of Parameters, Forcings (DLL) */
  initParms(initfunc, parms);
  isForcing = initForcings(flist);

  if (nout > 0 ) {
     xdytmp= (double *) R_alloc(n_eq, sizeof(double));
//.........这里部分代码省略.........

SEXP do_init_state (SEXP object, SEXP params, SEXP t0, SEXP nsim, SEXP gnsi)
{
  int nprotect = 0;
  SEXP Pnames, Snames;
  SEXP x = R_NilValue;
  int *dim;
  int npar, nrep, nvar, ns;
  int definit;
  int xdim[2];
  const char *dimnms[2] = {"variable","rep"};

  ns = *(INTEGER(AS_INTEGER(nsim)));
  PROTECT(params = as_matrix(params)); nprotect++;
  PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params))); nprotect++;
  dim = INTEGER(GET_DIM(params));
  npar = dim[0]; nrep = dim[1]; 

  if (ns % nrep != 0) 
    errorcall(R_NilValue,"in 'init.state': number of desired state-vectors 'nsim' is not a multiple of ncol('params')");

  definit = *(INTEGER(GET_SLOT(object,install("default.init"))));

  if (definit) {		// default initializer

    SEXP fcall, pat, repl, val, ivpnames, statenames;
    int *pidx, j, k;
    double *xp, *pp;
  
    PROTECT(pat = NEW_CHARACTER(1)); nprotect++;
    SET_STRING_ELT(pat,0,mkChar("\\.0$"));
    PROTECT(repl = NEW_CHARACTER(1)); nprotect++;
    SET_STRING_ELT(repl,0,mkChar(""));
    PROTECT(val = NEW_LOGICAL(1)); nprotect++;
    *(INTEGER(val)) = 1;
    
    // extract names of IVPs
    PROTECT(fcall = LCONS(val,R_NilValue)); nprotect++;
    SET_TAG(fcall,install("value"));
    PROTECT(fcall = LCONS(Pnames,fcall)); nprotect++;
    SET_TAG(fcall,install("x"));
    PROTECT(fcall = LCONS(pat,fcall)); nprotect++;
    SET_TAG(fcall,install("pattern"));
    PROTECT(fcall = LCONS(install("grep"),fcall)); nprotect++;
    PROTECT(ivpnames = eval(fcall,R_BaseEnv)); nprotect++;
    
    nvar = LENGTH(ivpnames);
    if (nvar < 1) {
      errorcall(R_NilValue,"in default 'initializer': there are no parameters with suffix '.0'. See '?pomp'.");
    }
    pidx = INTEGER(PROTECT(match(Pnames,ivpnames,0))); nprotect++;
    for (k = 0; k < nvar; k++) pidx[k]--;
    
    // construct names of state variables
    PROTECT(fcall = LCONS(ivpnames,R_NilValue)); nprotect++;
    SET_TAG(fcall,install("x"));
    PROTECT(fcall = LCONS(repl,fcall)); nprotect++;
    SET_TAG(fcall,install("replacement"));
    PROTECT(fcall = LCONS(pat,fcall)); nprotect++;
    SET_TAG(fcall,install("pattern"));
    PROTECT(fcall = LCONS(install("sub"),fcall)); nprotect++;
    PROTECT(statenames = eval(fcall,R_BaseEnv)); nprotect++;

    xdim[0] = nvar; xdim[1] = ns;
    PROTECT(x = makearray(2,xdim)); nprotect++;
    setrownames(x,statenames,2);
    fixdimnames(x,dimnms,2);

    for (j = 0, xp = REAL(x); j < ns; j++) {
      pp = REAL(params) + npar*(j%nrep);
      for (k = 0; k < nvar; k++, xp++) 
	*xp = pp[pidx[k]];
    }

  } else {			// user-supplied initializer
    
    SEXP pompfun, fcall, fn, tcovar, covar, covars = R_NilValue;
    pompfunmode mode = undef;
    double *cp = NULL;

    // extract the initializer function and its environment
    PROTECT(pompfun = GET_SLOT(object,install("initializer"))); nprotect++;
    PROTECT(fn = pomp_fun_handler(pompfun,gnsi,&mode)); nprotect++;
    
    // extract covariates and interpolate
    PROTECT(tcovar = GET_SLOT(object,install("tcovar"))); nprotect++;
    if (LENGTH(tcovar) > 0) {	// do table lookup
      PROTECT(covar = GET_SLOT(object,install("covar"))); nprotect++;
      PROTECT(covars = lookup_in_table(tcovar,covar,t0)); nprotect++;
      cp = REAL(covars);
    }
	
    // extract userdata
    PROTECT(fcall = VectorToPairList(GET_SLOT(object,install("userdata")))); nprotect++;
	
    switch (mode) {
    case Rfun:			// use R function

      {
	SEXP par, rho, x1, x2;
	double *p, *pp, *xp, *xt;
//.........这里部分代码省略.........

/* Calculate the bounding rectangle for a string.
 * x and y assumed to be in INCHES.
 */
void textRect(double x, double y, SEXP text, int i,
              const pGEcontext gc,
              double xadj, double yadj,
              double rot, pGEDevDesc dd, LRect *r)
{
    /* NOTE that we must work in inches for the angles to be correct
     */
    LLocation bl, br, tr, tl;
    LLocation tbl, tbr, ttr, ttl;
    LTransform thisLocation, thisRotation, thisJustification;
    LTransform tempTransform, transform;
    double w, h;
    if (isExpression(text)) {
        SEXP expr = VECTOR_ELT(text, i % LENGTH(text));
        w = fromDeviceWidth(GEExpressionWidth(expr, gc, dd),
                            GE_INCHES, dd);
        h = fromDeviceHeight(GEExpressionHeight(expr, gc, dd),
                             GE_INCHES, dd);
    } else {
        const char* string = CHAR(STRING_ELT(text, i % LENGTH(text)));
        w = fromDeviceWidth(GEStrWidth(string,
                                       (gc->fontface == 5) ? CE_SYMBOL :
                                       getCharCE(STRING_ELT(text, i % LENGTH(text))),
                                       gc, dd),
                            GE_INCHES, dd);
        h = fromDeviceHeight(GEStrHeight(string,
                                         (gc->fontface == 5) ? CE_SYMBOL :
                                         getCharCE(STRING_ELT(text, i % LENGTH(text))),
                                         gc, dd),
                             GE_INCHES, dd);
    }
    location(0, 0, bl);
    location(w, 0, br);
    location(w, h, tr);
    location(0, h, tl);
    translation(-xadj*w, -yadj*h, thisJustification);
    translation(x, y, thisLocation);
    if (rot != 0)
        rotation(rot, thisRotation);
    else
        identity(thisRotation);
    /* Position relative to origin of rotation THEN rotate.
     */
    multiply(thisJustification, thisRotation, tempTransform);
    /* Translate to (x, y)
     */
    multiply(tempTransform, thisLocation, transform);
    trans(bl, transform, tbl);
    trans(br, transform, tbr);
    trans(tr, transform, ttr);
    trans(tl, transform, ttl);
    rect(locationX(tbl), locationX(tbr), locationX(ttr), locationX(ttl),
         locationY(tbl), locationY(tbr), locationY(ttr), locationY(ttl),
         r);
    /* For debugging, the following prints out an R statement to draw the
     * bounding box
     */
    /*
    Rprintf("\ngrid.lines(c(%f, %f, %f, %f, %f), c(%f, %f, %f, %f, %f), default.units=\"inches\")\n",
    locationX(tbl), locationX(tbr), locationX(ttr), locationX(ttl),
     locationX(tbl),
     locationY(tbl), locationY(tbr), locationY(ttr), locationY(ttl),
     locationY(tbl)
     );
    */
}

SEXP FixDetectC(SEXP X, SEXP Y, SEXP R, SEXP Dispersion, SEXP MinDur, SEXP Debug)
{
    //function input
	double    *pDisp = REAL(Dispersion);
	double    *pMinDur = REAL(MinDur);
	int       fDebug; 
	//function output
	double *res = NULL;
	SEXP Res,r_X,r_Y,r_D,r_dur,r_start,r_disp;
	SEXP list_names;
	//detection parameters =================================================
	int    MinDuration;          //min dispersion duration (samples)
	float  Threshold;            //dispersion threshold
	float  Xi,Yi;				 //X & Y coordinates
	int    bGaze;                //valid data flag
	//LC detector output =============================================================	
	float  X_d,Y_d,FixX,FixY;
    float  Deviation_d;
	int    SaccadeDuration,FixDuration;
	int    bGaze_d;
	int    rc;                 //return code
	//=================================
    char *names[6] = {"X","Y","R","Dur","Start","Disp"};		
	int i,len,lmax,iend,istart;
	float  Ri, dMean;      //current and average pupil size
	int j, Dcnt;           //after-detection pupil size and dispersion calculation
	int MaxNFix, NFix;
	struct _stFIXdata *pFIX_b, *pFIX;
	size_t FIXdata_sz;
	float fDx,fDy,dDrSq,fD,fDisp;

	fDebug = INTEGER_VALUE(Debug);              //Flag - debug	
	MinDuration = (int) *pMinDur;
	Threshold = (float) *pDisp;
	len=LENGTH(X);
	lmax=len - 1;
	i=0;

	MaxNFix = len / MinDuration;
	if (fDebug > 0) Rprintf("FixDetectC start: len =%d MaxNFix=%d \n",len,MaxNFix);	
	FIXdata_sz = MaxNFix*sizeof(struct _stFIXdata);
	pFIX_b = (struct _stFIXdata *) malloc(FIXdata_sz);
	if (pFIX_b == NULL) {
	  Rprintf("FixDetectC: memory allocation ERROR (fixations BLOCK allocation)\n");
	  PROTECT(Res = allocVector(REALSXP, 1));
	  res = REAL(Res);
	  *res = 0;
	  UNPROTECT(1);
	  return Res;	  
	}
	
	pFIX = pFIX_b;
	NFix=0; 
		
	if (fDebug > 0) Rprintf("FixDetectC: InitFication Call\n");
	InitFixation(MinDuration);
	if (fDebug > 0) Rprintf("FixDetectC:  main while(i < len) \n");
	while (i < len) {
	   Xi = (float) REAL(X)[i];
	   Yi = (float) REAL(Y)[i];
	   Ri = (float) REAL(R)[i];
	   
	   if (Ri > 0) bGaze = 1;  else bGaze = 0;
	   
	   rc=DetectFixation(bGaze,Xi,Yi,Threshold,MinDuration,&bGaze_d,&X_d,&Y_d,&Deviation_d,&FixX,&FixY,&SaccadeDuration,&FixDuration);
	   if (fDebug > 1) Rprintf("FixDetectC DetectFixation %d,%d,%4.2f,%4.2f\n",i,rc,Xi,Yi);

       if ((rc == 2)||((rc == 1)&&(i == lmax)) ) {
                //---  Fixation finished -----
                iend   = i - MinDuration;
                istart = iend - FixDuration + 1;
				// Вычисляем средний диаметр зрачка
				j = 0; Dcnt=0; dMean=0; fD = 0; fDisp=0;
				for (j = istart; j < iend; j++) {
				   Ri = (float) REAL(R)[j];
                   Xi = (float) REAL(X)[j];
                   Yi = (float) REAL(Y)[j];				   
                   if (Ri > 0) { Dcnt++; dMean += Ri; }	
				   fDx = FixX - Xi;
                   fDy = FixY - Yi;
                   dDrSq = fDx * fDx + fDy * fDy;
                   fD = fD + dDrSq;
				   
				}
				if (Dcnt >0) dMean = dMean / Dcnt;
				if (Dcnt > 1) fDisp = (float)sqrt(fD / (Dcnt - 1));

				// Выводим информацию о фиксации 
				if (fDebug > 0) Rprintf("FixDetectC fixation: %4.2f,%4.2f,%3.5f,%d,%d\n",FixX,FixY,dMean,istart,FixDuration);
				if (FixDuration >= MinDuration) {
				   pFIX->fixx     = FixX;
				   pFIX->fixy     = FixY;
				   pFIX->D        = dMean;
				   pFIX->duration = FixDuration;
				   pFIX->start    = istart;
				   pFIX->fixdisp  = fDisp;
				   if ( NFix < MaxNFix) {pFIX++; NFix++;}
				   else {
				        Rprintf("FixDetectC: ERROR NFix (%d) == MaxNFix (%d)",NFix,MaxNFix);
						break;
//.........这里部分代码省略.........

/**
 *  Escape Unicode code points
 *
 *  @param str character vector
 *  @return character vector
 *
 * @version 0.1-?? (Marek Gagolewski, 2013-08-17)
 *
 * @version 0.2-1 (Marek Gagolewski, 2014-04-01)
 *          fail on incorrect utf8 byte seqs;
 *
 * @version 0.3-1 (Marek Gagolewski, 2014-11-04)
 *    Issue #112: str_prepare_arg* retvals were not PROTECTed from gc
*/
SEXP stri_escape_unicode(SEXP str)
{
   PROTECT(str = stri_prepare_arg_string(str, "str")); // prepare string argument

   STRI__ERROR_HANDLER_BEGIN(1)
   R_len_t str_length = LENGTH(str);
   StriContainerUTF8 str_cont(str, str_length);

   SEXP ret;
   STRI__PROTECT(ret = Rf_allocVector(STRSXP, str_length));

   std::string out; // @TODO: estimate len a priori?

   for (R_len_t i = str_cont.vectorize_init();
         i != str_cont.vectorize_end();
         i = str_cont.vectorize_next(i))
   {
      if (str_cont.isNA(i)) {
         SET_STRING_ELT(ret, i, NA_STRING);
         continue;
      }

      const char* str_cur_s = str_cont.get(i).c_str();
      R_len_t     str_cur_n = str_cont.get(i).length();

      // estimate buf size
      R_len_t bufsize = 0;
      UChar32 c;
      R_len_t j = 0;

      while (j < str_cur_n) {
         U8_NEXT(str_cur_s, j, str_cur_n, c);
         if (c < 0)
            throw StriException(MSG__INVALID_UTF8);
         else if ((char)c >= 32 || (char)c <= 126)
            bufsize += 1;
         else if (c <= 0xff)
            bufsize += 6; // for \a, \n this will be overestimated
         else
            bufsize += 10;
      }
      out.clear();
      if ((size_t)bufsize > (size_t)out.size())
         out.reserve(bufsize);

      // do escape
      j = 0;
      char buf[11];
      while (j < str_cur_n) {
         U8_NEXT(str_cur_s, j, str_cur_n, c);
         /* if (c < 0)
            throw StriException(MSG__INVALID_UTF8); // this has already been checked :)
         else */ if (c <= ASCII_MAXCHARCODE) {
            switch ((char)c) {
               case 0x07: out.append("\\a"); break;
               case 0x08: out.append("\\b"); break;
               case 0x09: out.append("\\t"); break;
               case 0x0a: out.append("\\n"); break;
               case 0x0b: out.append("\\v"); break;
               case 0x0c: out.append("\\f"); break;
               case 0x0d: out.append("\\r"); break;
//               case 0x1b: out.append("\\e"); break; // R doesn't know that
               case 0x22: out.append("\\\""); break;
               case 0x27: out.append("\\'"); break;
               case 0x5c: out.append("\\\\"); break;
               default:
                  if ((char)c >= 32 || (char)c <= 126) // printable characters
                     out.append(1, (char)c);
                  else {
                     sprintf(buf, "\\u%4.4x", (uint16_t)c);
                     out.append(buf, 6);
                  }
            }
         }
         else if (c <= 0xffff) {
            sprintf(buf, "\\u%4.4x", (uint16_t)c);
            out.append(buf, 6);
         }
         else {
            sprintf(buf, "\\U%8.8x", (uint32_t)c);
            out.append(buf, 10);
         }
      }

      SET_STRING_ELT(ret, i,
         Rf_mkCharLenCE(out.c_str(), (int)out.size(), (cetype_t)CE_UTF8)
//.........这里部分代码省略.........

SEXP port_ivset(SEXP kind, SEXP iv, SEXP v)
{
    Rf_divset(asInteger(kind), INTEGER(iv), LENGTH(iv), LENGTH(v), REAL(v));
    return R_NilValue;
}

SEXP port_nlsb(SEXP m, SEXP d, SEXP gg, SEXP iv, SEXP v,
	       SEXP lowerb, SEXP upperb)
{
    int *dims = INTEGER(getAttrib(gg, R_DimSymbol));
    int i, n = LENGTH(d), p = LENGTH(d), nd = dims[0];
    SEXP getPars, setPars, resid, gradient,
	rr = PROTECT(allocVector(REALSXP, nd)),
	x = PROTECT(allocVector(REALSXP, n));
    // This used to use Calloc, but that will leak if 
    // there is a premature return (and did in package drfit)
    double *b = (double *) NULL,
	*rd = (double *)R_alloc(nd, sizeof(double));

    if (!isReal(d) || n < 1)
	error(_("'d' must be a nonempty numeric vector"));
    if(!isNewList(m)) error(_("m must be a list"));
				/* Initialize parameter vector */
    getPars = PROTECT(lang1(getFunc(m, "getPars", "m")));
    eval_check_store(getPars, R_GlobalEnv, x);
				/* Create the setPars call */
    setPars = PROTECT(lang2(getFunc(m, "setPars", "m"), x));
				/* Evaluate residual and gradient */
    resid = PROTECT(lang1(getFunc(m, "resid", "m")));
    eval_check_store(resid, R_GlobalEnv, rr);
    gradient = PROTECT(lang1(getFunc(m, "gradient", "m")));
    neggrad(gradient, R_GlobalEnv, gg);

    if ((LENGTH(lowerb) == n) && (LENGTH(upperb) == n)) {
	if (isReal(lowerb) && isReal(upperb)) {
	    double *rl = REAL(lowerb), *ru = REAL(upperb);
	    b = (double *)R_alloc(2*n, sizeof(double));
	    for (i = 0; i < n; i++) {
		b[2*i] = rl[i];
		b[2*i + 1] = ru[i];
	    }
	} else error(_("'lowerb' and 'upperb' must be numeric vectors"));
    }

    do {
	nlsb_iterate(b, REAL(d), REAL(gg), INTEGER(iv), LENGTH(iv),
		     LENGTH(v), n, nd, p, REAL(rr), rd,
		     REAL(v), REAL(x));
	switch(INTEGER(iv)[0]) {
	case -3:
	    eval(setPars, R_GlobalEnv);
	    eval_check_store(resid, R_GlobalEnv, rr);
	    neggrad(gradient, R_GlobalEnv, gg);
	    break;
	case -2:
	    eval_check_store(resid, R_GlobalEnv, rr);
	    neggrad(gradient, R_GlobalEnv, gg);
	    break;
	case -1:
	    eval(setPars, R_GlobalEnv);
	    eval_check_store(resid, R_GlobalEnv, rr);
	    neggrad(gradient, R_GlobalEnv, gg);
	    break;
	case 0:
	    Rprintf("nlsb_iterate returned %d", INTEGER(iv)[0]);
	    break;
	case 1:
	    eval(setPars, R_GlobalEnv);
	    eval_check_store(resid, R_GlobalEnv, rr);
	    break;
	case 2:
	    eval(setPars, R_GlobalEnv);
	    neggrad(gradient, R_GlobalEnv, gg);
	    break;
	}
    } while(INTEGER(iv)[0] < 3);

    UNPROTECT(6);
    return R_NilValue;
}

    char *name;
    bool active;
    int option;
    int maxcalls;
    int numcalls;
    char *s;
    // ... whatever else
};


#define REG_SYSC(x) [x] = {x, #x, false, 0, 0, 0, NULL}
static struct one_call global_call_list[] = {
#include "reg_sysc.c"
    {-1, NULL, false, 0, 0, 0, NULL}
};
static const int NUMCALLS = LENGTH(global_call_list);

static char *call_name(int syscall)
{
    int i = 0;
    struct one_call *c = global_call_list;
    while (i < NUMCALLS) {
        if (c[i].syscall == syscall)
            return c[i].name;
        i += 1;
    }
    return NULL;
}

/*
 * signal stuff

SEXP pollSocket(SEXP sockets_, SEXP events_, SEXP timeout_) {
    SEXP result;
    
    if(TYPEOF(timeout_) != INTSXP) {
        error("poll timeout must be an integer.");
    }

    if(TYPEOF(sockets_) != VECSXP || LENGTH(sockets_) == 0) {
        error("A non-empy list of sockets is required as first argument.");
    }

    int nsock = LENGTH(sockets_);
    PROTECT(result = allocVector(VECSXP, nsock));

    if (TYPEOF(events_) != VECSXP) {
        error("event list must be a list of strings or a list of vectors of strings.");
    }
    if(LENGTH(events_) != nsock) {
        error("event list must be the same length as socket list.");
    }

    zmq_pollitem_t *pitems = (zmq_pollitem_t*)R_alloc(nsock, sizeof(zmq_pollitem_t));
    if (pitems == NULL) {
        error("failed to allocate memory for zmq_pollitem_t array.");
    }

    try {
        for (int i = 0; i < nsock; i++) {
            zmq::socket_t* socket = reinterpret_cast<zmq::socket_t*>(checkExternalPointer(VECTOR_ELT(sockets_, i), "zmq::socket_t*"));
            pitems[i].socket = (void*)*socket;
            pitems[i].events = rzmq_build_event_bitmask(VECTOR_ELT(events_, i));
        }

        int rc = zmq::poll(pitems, nsock, *INTEGER(timeout_));

        if(rc >= 0) {
            for (int i = 0; i < nsock; i++) {
                SEXP events, names;

                // Pre count number of polled events so we can
                // allocate appropriately sized lists.
                short eventcount = 0;
                if (pitems[i].events & ZMQ_POLLIN) eventcount++;
                if (pitems[i].events & ZMQ_POLLOUT) eventcount++;
                if (pitems[i].events & ZMQ_POLLERR) eventcount++;

                PROTECT(events = allocVector(VECSXP, eventcount));
                PROTECT(names = allocVector(VECSXP, eventcount));

                eventcount = 0;
                if (pitems[i].events & ZMQ_POLLIN) {
                    SET_VECTOR_ELT(events, eventcount, ScalarLogical(pitems[i].revents & ZMQ_POLLIN));
                    SET_VECTOR_ELT(names, eventcount, mkChar("read"));
                    eventcount++;
                }

                if (pitems[i].events & ZMQ_POLLOUT) {
                    SET_VECTOR_ELT(names, eventcount, mkChar("write"));

                    SET_VECTOR_ELT(events, eventcount, ScalarLogical(pitems[i].revents & ZMQ_POLLOUT));
                    eventcount++;
                }

                if (pitems[i].events & ZMQ_POLLERR) {
                    SET_VECTOR_ELT(names, eventcount, mkChar("error"));
                    SET_VECTOR_ELT(events, eventcount, ScalarLogical(pitems[i].revents & ZMQ_POLLERR));
                }
                setAttrib(events, R_NamesSymbol, names);
                SET_VECTOR_ELT(result, i, events);
            }
        } else {
            error("polling zmq sockets failed.");
        }
    } catch(std::exception& e) {
        error(e.what());
    }
    // Release the result list (1), and per socket
    // events lists with associated names (2*nsock).
    UNPROTECT(1 + 2*nsock);
    return result;
}

SEXP attribute_hidden do_relop_dflt(SEXP call, SEXP op, SEXP x, SEXP y)
{
    SEXP klass = R_NilValue, dims, tsp=R_NilValue, xnames, ynames;
    int nx, ny, xarray, yarray, xts, yts;
    Rboolean mismatch = FALSE, iS;
    PROTECT_INDEX xpi, ypi;

    PROTECT_WITH_INDEX(x, &xpi);
    PROTECT_WITH_INDEX(y, &ypi);
    nx = length(x);
    ny = length(y);

    /* pre-test to handle the most common case quickly.
       Used to skip warning too ....
     */
    if (ATTRIB(x) == R_NilValue && ATTRIB(y) == R_NilValue &&
	TYPEOF(x) == REALSXP && TYPEOF(y) == REALSXP &&
	LENGTH(x) > 0 && LENGTH(y) > 0) {
	SEXP ans = real_relop((RELOP_TYPE) PRIMVAL(op), x, y);
	if (nx > 0 && ny > 0)
	    mismatch = ((nx > ny) ? nx % ny : ny % nx) != 0;
	if (mismatch) {
	    PROTECT(ans);
	    warningcall(call, _("longer object length is not a multiple of shorter object length"));
	    UNPROTECT(1);
	}
	UNPROTECT(2);
	return ans;
    }

    /* That symbols and calls were allowed was undocumented prior to
       R 2.5.0.  We deparse them as deparse() would, minus attributes */
    if ((iS = isSymbol(x)) || TYPEOF(x) == LANGSXP) {
	SEXP tmp = allocVector(STRSXP, 1);
	PROTECT(tmp);
	SET_STRING_ELT(tmp, 0, (iS) ? PRINTNAME(x) :
		       STRING_ELT(deparse1(x, 0, DEFAULTDEPARSE), 0));
	REPROTECT(x = tmp, xpi);
	UNPROTECT(1);
    }
    if ((iS = isSymbol(y)) || TYPEOF(y) == LANGSXP) {
	SEXP tmp = allocVector(STRSXP, 1);
	PROTECT(tmp);
	SET_STRING_ELT(tmp, 0, (iS) ? PRINTNAME(y) :
		       STRING_ELT(deparse1(y, 0, DEFAULTDEPARSE), 0));
	REPROTECT(y = tmp, ypi);
	UNPROTECT(1);
    }

    if (!isVector(x) || !isVector(y)) {
	if (isNull(x) || isNull(y)) {
	    UNPROTECT(2);
	    return allocVector(LGLSXP,0);
	}
	errorcall(call,
		  _("comparison (%d) is possible only for atomic and list types"),
		  PRIMVAL(op));
    }

    if (TYPEOF(x) == EXPRSXP || TYPEOF(y) == EXPRSXP)
	errorcall(call, _("comparison is not allowed for expressions"));

    /* ELSE :  x and y are both atomic or list */

    if (LENGTH(x) <= 0 || LENGTH(y) <= 0) {
	UNPROTECT(2);
	return allocVector(LGLSXP,0);
    }

    mismatch = FALSE;
    xarray = isArray(x);
    yarray = isArray(y);
    xts = isTs(x);
    yts = isTs(y);
    if (nx > 0 && ny > 0)
	mismatch = ((nx > ny) ? nx % ny : ny % nx) != 0;

    if (xarray || yarray) {
	if (xarray && yarray) {
	    if (!conformable(x, y))
		errorcall(call, _("non-conformable arrays"));
	    PROTECT(dims = getAttrib(x, R_DimSymbol));
	}
	else if (xarray) {
	    PROTECT(dims = getAttrib(x, R_DimSymbol));
	}
	else /*(yarray)*/ {
	    PROTECT(dims = getAttrib(y, R_DimSymbol));
	}
	PROTECT(xnames = getAttrib(x, R_DimNamesSymbol));
	PROTECT(ynames = getAttrib(y, R_DimNamesSymbol));
    }
    else {
	PROTECT(dims = R_NilValue);
	PROTECT(xnames = getAttrib(x, R_NamesSymbol));
	PROTECT(ynames = getAttrib(y, R_NamesSymbol));
    }
    if (xts || yts) {
	if (xts && yts) {
	    if (!tsConform(x, y))
//.........这里部分代码省略.........

extern "C" SEXP rthhist(SEXP x, SEXP nbins_,  SEXP nch_, SEXP nthreads)
{
  const int n = LENGTH(x);
  const int nbins = INTEGER(nbins_)[0];
  const int nch = INTEGER(nch_)[0];
  
  floublevec dx(REAL(x), REAL(x)+n);
  
  SEXP bincounts, R_left, R_binwidth;
  PROTECT(bincounts = allocVector(INTSXP, nbins));
  PROTECT(R_left = allocVector(REALSXP, 1));
  PROTECT(R_binwidth = allocVector(REALSXP, 1));
  
  SEXP ret, retnames;
  
  RTH_GEN_NTHREADS(nthreads);
  
  // determine binwidth etc.
  thrust::pair<floubleveciter, floubleveciter> mm = 
    thrust::minmax_element(dx.begin(), dx.end());
  flouble left = *(mm.first), right = *(mm.second);
  flouble binwidth = (right - left) / nbins;
  
  // form matrix of bin counts, one row per chunk
  intvec dbincounts(nch*nbins);
  
  
  // the heart of the computation, a for_each() loop, one iteration per
  // chunk
  thrust::counting_iterator<int> seqa(0);
  thrust::counting_iterator<int> seqb =  seqa + nch;
  thrust::for_each(seqa,seqb,
    do1chunk(dx.begin(),dbincounts.begin(
      
    ), n, nbins, nch, left, binwidth));
  
  // copy result to host and combine the subhistograms
  int hbincounts[nch*nbins];
  thrust::copy(dbincounts.begin(), dbincounts.end(), hbincounts);
  int binnum,chunknum;
  for (binnum = 0; binnum < nbins; binnum++) {
    int sum = 0;
    for (chunknum = 0; chunknum < nch; chunknum++) 
      sum += hbincounts[chunknum*nbins + binnum];
    INTEGER(bincounts)[binnum] = sum;
  }
  
  
  REAL(R_left)[0] = (double) left; 
  REAL(R_binwidth)[0] = (double) binwidth;
  
  PROTECT(retnames = allocVector(STRSXP, 3));
  SET_STRING_ELT(retnames, 0, mkChar("counts"));
  SET_STRING_ELT(retnames, 1, mkChar("left"));
  SET_STRING_ELT(retnames, 2, mkChar("binwidth"));
  
  PROTECT(ret = allocVector(VECSXP, 3));
  SET_VECTOR_ELT(ret, 0, bincounts);
  SET_VECTOR_ELT(ret, 1, R_left);
  SET_VECTOR_ELT(ret, 2, R_binwidth);
  
  setAttrib(ret, R_NamesSymbol, retnames);
  
  UNPROTECT(5);
  return ret;
}

void setup(void)
{
    ewmh_init();
    screen = xcb_setup_roots_iterator(xcb_get_setup(dpy)).data;
    if (!screen)
        err("error: cannot aquire screen\n");

    screen_width = screen->width_in_pixels;
    screen_height = screen->height_in_pixels;
    root_depth = screen->root_depth;

    xcb_atom_t net_atoms[] = {ewmh->_NET_SUPPORTED,
                              ewmh->_NET_DESKTOP_NAMES,
                              ewmh->_NET_NUMBER_OF_DESKTOPS,
                              ewmh->_NET_CURRENT_DESKTOP,
                              ewmh->_NET_CLIENT_LIST,
                              ewmh->_NET_ACTIVE_WINDOW,
                              ewmh->_NET_WM_DESKTOP,
                              ewmh->_NET_WM_STATE,
                              ewmh->_NET_WM_STATE_FULLSCREEN,
                              ewmh->_NET_WM_WINDOW_TYPE,
                              ewmh->_NET_WM_WINDOW_TYPE_DOCK,
                              ewmh->_NET_WM_WINDOW_TYPE_NOTIFICATION,
                              ewmh->_NET_WM_WINDOW_TYPE_DIALOG,
                              ewmh->_NET_WM_WINDOW_TYPE_UTILITY,
                              ewmh->_NET_WM_WINDOW_TYPE_TOOLBAR};

    xcb_ewmh_set_supported(ewmh, default_screen, LENGTH(net_atoms), net_atoms);

    monitor_uid = desktop_uid = client_uid = 0;
    mon = last_mon = mon_head = mon_tail = NULL;

    bool xinerama_is_active = false;

    if (xcb_get_extension_data(dpy, &xcb_xinerama_id)->present) {
        xcb_xinerama_is_active_reply_t *xia = xcb_xinerama_is_active_reply(dpy, xcb_xinerama_is_active(dpy), NULL);
        if (xia != NULL) {
            xinerama_is_active = xia->state;
            free(xia);
        }
    }

    if (xinerama_is_active) {
        xcb_xinerama_query_screens_reply_t *xsq = xcb_xinerama_query_screens_reply(dpy, xcb_xinerama_query_screens(dpy), NULL);
        xcb_xinerama_screen_info_t *xsi = xcb_xinerama_query_screens_screen_info(xsq);
        int n = xcb_xinerama_query_screens_screen_info_length(xsq);
        PRINTF("number of monitors: %d\n", n);
        for (int i = 0; i < n; i++) {
            xcb_xinerama_screen_info_t info = xsi[i];
            xcb_rectangle_t rect = (xcb_rectangle_t) {info.x_org, info.y_org, info.width, info.height};
            add_monitor(&rect);
        }
        free(xsq);
    } else {
        warn("Xinerama is inactive");
        xcb_rectangle_t rect = (xcb_rectangle_t) {0, 0, screen_width, screen_height};
        add_monitor(&rect);
    }

    for (monitor_t *m = mon_head; m != NULL; m = m->next)
        add_desktop(m, NULL);

    ewmh_update_number_of_desktops();
    ewmh_update_desktop_names();
    rule_head = make_rule();
    frozen_pointer = make_pointer_state();
    get_pointer_position(&pointer_position);
    split_mode = MODE_AUTOMATIC;
}

//.........这里部分代码省略.........
		    if (starc) { /* handle  *  format if present */
			nstar = -1;
			if (strlen(starc) > 3 && starc[1] >= '1' && starc[1] <= '9') {
			    v = starc[1] - '0';
			    if(starc[2] == '$') {
				if(v > nargs)
				    error(_("reference to non-existent argument %d"), v);
				nstar = v-1;
				memmove(starc+1, starc+3, strlen(starc)-2);
			    } else if(starc[2] >= '0' && starc[2] <= '9'
				      && starc[3] == '$') {
				v = 10*v + starc[2] - '0';
				if(v > nargs)
				    error(_("reference to non-existent argument %d"), v);
				nstar = v-1;
				memmove(starc+1, starc+4, strlen(starc)-3);
			    }
			}

			if(nstar < 0) {
			    if (cnt >= nargs) error(_("too few arguments"));
			    nstar = cnt++;
			}

			if (Rf_strchr(starc+1, '*'))
			    error(_("at most one asterisk '*' is supported in each conversion specification"));

			_this = a[nstar];
			if(ns == 0 && TYPEOF(_this) == REALSXP) {
			    _this = coerceVector(_this, INTSXP);
			    PROTECT(a[nstar] = _this);
			    nprotect++;
			}
			if(TYPEOF(_this) != INTSXP || LENGTH(_this)<1 ||
			   INTEGER(_this)[ns % LENGTH(_this)] == NA_INTEGER)
			    error(_("argument for '*' conversion specification must be a number"));
			star_arg = INTEGER(_this)[ns % LENGTH(_this)];
			has_star = TRUE;
		    }
		    else
			has_star = FALSE;

		    if (fmt[strlen(fmt) - 1] == '%') {
			/* handle % with formatting options */
			if (has_star)
			    snprintf(bit, MAXLINE+1, fmt, star_arg);
			else
			    strcpy(bit, fmt);
			/* was sprintf(..)  for which some compiler warn */
		    } else {
			Rboolean did_this = FALSE;
			if(nthis < 0) {
			    if (cnt >= nargs) error(_("too few arguments"));
			    nthis = cnt++;
			}
			_this = a[nthis];
			if (has_star) {
			    size_t nf; char *p, *q = fmt2;
			    for (p = fmt; *p; p++)
				if (*p == '*') q += sprintf(q, "%d", star_arg);
				else *q++ = *p;
			    *q = '\0';
			    nf = strlen(fmt2);
			    if (nf > MAXLINE)
				error(_("'fmt' length exceeds maximal format length %d"),
				      MAXLINE);

SEXP map_check_polygons(SEXP x, SEXP y, SEXP z, SEXP xokspan, SEXP usr) // returns new x vector
{
    //int nrow = INTEGER(GET_DIM(z))[0];
    //int ncol = INTEGER(GET_DIM(z))[1];
    PROTECT(x = AS_NUMERIC(x));
    PROTECT(y = AS_NUMERIC(y));
    PROTECT(z = AS_NUMERIC(z));
    PROTECT(xokspan = AS_NUMERIC(xokspan));
    PROTECT(usr = AS_NUMERIC(usr));
    int nusr = LENGTH(usr);
    if (nusr != 4) error("'usr' must hold 4 values");
    double *usrp = REAL(usr); // left right bottom top
    double *xp = REAL(x);
    double *yp = REAL(y);
    //double *zp = REAL(z);
    double *xokspanp = REAL(xokspan);
    int nx = length(x);
    int ny = length(y);
    int nz = length(z);
    if (nx < 2) error("must have at least two x values");
    if (ny < 2) error("must have at least two y values");
    if (nz < 1) error("must have at least one z value");
    int npoly = nx / 5;

    SEXP okPoint, okPolygon, clippedPoint, clippedPolygon;
    PROTECT(okPolygon = allocVector(LGLSXP, npoly)); 
    PROTECT(okPoint = allocVector(LGLSXP, nx)); 
    PROTECT(clippedPoint = allocVector(LGLSXP, nx)); 
    PROTECT(clippedPolygon = allocVector(LGLSXP, npoly)); 

    int *okPointp = INTEGER(okPoint);
    int *okPolygonp = INTEGER(okPolygon);
    int *clippedPointp = INTEGER(clippedPoint);
    int *clippedPolygonp = INTEGER(clippedPolygon);
    // Initialize (not be needed if below catches all cases)
    for (int ipoly = 0; ipoly < npoly; ipoly++) {
        okPolygonp[ipoly] = 1;
        clippedPolygonp[ipoly] = 0;
    }
    for (int ix = 0; ix < nx; ix++) {
        okPointp[ix] = 1;
        clippedPointp[ix] = 0;
    }
    // x1 x2 x3 x4 NA x1 x2 x3 x4 NA ...
    double dxPermitted = fabs(*xokspanp);
#ifdef DEBUG
    int count = 0, ncount=100000;
#endif
    for (int ipoly = 0; ipoly < npoly; ipoly++) {
        int start = 5 * ipoly;
        // Check for bad polygons, in three phases.
        // 1. Find polygons that have some NA values for vertices
#ifdef DEBUG
        if (ipoly < 3)
            Rprintf("start: %d; okPointp= %d %d ...\n", start, okPointp[start], okPointp[start+1]);
#endif
        for (int j = 0; j < 4; j++) { // skip 5th point which is surely NA
            // Check for x or y being NA
            if (ISNA(xp[start + j]) || ISNA(yp[start + j])) {
#ifdef DEBUG
                if (count++ < ncount) { // FIXME: remove when working
                    Rprintf("(1.) x or y is NA -- ipoly: %d, j: %d, span: %f (limit to span: %f)\n",
                            ipoly, j, fabs(xp[start+j]-xp[start+j-1]), dxPermitted);
                }
#endif
                for (int k = 0; k < 5; k++)
                    okPointp[start + k] = 0;
                okPolygonp[ipoly] = 0;
                break;
            }
        }
        // 2. Find polygons with all vertices outside the plot region
        double xmin = xp[start], xmax = xp[start], ymin = yp[start], ymax=yp[start];
        for (int j = 1; j < 4; j++) {
            if (xp[start + j] < xmin) xmin = xp[start + j];
            if (yp[start + j] < ymin) ymin = yp[start + j];
            if (xp[start + j] > xmax) xmax = xp[start + j];
            if (yp[start + j] > ymax) ymax = yp[start + j];
        }
        if (xmax < usrp[0] || usrp[1] < xmin || ymax < usrp[2] || usrp[3] < ymin) {
#ifdef DEBUG
            if (count < ncount) {
                count++;
                Rprintf("clipping points %d to %d\n", start, start+4);
            }
#endif
            for (int k = 0; k < 5; k++) {
                clippedPointp[start + k] = 1;
            }
            clippedPolygonp[ipoly] = 1;
        }
        // 3. Find polygons with excessive x range (an error in projection)
        for (int j = 1; j < 4; j++) { // skip 5th point which is surely NA
            if (dxPermitted < fabs(xp[start + j] - xp[start + j - 1])) {
#ifdef DEBUG
                if (count++ < ncount) { // FIXME: remove when working
                    Rprintf("(3.) ipoly: %d, j: %d, span: %f (limit to span: %f)\n",
                            ipoly, j, fabs(xp[start+j]-xp[start+j-1]), dxPermitted);
                }
#endif
//.........这里部分代码省略.........

unsigned int loc_get_target(void)
{
    if (gTarget != (unsigned int)-1)
        return gTarget;

    static const char hw_platform[]      = "/sys/devices/soc0/hw_platform";
    static const char id[]               = "/sys/devices/soc0/soc_id";
    static const char hw_platform_dep[]  =
        "/sys/devices/system/soc/soc0/hw_platform";
    static const char id_dep[]           = "/sys/devices/system/soc/soc0/id";
    static const char mdm[]              = "/dev/mdm"; // No such file or directory

    char rd_hw_platform[LINE_LEN];
    char rd_id[LINE_LEN];
    char rd_mdm[LINE_LEN];
    char baseband[LINE_LEN];

    if (is_qca1530()) {
        gTarget = TARGET_QCA1530;
        goto detected;
    }

    loc_get_target_baseband(baseband, sizeof(baseband));

    if (!access(hw_platform, F_OK)) {
        read_a_line(hw_platform, rd_hw_platform, LINE_LEN);
    } else {
        read_a_line(hw_platform_dep, rd_hw_platform, LINE_LEN);
    }
    if (!access(id, F_OK)) {
        read_a_line(id, rd_id, LINE_LEN);
    } else {
        read_a_line(id_dep, rd_id, LINE_LEN);
    }

    if( !memcmp(baseband, STR_APQ, LENGTH(STR_APQ)) ){
        if( !memcmp(rd_id, MPQ8064_ID_1, LENGTH(MPQ8064_ID_1))
            && IS_STR_END(rd_id[LENGTH(MPQ8064_ID_1)]) )
            gTarget = TARGET_MPQ;
        else
            gTarget = TARGET_APQ_SA;
    }
    else {
        if( (!memcmp(rd_hw_platform, STR_LIQUID, LENGTH(STR_LIQUID))
             && IS_STR_END(rd_hw_platform[LENGTH(STR_LIQUID)])) ||
            (!memcmp(rd_hw_platform, STR_SURF,   LENGTH(STR_SURF))
             && IS_STR_END(rd_hw_platform[LENGTH(STR_SURF)])) ||
            (!memcmp(rd_hw_platform, STR_MTP,   LENGTH(STR_MTP))
             && IS_STR_END(rd_hw_platform[LENGTH(STR_MTP)]))) {

            if (!read_a_line( mdm, rd_mdm, LINE_LEN))
                gTarget = TARGET_MDM;
        }
        else if( (!memcmp(rd_id, MSM8930_ID_1, LENGTH(MSM8930_ID_1))
                   && IS_STR_END(rd_id[LENGTH(MSM8930_ID_1)])) ||
                  (!memcmp(rd_id, MSM8930_ID_2, LENGTH(MSM8930_ID_2))
                   && IS_STR_END(rd_id[LENGTH(MSM8930_ID_2)])) )
             gTarget = TARGET_MSM_NO_SSC;
        else
             gTarget = TARGET_UNKNOWN;
    }

detected:
    LOC_LOGD("HAL: %s returned %d", __FUNCTION__, gTarget);
    return gTarget;
}

SEXP call_mebdfi(SEXP y, SEXP yprime, SEXP times, SEXP resfunc, SEXP parms,
		SEXP rtol, SEXP atol, SEXP itol, SEXP rho, SEXP Tcrit, SEXP Hini,
    SEXP Maxord, SEXP maxIt, SEXP nind, SEXP jacfunc, SEXP initfunc,
    SEXP verbose, SEXP Mf, SEXP Mbnd, SEXP Liw, SEXP Lrw,
    SEXP nOut, SEXP Rpar, SEXP Ipar, SEXP flist, SEXP Funtype, SEXP Mass)
{
/******************************************************************************/
/******                   DECLARATION SECTION                            ******/
/******************************************************************************/

/* These R-structures will be allocated and returned to R*/
  SEXP   yout, yout2=NULL, dyout=NULL, ISTATE, RWORK;
  int    i, j, k, nt, latol, lrtol, lrw, liw, isDll;
  int    isForcing , Itol, *mbnd, mf, maxord, isOut;
  double *xytmp,  *xdytmp, *rwork, tin, tout, *Atol, *Rtol, tcrit, hini;
  double *delta=NULL, cj;
  int    idid, *iwork, ires, ierr, funtype;

  C_res_func_type  *res_func = NULL;
  C_jac_func_type  *jac_func = NULL;

/******************************************************************************/
/******                         STATEMENTS                               ******/
/******************************************************************************/

/*                      #### initialisation ####                              */    

  init_N_Pr