C++ R_CheckUserInterrupt函数代码示例

SEXP as_bitstring_integer64(SEXP x_, SEXP ret_){
  int i, n = LENGTH(ret_);
  long long * x = (long long *) REAL(x_);
  unsigned long long mask;
  long long v;
  static char buff[NCHARS_BITS_INTEGER64];
  char * str;
  for(i=0; i<n; i++){
	v = x[i];
	str = buff;
	mask = LEFTBIT_INTEGER64;
    while (mask){
        if (v & mask)
              *str = '1';
          else 
              *str = '0';
        str++;
        mask >>= 1;
    }
    *str = 0;
    SET_STRING_ELT(ret_, i, mkChar(buff)); 
    R_CheckUserInterrupt();
  }
  return ret_;
}

void omxGlobal::reportProgress(const char *context, FitContext *fc)
{
	if (omx_absolute_thread_num() != 0) {
		mxLog("omxGlobal::reportProgress called in a thread context (report this bug to developers)");
		return;
	}

	R_CheckUserInterrupt();

	time_t now = time(0);
	if (Global->maxSeconds > 0 && now > Global->startTime + Global->maxSeconds && !Global->timedOut) {
		Global->timedOut = true;
		Rf_warning("Time limit of %d minutes %d seconds exceeded",
			   Global->maxSeconds/60, Global->maxSeconds % 60);
	}
	if (silent || now - lastProgressReport < 1 || fc->getGlobalComputeCount() == previousComputeCount) return;

	lastProgressReport = now;

	std::string str;
	if (previousReportFit == 0.0 || previousReportFit == fc->fit) {
		str = string_snprintf("%s %d %.6g",
				      context, fc->getGlobalComputeCount(), fc->fit);
	} else {
		str = string_snprintf("%s %d %.6g %.4g",
				      context, fc->getGlobalComputeCount(), fc->fit, fc->fit - previousReportFit);
	}

	reportProgressStr(str.c_str());
	previousReportLength = str.size();
	previousReportFit = fc->fit;
	previousComputeCount = fc->getGlobalComputeCount();
}

void MBC_MCMC_loop(ERGMM_MCMC_Model *model, ERGMM_MCMC_Priors *prior,
                   ERGMM_MCMC_MCMCState *cur, ERGMM_MCMC_MCMCSettings *setting, ERGMM_MCMC_ROutput *outlists)
{
    unsigned int pos=0;
    unsigned int iter, total_iters = setting->sample_size*setting->interval;

    //  Rprintf("Started MCMC loop.\n");
    /* Note that indexing here starts with 1.
       It can be thought of as follows:
       At the end of the updates, we have made iter complete MCMC updates. */
    for(iter=1; iter<=total_iters; iter++) {

        R_CheckUserInterrupt(); // So that CTRL-C can interrupt the run.

        ERGMM_MCMC_CV_up(model,prior,cur);
        ERGMM_MCMC_logp_Z(model, cur->state);

        // If we have a new MLE (actually, the highest likelihood encountered to date), store it.
        if( cur->state->lpZ > outlists->lpZ[0] ) MBC_MCMC_store_iteration(0,model,cur->state,setting,outlists);
        if( cur->state->lpZ + cur->state->lpLV > outlists->lpZ[1] + outlists->lpLV[1] ) MBC_MCMC_store_iteration(1,model,cur->state,setting,outlists);

        /* every interval save the results */
        if((iter % setting->interval) == 0) {
            pos = (iter/setting->interval-1)+MBC_OUTLISTS_RESERVE;

            // Store current iteration.
            MBC_MCMC_store_iteration(pos, model, cur->state, setting, outlists);

        }

    } // end main MCMC loop

    return;

}

void mqmscan(int Nind, int Nmark,int Npheno,int **Geno,int **Chromo, double **Dist, double **Pheno, int **Cofactors, int Backwards, int RMLorML,double Alfa,
             int Emiter, double Windowsize,double Steps, double Stepmi,double Stepma,int NRUN,int out_Naug,int **INDlist, double **QTL, int re_estimate,
             RqtlCrossType rqtlcrosstype,int domi,int verbose){
  int cof_cnt=0;
  MQMMarkerMatrix markers = newMQMMarkerMatrix(Nmark+1,Nind);
  cvector cofactor        = newcvector(Nmark);
  vector mapdistance      = newvector(Nmark);

  MQMCrossType crosstype = determine_MQMCross(Nmark,Nind,(const int **)Geno,rqtlcrosstype);

  change_coding(&Nmark, &Nind, Geno, markers, crosstype); // Change all the markers from R/qtl format to MQM internal

  for (int i=0; i< Nmark; i++) {
    mapdistance[i] = POSITIONUNKNOWN;  // Mapdistances
    mapdistance[i] = Dist[0][i];
    cofactor[i]    = MNOCOF;           // Cofactors
    if (Cofactors[0][i] == 1) {
      cofactor[i] = MCOF;              // Set cofactor
      cof_cnt++;
    }
    if (Cofactors[0][i] == 2) {
      cofactor[i] = MSEX;
      cof_cnt++;
    }
    if (cof_cnt+10 > Nind){ fatal("Setting %d cofactors would leave less than 10 degrees of freedom.\n", cof_cnt); }
  }

  char reestimate = 'y';
  if(re_estimate == 0) reestimate = 'n';

  if (crosstype != CF2) {  // Determine what kind of cross we have
    if (verbose==1) Rprintf("INFO: Dominance setting ignored (setting dominance to 0)\n"); // Update dominance accordingly
    domi = 0;
  }

  bool dominance=false;
  if(domi != 0){ dominance=true; }

  //WE HAVE EVERYTHING START WITH MAIN SCANNING FUNCTION
  analyseF2(Nind, &Nmark, &cofactor, (MQMMarkerMatrix)markers, Pheno[(Npheno-1)], Backwards, QTL, &mapdistance, Chromo, NRUN, RMLorML, Windowsize,
            Steps, Stepmi, Stepma, Alfa, Emiter, out_Naug, INDlist, reestimate, crosstype, dominance, verbose);

  if (re_estimate) {
    if (verbose==1) Rprintf("INFO: Sending back the re-estimated map used during the MQM analysis\n");
    for (int i=0; i< Nmark; i++) {
      Dist[0][i] = mapdistance[i];
    }
  }
  if (Backwards) {
    if (verbose==1) Rprintf("INFO: Sending back the model\n");
    for (int i=0; i< Nmark; i++) { Cofactors[0][i] = cofactor[i]; }
  }

  if(verbose) Rprintf("INFO: All done in C returning to R\n");
  #ifndef STANDALONE
    R_CheckUserInterrupt(); /* check for ^C */
    R_FlushConsole();
  #endif
  return;
}  /* end of function mqmscan */

SEXP as_integer64_bitstring(SEXP x_, SEXP ret_){
  Rboolean naflag = FALSE;
  int i, k, l, n = LENGTH(x_);
  long long * ret = (long long *) REAL(ret_);
  unsigned long long mask;
  long long v;
  const char * str;
  for(i=0; i<n; i++){
    str = CHAR(STRING_ELT(x_, i));
    l = strlen(str);
    if (l>BITS_INTEGER64){
      ret[i] = NA_INTEGER64;
      naflag = TRUE;
      break;
    }
    mask = 1;
    v = 0;
    for (k=l-1; k>=0; k--){
      if (str[k] != '0' &&  str[k] != ' '){
        v |= mask;
      }
      mask <<= 1;
    }
    ret[i] = v;
    R_CheckUserInterrupt();
  }
  if (naflag)warning(BITSTRING_OVERFLOW_WARNING);
  return ret_;
}

/**********************************************************************
 * runningratio
 *
 * Take sum(numerator)/sum(denominator) in sliding window
 *
 * We assume that pos and resultpos are sorted (lo to high)
 **********************************************************************/
void runningratio(int n, double *pos, double *numerator, double *denominator,
		  int n_result, double *resultpos, double *result, double window)
{
  int lo, ns;
  int i, j;
  double top, bottom;
  
  window /= 2.0;

  lo=0; 
  for(i=0; i<n_result; i++) {

    R_CheckUserInterrupt(); /* check for ^C */

    top = bottom = 0.0;  ns=0;
    for(j=lo; j<n; j++) {
      if(pos[j] < resultpos[i]-window) lo = j+1;
      else if(pos[j] > resultpos[i]+window) break;
      else {
	top += numerator[j];
	bottom += denominator[j];
	ns++;
      }
    }

    if(ns==0) result[i] = NA_REAL;
    else result[i] = (top / bottom);

  }

}

static int sock_recv(rsconn_t *c, void *buf, int len) {
    char* cb = (char*) buf;
    if (c->intr && c->s != -1) {
	closesocket(c->s);
	c->s = -1;
	IOerr(c, "previous operation was interrupted, connection aborted");
    }
    while (len > 0) {
	int n = recv(c->s, cb, len, 0);
	/* fprintf(stderr, "sock_recv(%d) = %d [errno=%d]\n", len, n, errno); */
	/* bail out only on non-timeout errors */
	if (n == -1 && errno != EAGAIN && errno != EWOULDBLOCK)
	    return -1;
	if (n == 0)
	    break;
	if (n > 0) {
	    cb += n;
	    len -= n;
	}
	if (len) {
	    c->intr = 1;
	    R_CheckUserInterrupt();
	    c->intr = 0;
	}
    }
    return (int) (cb - (char*)buf);
}

// Compute minimal distances in one direction
void distmap_onesided(int right2left) {
  int i,j,k;

  // initialize vj
  for (i=0;i<height;i++) vj[i]=-1;

  for (j=0;j<width;j++) {
    // compute vj, knowing v(j-1)
    for (i=0;i<height;i++) {
      if (vj[i]<j) {
	k=j;
	if (right2left)	while (k<width) if (a[k+i*width]!=0) k++; else break;
	else while (k<width) if (a[width-1-k+i*width]!=0) k++; else break;
	if (k==width) vj[i]=INT_MAX;
	else vj[i]=k;
      }
    }

    if (right2left) find_ndist(0,height-1,0,height-1,j);
    else {
      for (i=0;i<height;i++) if (vj[i]!=INT_MAX) vj[i]=width-1-vj[i];
      find_ndist(0,height-1,0,height-1,width-1-j);
      for (i=0;i<height;i++) if (vj[i]!=INT_MAX) vj[i]=width-1-vj[i];
    }

    // check for user interruption
    R_CheckUserInterrupt();
  }
}

SEXP hankelize_multi(SEXP U, SEXP V) {
  double *rU = REAL(U), *rV = REAL(V), *rF;
  R_len_t L, K, N, i, count;
  SEXP F;
  int *dimu, *dimv;

  /* Calculate length of inputs and output */
  dimu = INTEGER(getAttrib(U, R_DimSymbol));
  dimv = INTEGER(getAttrib(V, R_DimSymbol));
  L = dimu[0]; K = dimv[0];
  if ((count = dimu[1]) != dimv[1])
    error("Both 'U' and 'V' should have equal number of columns");
  N = K + L - 1;

  /* Allocate buffer for output */
  PROTECT(F = allocMatrix(REALSXP, N, count));
  rF = REAL(F);

  /* Perform the actual hankelization */
  for (i = 0; i < count; ++i) {
    R_CheckUserInterrupt();
    /* TODO: nice target for OpenMP stuff */
    hankelize(rF+i*N, rU+i*L, rV+i*K, L, K);
  }

  UNPROTECT(1);
  return F;
}

/**********************************************************************
 * step_bci
 *
 * Calculate transition probabilities (for all intervals) for
 * the Stahl model
 **********************************************************************/
void step_bci(int n_mar, int n_states, double ***tm, double *d,
              int m, double p, int maxit, double tol)
{
    int i, v1, v2;
    double *the_distinct_tm;
    double *fms_bci_result;
    double lambda1, lambda2, rfp;

    allocate_double(2*m+1, &fms_bci_result);
    allocate_double(3*m+2, &the_distinct_tm);

    for(i=0; i<n_mar-1; i++) {
        R_CheckUserInterrupt(); /* check for ^C */

        lambda1 = d[i]*(1-p)*(double)(m+1)*2.0;
        lambda2 = d[i]*p*2.0;
        rfp = 0.5*(1.0 - exp(-lambda2));

        fms_bci(lambda1, fms_bci_result, m, tol, maxit);
        distinct_tm_bci(lambda1, the_distinct_tm, m, fms_bci_result);

        for(v1=0; v1<n_states; v1++) {
            for(v2=0; v2<n_states; v2++) {
                tm[v1][v2][i] = tm_bci(v1, v2, the_distinct_tm, m);
                if(p > 0)
                    tm[v1][v2][i] = (1.0-rfp)*tm[v1][v2][i] +
                        rfp*tm_bci(v1, (v2+m+1) % (2*m+2), the_distinct_tm, m);
                tm[v1][v2][i] = log(tm[v1][v2][i]);
            }
        }
    }
}

void R_info(int *n_ind, int *n_pos, int *n_gen, double *genoprob, 
	    double *info1, double *info2, int *which)
{
  int i, j, k;
  double ***Genoprob, p, s, ss;

  reorg_genoprob(*n_ind, *n_pos, *n_gen, genoprob, &Genoprob);

  for(i=0; i< *n_pos; i++) {
    R_CheckUserInterrupt(); /* check for ^C */

    info1[i] = info2[i] = 0.0;
    for(j=0; j< *n_ind; j++) {
      s=ss=0.0;
      for(k=0; k< *n_gen; k++) {
	p = Genoprob[k][i][j];
	if(*which != 1) if(p > 0) info1[i] += p*log(p);
	if(*which != 0) { s += p*(double)k; ss += p*(double)(k*k); }
      }
      if(*which != 0) info2[i] += (ss - s*s);
    }
    if(*which != 1) info1[i] /= (double)(*n_ind);
    if(*which != 0) info2[i] /= (double)(*n_ind);
  }
}

static void swapcount(int *m, int *nr, int *nc, int *thin)
{
    int k, ij[4], changed, pm[4] = {1, -1, -1, 1} ;
    int sm[4], ev;
    size_t intcheck;

    /* GetRNGstate(); */

    changed = 0;
    intcheck = 0;
    while (changed < *thin) {
	/* Select a random 2x2 matrix*/
	get2x2((*nr) * (*nc) - 1, *nr, ij);
	for (k = 0; k < 4; k ++)
	    sm[k] = m[ij[k]];
	/* The largest value that can be swapped */
	ev = isDiagFill(sm);
 	if (ev != 0) { 
		for (k = 0; k < 4; k++)
			m[ij[k]] += pm[k]*ev;
		changed++;
	}
	if (intcheck % 10000 == 9999)
	    R_CheckUserInterrupt();
	intcheck++;
    }

    /* PutRNGstate(); */
}

static double
csignrank(int k, int n)
{
    int c, u, j;

#ifndef MATHLIB_STANDALONE
    R_CheckUserInterrupt();
#endif

    u = n * (n + 1) / 2;
    c = (u / 2);

    if (k < 0 || k > u)
	return 0;
    if (k > c)
	k = u - k;

    if (n == 1)
        return 1.;
    if (w[0] == 1.)
        return w[k];

    w[0] = w[1] = 1.;
    for(j = 2; j < n+1; ++j) {
        int i, end = imin2(j*(j+1)/2, c);
	for(i = end; i >= j; --i)
	    w[i] += w[i-j];
    }

    return w[k];
}

/**********************************************************************
 * 
 * convertMWril    Convert simulated RIL genotypes using genotypes in founders
 *                 (and the cross types).  [for a single chr]
 *
 * n_ril     Number of RILs to simulate
 * n_mar     Number of markers
 * n_str     Number of founder strains
 *
 * Parents   SNP data for the founder strains [dim n_mar x n_str]
 * 
 * Geno      On entry, the detailed genotype data; on exit, the 
 *           SNP data written bitwise. [dim n_ril x n_mar]
 * 
 * Crosses   The crosses [n_ril x n_str]
 *
 * all_snps  0/1 indicator of whether all parent genotypes are snps
 *
 * error_prob  Genotyping error probability (used only if all_snps==1)
 *
 * Errors      Error indicators
 *
 **********************************************************************/
void convertMWril(int n_ril, int n_mar, int n_str, 
		  int **Parents, int **Geno, int **Crosses,
		  int all_snps, double error_prob, int **Errors)
{
  int i, j, k, temp;

  for(i=0; i<n_ril; i++) {
    R_CheckUserInterrupt(); /* check for ^C */

    for(j=0; j<n_mar; j++) {

      if(Geno[j][i] < 1 || Geno[j][i] > n_str) {
	if(Geno[j][i] > n_str) 
	  warning("Error in RIL genotype (%d): line %d at marker %d\n", Geno[j][i], i+1, j+1);
	Geno[j][i] = 0;
      }
      else {
	temp = Parents[Geno[j][i]-1][j]; /* SNP genotype of RIL i at marker j */

	if(all_snps && unif_rand() < error_prob) { /* make it an error */
	  temp = 1 - temp;
	  Errors[j][i] = 1;
	}

	Geno[j][i] = 0;
	for(k=0; k<n_str; k++) 
	  if(temp == Parents[Crosses[k][i]-1][j]) 
	    Geno[j][i] += (1 << k);
      }

    }
  }
}

void fitqtl_hk_binary(int n_ind, int n_qtl, int *n_gen, double ***Genoprob,
		      double **Cov, int n_cov, 
		      int *model, int n_int, double *pheno, int get_ests,
		      double *lod, int *df, double *ests, double *ests_covar,
		      double *design_mat, double tol, int maxit) 
{

  /* create local variables */
  int i, j, n_qc, itmp; /* loop variants and temp variables */
  double llik, llik0;
  double *dwork, **Ests_covar;
  int *iwork, sizefull;

  /* initialization */
  sizefull = 1;

  /* calculate the dimension of the design matrix for full model */
  n_qc = n_qtl+n_cov; /* total number of QTLs and covariates */
  /* for additive QTLs and covariates*/
  for(i=0; i<n_qc; i++)
    sizefull += n_gen[i];
  /* for interactions, loop thru all interactions */
  for(i=0; i<n_int; i++) { 
    for(j=0, itmp=1; j<n_qc; j++) {
      if(model[i*n_qc+j])
	itmp *= n_gen[j];
    }
    sizefull += itmp; 
  }

  /* reorganize Ests_covar for easy use later */
  /* and make space for estimates and covariance matrix */
  if(get_ests) 
    reorg_errlod(sizefull, sizefull, ests_covar, &Ests_covar);

  /* allocate memory for working arrays, total memory is
     sizefull*n_ind+6*n_ind+4*sizefull for double array, 
     and sizefull for integer array.
     All memory will be allocated one time and split later */
  dwork = (double *)R_alloc(sizefull*n_ind+6*n_ind+4*sizefull,
			    sizeof(double));
  iwork = (int *)R_alloc(sizefull, sizeof(int));


  /* calculate null model log10 likelihood */
  llik0 = nullLODbin(pheno, n_ind);

  R_CheckUserInterrupt(); /* check for ^C */

  /* fit the model */
  llik = galtLODHKbin(pheno, n_ind, n_gen, n_qtl, Genoprob,
		      Cov, n_cov, model, n_int, dwork, iwork, 
		      sizefull, get_ests, ests, Ests_covar,
		      design_mat, tol, maxit);

  *lod = llik - llik0;

  /* degree of freedom equals to the number of columns of x minus 1 (mean) */
  *df = sizefull - 1;
}