C++ R_alloc函数代码示例

SEXP coxfit6(SEXP maxiter2,  SEXP time2,   SEXP status2, 
	     SEXP covar2,    SEXP offset2, SEXP weights2,
	     SEXP strata2,   SEXP method2, SEXP eps2, 
	     SEXP toler2,    SEXP ibeta,    SEXP doscale2) {
    int i,j,k, person;
    
    double **covar, **cmat, **imat;  /*ragged arrays */
    double  wtave;
    double *a, *newbeta;
    double *a2, **cmat2;
    double *scale;
    double  denom=0, zbeta, risk;
    double  temp, temp2;
    int     ndead;  /* actually, the sum of their weights */
    double  newlk=0;
    double  dtime, d2;
    double  deadwt;  /*sum of case weights for the deaths*/
    double  efronwt; /* sum of weighted risk scores for the deaths*/
    int     halving;    /*are we doing step halving at the moment? */
    int     nrisk;   /* number of subjects in the current risk set */
 
    /* copies of scalar input arguments */
    int     nused, nvar, maxiter;
    int     method;
    double  eps, toler;
    int doscale;

    /* vector inputs */
    double *time, *weights, *offset;
    int *status, *strata;
    
    /* returned objects */
    SEXP imat2, means2, beta2, u2, loglik2;
    double *beta, *u, *loglik, *means;
    SEXP sctest2, flag2, iter2;
    double *sctest;
    int *flag, *iter;
    SEXP rlist, rlistnames;
    int nprotect;  /* number of protect calls I have issued */

    /* get local copies of some input args */
    nused = LENGTH(offset2);
    nvar  = ncols(covar2);
    method = asInteger(method2);
    maxiter = asInteger(maxiter2);
    eps  = asReal(eps2);     /* convergence criteria */
    toler = asReal(toler2);  /* tolerance for cholesky */
    doscale = asInteger(doscale2);

    time = REAL(time2);
    weights = REAL(weights2);
    offset= REAL(offset2);
    status = INTEGER(status2);
    strata = INTEGER(strata2);
    
    /*
    **  Set up the ragged arrays and scratch space
    **  Normally covar2 does not need to be duplicated, even though
    **  we are going to modify it, due to the way this routine was
    **  was called.  In this case NAMED(covar2) will =0
    */
    nprotect =0;
    if (NAMED(covar2)>0) {
	PROTECT(covar2 = duplicate(covar2)); 
	nprotect++;
	}
    covar= dmatrix(REAL(covar2), nused, nvar);

    PROTECT(imat2 = allocVector(REALSXP, nvar*nvar)); 
    nprotect++;
    imat = dmatrix(REAL(imat2),  nvar, nvar);
    a = (double *) R_alloc(2*nvar*nvar + 4*nvar, sizeof(double));
    newbeta = a + nvar;
    a2 = newbeta + nvar;
    scale = a2 + nvar;
    cmat = dmatrix(scale + nvar,   nvar, nvar);
    cmat2= dmatrix(scale + nvar +nvar*nvar, nvar, nvar);

    /* 
    ** create output variables
    */ 
    PROTECT(beta2 = duplicate(ibeta));
    beta = REAL(beta2);
    PROTECT(means2 = allocVector(REALSXP, nvar));
    means = REAL(means2);
    PROTECT(u2 = allocVector(REALSXP, nvar));
    u = REAL(u2);
    PROTECT(loglik2 = allocVector(REALSXP, 2)); 
    loglik = REAL(loglik2);
    PROTECT(sctest2 = allocVector(REALSXP, 1));
    sctest = REAL(sctest2);
    PROTECT(flag2 = allocVector(INTSXP, 1));
    flag = INTEGER(flag2);
    PROTECT(iter2 = allocVector(INTSXP, 1));
    iter = INTEGER(iter2);
    nprotect += 7;

    /*
    ** Subtract the mean from each covar, as this makes the regression
    **  much more stable.
//.........这里部分代码省略.........

static SEXP readRegistryKey(HKEY hkey, int depth, int view)
{
    int i, k = 0, size0, *indx;
    SEXP ans, nm, ans0, nm0, tmp, sind;
    DWORD res, nsubkeys, maxsubkeylen, nval, maxvalnamlen, size;
    wchar_t *name;
    HKEY sub;
    REGSAM acc = KEY_READ;

    if (depth <= 0) return mkString("<subkey>");

    if(view == 2) acc |= KEY_WOW64_32KEY;
    else if(view == 3) acc |= KEY_WOW64_64KEY;

    res = RegQueryInfoKey(hkey, NULL, NULL, NULL,
			  &nsubkeys, &maxsubkeylen, NULL, &nval,
			  &maxvalnamlen, NULL, NULL, NULL);
    if (res != ERROR_SUCCESS)
	error("RegQueryInfoKey error code %d: '%s'", (int) res,
	      formatError(res));
    size0 = max(maxsubkeylen, maxvalnamlen) + 1;
    name = (wchar_t *) R_alloc(size0, sizeof(wchar_t));
    PROTECT(ans = allocVector(VECSXP, nval + nsubkeys));
    PROTECT(nm = allocVector(STRSXP, nval+ nsubkeys));
    if (nval > 0) {
	PROTECT(ans0 = allocVector(VECSXP, nval));
	PROTECT(nm0 = allocVector(STRSXP, nval));
	for (i = 0; i < nval; i++) {
	    size = size0;
	    res  = RegEnumValueW(hkey, i, (LPWSTR) name, &size,
				 NULL, NULL, NULL, NULL);
	    if (res != ERROR_SUCCESS) break;
	    SET_VECTOR_ELT(ans0, i, readRegistryKey1(hkey, name));
	    SET_STRING_ELT(nm0, i, mkCharUcs(name));
	}
	/* now sort by name */
	PROTECT(sind = allocVector(INTSXP, nval));  indx = INTEGER(sind);
	for (i = 0; i < nval; i++) indx[i] = i;
	orderVector1(indx, nval, nm0, TRUE, FALSE, R_NilValue);
	for (i = 0; i < nval; i++, k++) {
	    SET_VECTOR_ELT(ans, k, VECTOR_ELT(ans0, indx[i]));
	    if (LENGTH(tmp = STRING_ELT(nm0, indx[i])))
	    	SET_STRING_ELT(nm, k, tmp);
	    else
	    	SET_STRING_ELT(nm, k, mkChar("(Default)"));
	}
	UNPROTECT(3);
    }
    if (nsubkeys > 0) {
	PROTECT(ans0 = allocVector(VECSXP, nsubkeys));
	PROTECT(nm0 = allocVector(STRSXP, nsubkeys));
	for (i = 0; i < nsubkeys; i++) {
	    size = size0;
	    res = RegEnumKeyExW(hkey, i, (LPWSTR) name, &size,
				NULL, NULL, NULL, NULL);
	    if (res != ERROR_SUCCESS) break;
	    res = RegOpenKeyExW(hkey, (LPWSTR) name, 0, acc, &sub);
	    if (res != ERROR_SUCCESS) break;
	    SET_VECTOR_ELT(ans0, i, readRegistryKey(sub, depth-1, view));
	    SET_STRING_ELT(nm0, i, mkCharUcs(name));
	    RegCloseKey(sub);
	}
	/* now sort by name */
	PROTECT(sind = allocVector(INTSXP, nsubkeys));  indx = INTEGER(sind);
	for (i = 0; i < nsubkeys; i++) indx[i] = i;
	orderVector1(indx, nsubkeys, nm0, TRUE, FALSE, R_NilValue);
	for (i = 0; i < nsubkeys; i++, k++) {
	    SET_VECTOR_ELT(ans, k, VECTOR_ELT(ans0, indx[i]));
	    SET_STRING_ELT(nm, k, STRING_ELT(nm0, indx[i]));
	}
	UNPROTECT(3);
    }
    setAttrib(ans, R_NamesSymbol, nm);
    UNPROTECT(2);
    return ans;
}

SEXP call_zvode(SEXP y, SEXP times, SEXP derivfunc, SEXP parms, SEXP rtol,
		SEXP atol, SEXP rho, SEXP tcrit, SEXP jacfunc, SEXP initfunc, 
		SEXP iTask, SEXP rWork, SEXP iWork, SEXP jT, SEXP nOut,
    SEXP lZw, SEXP lRw, SEXP lIw, SEXP Rpar, SEXP Ipar, SEXP flist)
    
{
/******************************************************************************/
/******                   DECLARATION SECTION                            ******/
/******************************************************************************/

  int    i, j, k, nt, latol, lrtol, lrw, liw, lzw;
  double tin, tout, *Atol, *Rtol, ss;
  int    neq, itol, itask, istate, iopt, jt, //mflag, 
         is, isDll, isForcing;
  Rcomplex  *xytmp, *dy = NULL, *zwork;
  int    *iwork, it, ntot, nout;   
  double *rwork;  
  C_zderiv_func_type *zderiv_func;
  C_zjac_func_type   *zjac_func = NULL;

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

  lock_solver(); /* prevent nested call of solvers that have global variables */

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

  //init_N_Protect();
  long int old_N_Protect = save_N_Protected();  

  jt = INTEGER(jT)[0];        
  neq = LENGTH(y);
  nt = LENGTH(times);

  nout  = INTEGER(nOut)[0];
  
/* The output:
    zout and ipar are used to pass output variables (number set by nout)
    followed by other input (e.g. forcing functions) provided 
    by R-arguments rpar, ipar
    ipar[0]: number of output variables, ipar[1]: length of rpar, 
    ipar[2]: length of ipar */

/* is function a dll ?*/
  if (inherits(derivfunc, "NativeSymbol")) {
    isDll = 1;
  } else {
    isDll = 0;
  }

/* initialise output for Complex variables ... */
  initOutComplex(isDll, &nout, &ntot, neq, nOut, Rpar, Ipar);

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

  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];

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

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

  /* global variable */
  //timesteps = (double *) R_alloc(2, sizeof(double));
  for (j=0; j<2; j++) timesteps[j] = 0.;

  lzw = INTEGER(lZw)[0];
  zwork = (Rcomplex *) R_alloc(lzw, sizeof(Rcomplex));

  /* initialise global R-variables... */
  
  PROTECT(cY = allocVector(CPLXSXP , neq) )       ;incr_N_Protect();        
  PROTECT(YOUT = allocMatrix(CPLXSXP,ntot+1,nt))  ;incr_N_Protect();
  
  /**************************************************************************/
  /****** Initialization of Parameters and Forcings (DLL functions)    ******/
  /**************************************************************************/
  initParms(initfunc, parms);
  isForcing = initForcings(flist);

/* pointers to functions zderiv_func and zjac_func, passed to the FORTRAN subroutine */

  if (isDll == 1) { /* DLL address passed to FORTRAN */
    zderiv_func = (C_zderiv_func_type *) R_ExternalPtrAddr(derivfunc);
    /* no need to communicate with R - but output variables set here */      
    if (isOut) {
      dy = (Rcomplex *) R_alloc(neq, sizeof(Rcomplex));
//.........这里部分代码省略.........

int  oncentb(int *xrows, int *xcols, double *x, int *ncenters,
	   double *centers, int *cluster, int *clustersize,
	   int *dist,int *iter,int *itermax, double *par, int *verbose)
{
  int k, m, n, chang, a ,seira, minn;
  double e, h, l, aa, i,ermin,serror, mindist;
  double *dista;
  int  *ordd;

  dista = (double *) R_alloc(*ncenters, sizeof(double));
  ordd = (int *) R_alloc(*ncenters, sizeof(int));
  
  ermin=0.0;
  serror=0.0;


  for(k=0; k<*xrows; k++){
    
    for(m=0; m<*ncenters; m++){
         dista[m]=0.0;
    }
    
     for(m=0; m<*ncenters; m++){
      for(n=0;n<*xcols;n++){
	if(*dist == 0){
	  dista[m] += (x[k+(*xrows)*n] - centers[m +(*ncenters)*n])*(x[k+(*xrows)*n] - centers[m + (*ncenters)*n]); 
	}
	else if(*dist ==1){
          dista[m] += fabs(x[k+(*xrows)*n] - centers[m + (*ncenters)*n]);
	}
      }
     }
     
      /*ORDERING ACCORDING TO THE DISTANCE*/   
      for (m=0;m<*ncenters;m++){
	 ordd[m]=m; 
      }
      chang=1;
       while(chang!=0){
	     chang=0;
	for (m=0;m<(*ncenters-1);m++){
	
	  if (dista[m]> dista[m+1]){
	     aa=dista[m];
	     dista[m]=dista[m+1];
	     dista[m+1]=aa;
	    a=ordd[m];
	    ordd[m]=ordd[m+1];
	    ordd[m+1]=a;
	    chang=chang+1;
	  }
         
	}
       }
       
   
      /*NEW CENTERS*/
      for (m=0;m<*ncenters;m++){
	 seira=ordd[m];
	 /*printf("m:%d\n....ord:%d\n",m,seira  );*/
	  i=(double)(((*iter)-1)*(*xrows)+k)/((*itermax)*(*xrows));
          e=par[0]*pow(par[1]/par[0],i);
          l=par[2]*pow(par[3]/par[2],i);
          h=exp(-(double)m/l);
	  /*	  printf("par: %f,  %f, i:%f\n", par[0], par[2], i);
	  printf("m: %i, seira: %i\n", m,seira); */
	  for (n=0; n<*xcols;n++){
	    centers[seira+(*ncenters)*n]+=e*h*(x[k+(*xrows)*n]-centers[seira+(*ncenters)*n]); 
	  }
	  /*	  printf("\n");*/
      }
  }

  for (k=0;k<*xrows;k++){
      mindist=0.0;/*just to avoid compiling warnings*/
      minn=0; /*the same reason*/
      for (m=0;m<*ncenters;m++){
	  dista[m] = 0.0;
	  for(n=0;n<*xcols;n++){
	      if(*dist == 0){
		  dista[m] += (x[k+(*xrows)*n] -
			       centers[m+(*ncenters)*n])*
		      (x[k+(*xrows)*n] - centers[m + (*ncenters)*n]); 
	      }
	      else if(*dist ==1){
		  dista[m] += fabs(x[k+(*xrows)*n] -
				   centers[m + (*ncenters)*n]);
	      }
	  }
	  if (m == 0)
	  {
	      mindist = dista[0]; minn = 0;
	  }
	  else
	  {
	      if (dista[m] < mindist)
	      {
		  mindist = dista[m];
		  minn = m;
	      }
//.........这里部分代码省略.........

SEXP attribute_hidden do_fft(SEXP call, SEXP op, SEXP args, SEXP env)
{
    SEXP z, d;
    int i, inv, maxf, maxmaxf, maxmaxp, maxp, n, ndims, nseg, nspn;
    double *work;
    int *iwork;

    checkArity(op, args);

    z = CAR(args);

    switch (TYPEOF(z)) {
    case INTSXP:
    case LGLSXP:
    case REALSXP:
	z = coerceVector(z, CPLXSXP);
	break;
    case CPLXSXP:
	if (NAMED(z)) z = duplicate(z);
	break;
    default:
	error(_("non-numeric argument"));
    }
    PROTECT(z);

    /* -2 for forward transform, complex values */
    /* +2 for backward transform, complex values */

    inv = asLogical(CADR(args));
    if (inv == NA_INTEGER || inv == 0)
	inv = -2;
    else
	inv = 2;

    if (LENGTH(z) > 1) {
	if (isNull(d = getAttrib(z, R_DimSymbol))) {  /* temporal transform */
	    n = length(z);
	    fft_factor(n, &maxf, &maxp);
	    if (maxf == 0)
		error(_("fft factorization error"));
	    work = (double*)R_alloc(4 * maxf, sizeof(double));
	    iwork = (int*)R_alloc(maxp, sizeof(int));
	    fft_work(&(COMPLEX(z)[0].r), &(COMPLEX(z)[0].i),
		     1, n, 1, inv, work, iwork);
	}
	else {					     /* spatial transform */
	    maxmaxf = 1;
	    maxmaxp = 1;
	    ndims = LENGTH(d);
	    /* do whole loop just for error checking and maxmax[fp] .. */
	    for (i = 0; i < ndims; i++) {
		if (INTEGER(d)[i] > 1) {
		    fft_factor(INTEGER(d)[i], &maxf, &maxp);
		    if (maxf == 0)
			error(_("fft factorization error"));
		    if (maxf > maxmaxf)
			maxmaxf = maxf;
		    if (maxp > maxmaxp)
			maxmaxp = maxp;
		}
	    }
	    work = (double*)R_alloc(4 * maxmaxf, sizeof(double));
	    iwork = (int*)R_alloc(maxmaxp, sizeof(int));
	    nseg = LENGTH(z);
	    n = 1;
	    nspn = 1;
	    for (i = 0; i < ndims; i++) {
		if (INTEGER(d)[i] > 1) {
		    nspn *= n;
		    n = INTEGER(d)[i];
		    nseg /= n;
		    fft_factor(n, &maxf, &maxp);
		    fft_work(&(COMPLEX(z)[0].r), &(COMPLEX(z)[0].i),
			     nseg, n, nspn, inv, work, iwork);
		}
	    }
	}
    }
    UNPROTECT(1);
    return z;
}

/* A function for reading compact format files.
 * Use with the .Call interface function.
 * Written by Mikko Korpela
 */
SEXP rcompact(SEXP filename){
    char field_id, line[LINE_LENGTH], mplier_str[MPLIER_LENGTH], *found1,
        *found2, *found_leftpar, *found_dot, *found_rightpar, *found_tilde,
        *id_start, *old_point, *point, *point2, *endp, *tmp_name,
        *tmp_comment;
    int i, j, n, first_yr, last_yr, id_length, exponent,
	n_repeats, field_width, n_x_w, n_lines, remainder, idx,
	this_last, *i_first, *i_last;
    long int precision;
    size_t idx2;
    Rboolean n_found, divide;
    long long int read_int;
    double read_double, mplier, *r_mplier, *r_data;
    FILE *f;
    SEXP result, series_id, series_first, series_last, series_mplier,
	series_data, project_comments;
    rwlnode first, *this;
    commentnode comment_first, *comment_this;
    double divisor = 1; /* assign a value to avoid compiler nag */
    int n_content = 0;
    int n_comments = 0;
    Rboolean early_eof = FALSE;

    /* Open the file for reading */
    const char *fname = CHAR(STRING_ELT(filename, 0));
    f = fopen(fname, "r");
    if(f == NULL)
	error(_("Could not open file %s for reading"), fname);

    this = &first;      /* current rwlnode */
    comment_this = &comment_first; /* current commentnode */
    n = 0;              /* number of series */
    first_yr = R_INT_MAX; /* the first year in all data */
    last_yr = R_INT_MIN;  /* the last year in all data */

    /* Each round of the loop reads a header line,
     * then the data lines of the corresponding series
     */
    while(fgets_eol(line, &n_content, LINE_LENGTH, f) != NULL){
	/* In the beginning of the file, if no ~ is found, we assume
	 * the line is a comment. This is the same approach as in the
	 * TRiCYCLE program.
	 */
	while(strchr(line, '~') == NULL){
	    if(n_content > 0){ /* Skip empty lines */
		if(n_comments == R_INT_MAX)
		    error(_("Number of comments exceeds integer range"));
		++n_comments;
		tmp_comment = (char *) R_alloc(n_content+1, sizeof(char));
		strncpy(tmp_comment, line, n_content);
		tmp_comment[n_content] = '\0'; /* Null termination */
		comment_this->text = tmp_comment;
		comment_this->next =
		    (commentnode *) R_alloc(1, sizeof(commentnode));
		comment_this = comment_this->next;
	    }
	    if(fgets_eol(line, &n_content, LINE_LENGTH, f) == NULL){
		early_eof = TRUE;
		break;
	    }
	}
	if(early_eof == TRUE)
	    break;

	if(n == R_INT_MAX)
	    error(_("Number of series exceeds integer range"));

	/* A simple check to point out too long header
	 * lines. Generally, if one line is too long, this function
	 * will probably be unable to parse the next line. In that
	 * case, finding the faulty line may be of some value. Of
	 * course, if the input is generated by some program, lines
	 * are expected to be short enough. Data edited by hand may be
	 * a different case.
	 */
	if(n_content > CONTENT_LENGTH){
	    fclose(f);
	    error(_("Series %d: Header line is too long (max length %d)"),
		  n+1, CONTENT_LENGTH);
	}
	n_found = FALSE;

	/* Find the first '=' character (N or I field) */
	found1 = strchr(line, '=');
	/* Not a header line, not a valid file */
	if(found1 == NULL){
	    fclose(f);
	    error(_("Series %d: No '=' found when header line was expected"),
		  n+1);
	}
	if(found1 == line){
	    fclose(f);
	    error(_("Series %d: No room for number before first '='"), n+1);
	}

	/* Convert the part left of the first '=' to an integer */
//.........这里部分代码省略.........

SEXP nlm(SEXP call, SEXP op, SEXP args, SEXP rho)
{
    SEXP value, names, v, R_gradientSymbol, R_hessianSymbol;

    double *x, *typsiz, fscale, gradtl, stepmx,
	steptol, *xpls, *gpls, fpls, *a, *wrk, dlt;

    int code, i, j, k, itnlim, method, iexp, omsg, msg,
	n, ndigit, iagflg, iahflg, want_hessian, itncnt;


/* .Internal(
 *	nlm(function(x) f(x, ...), p, hessian, typsize, fscale,
 *	    msg, ndigit, gradtol, stepmax, steptol, iterlim)
 */
    function_info *state;

    args = CDR(args);
    PrintDefaults();

    state = (function_info *) R_alloc(1, sizeof(function_info));

    /* the function to be minimized */

    v = CAR(args);
    if (!isFunction(v))
	error(_("attempt to minimize non-function"));
    PROTECT(state->R_fcall = lang2(v, R_NilValue));
    args = CDR(args);

    /* `p' : inital parameter value */

    n = 0;
    x = fixparam(CAR(args), &n);
    args = CDR(args);

    /* `hessian' : H. required? */

    want_hessian = asLogical(CAR(args));
    if (want_hessian == NA_LOGICAL) want_hessian = 0;
    args = CDR(args);

    /* `typsize' : typical size of parameter elements */

    typsiz = fixparam(CAR(args), &n);
    args = CDR(args);

    /* `fscale' : expected function size */

    fscale = asReal(CAR(args));
    if (ISNA(fscale)) error(_("invalid NA value in parameter"));
    args = CDR(args);

    /* `msg' (bit pattern) */
    omsg = msg = asInteger(CAR(args));
    if (msg == NA_INTEGER) error(_("invalid NA value in parameter"));
    args = CDR(args);

    ndigit = asInteger(CAR(args));
    if (ndigit == NA_INTEGER) error(_("invalid NA value in parameter"));
    args = CDR(args);

    gradtl = asReal(CAR(args));
    if (ISNA(gradtl)) error(_("invalid NA value in parameter"));
    args = CDR(args);

    stepmx = asReal(CAR(args));
    if (ISNA(stepmx)) error(_("invalid NA value in parameter"));
    args = CDR(args);

    steptol = asReal(CAR(args));
    if (ISNA(steptol)) error(_("invalid NA value in parameter"));
    args = CDR(args);

    /* `iterlim' (def. 100) */
    itnlim = asInteger(CAR(args));
    if (itnlim == NA_INTEGER) error(_("invalid NA value in parameter"));

    state->R_env = rho;

    /* force one evaluation to check for the gradient and hessian */
    iagflg = 0;			/* No analytic gradient */
    iahflg = 0;			/* No analytic hessian */
    state->have_gradient = 0;
    state->have_hessian = 0;
    R_gradientSymbol = install("gradient");
    R_hessianSymbol = install("hessian");

    /* This vector is shared with all subsequent calls */
    v = allocVector(REALSXP, n);
    for (i = 0; i < n; i++) REAL(v)[i] = x[i];
    SETCADR(state->R_fcall, v);
    SET_NAMED(v, 2); // in case the functions try to alter it
    value = eval(state->R_fcall, state->R_env);

    v = getAttrib(value, R_gradientSymbol);
    if (v != R_NilValue) {
	if (LENGTH(v) == n && (isReal(v) || isInteger(v))) {
	    iagflg = 1;
	    state->have_gradient = 1;
//.........这里部分代码省略.........

int_node* gen_int_node(int_vec* unit)
{  int_node *newnode = (int_node*) R_alloc(1,sizeof(*newnode));
   newnode -> value = unit;
   newnode -> next = NULL;
   return newnode;
}

int* gen_int_array(int la)
{  if(la > 0) 
     return((int*) R_alloc(la,sizeof(int)));
   else return NULL;
}

static void predictions(char **casev,
                        char **namesv,
                        char **treev,
                        char **rulesv,
                        char **costv,
                        int *predv,  /* XXX predictions are character */
			double *confidencev,
			int *trials,
                        char **outputv)
{
    int val;  /* Used by setjmp/longjmp for implementing rbm_exit */

    // Announce ourselves for testing
    // Rprintf("predictions called\n");

    // Initialize the globals
    initglobals();

    // Handles the strbufv data structure
    rbm_removeall();

    // XXX Should this be controlled via an option?
    // Rprintf("Calling setOf\n");
    setOf();

    STRBUF *sb_cases = strbuf_create_full(*casev, strlen(*casev));
	if (rbm_register(sb_cases, "undefined.cases", 0) < 0) {
		error("undefined.cases already exists");
	}

    STRBUF *sb_names = strbuf_create_full(*namesv, strlen(*namesv));
	if (rbm_register(sb_names, "undefined.names", 0) < 0) {
		error("undefined.names already exists");
	}

    if (strlen(*treev)) {
	STRBUF *sb_treev = strbuf_create_full(*treev, strlen(*treev));
	/* XXX should sb_treev be copied? */
	if (rbm_register(sb_treev, "undefined.tree", 0) < 0) {
	    error("undefined.tree already exists");
	}
    } else if (strlen(*rulesv))  {
	STRBUF *sb_rulesv = strbuf_create_full(*rulesv, strlen(*rulesv));
	/* XXX should sb_rulesv be copied? */
	if (rbm_register(sb_rulesv, "undefined.rules", 0) < 0) {
	    error("undefined.rules already exists");
	}
	setrules(1);
    } else {
	error("either a tree or rules must be provided");
    }

    // Create a strbuf using *costv and register it as "undefined.costs"
    if (strlen(*costv) > 0) {
        // Rprintf("registering cost matrix: %s", *costv);
        STRBUF *sb_costv = strbuf_create_full(*costv, strlen(*costv));
        // XXX should sb_costv be copied?
	    if (rbm_register(sb_costv, "undefined.costs", 0) < 0) {
		    error("undefined.cost already exists");
	    }
    } else {
        // Rprintf("no cost matrix to register\n");
    }

    /*
     * We need to initialize rbm_buf before calling any code that
     * might call exit/rbm_exit.
     */
    if ((val = setjmp(rbm_buf)) == 0) {
        // Real work is done here
        // Rprintf("\n\nCalling rpredictmain\n");
        rpredictmain(trials ,predv ,confidencev);

        // Rprintf("predict finished\n\n");
    } else {
        Rprintf("predict code called exit with value %d\n\n", val - JMP_OFFSET);
    }

    // Close file object "Of", and return its contents via argument outputv
    char *outputString = closeOf();
    char *output = R_alloc(strlen(outputString) + 1, 1);
    strcpy(output, outputString);
    *outputv = output;

    // We reinitialize the globals on exit out of general paranoia
    initglobals();
}

SEXP attribute_hidden do_edit(SEXP call, SEXP op, SEXP args, SEXP rho)
{
    int   i, rc;
    ParseStatus status;
    SEXP  x, fn, envir, ed, src, srcfile, Rfn;
    char *filename, *editcmd;
    const char *cmd;
    const void *vmaxsave;
    FILE *fp;
#ifdef Win32
    SEXP ti;
    char *title;
#endif

	checkArity(op, args);

    vmaxsave = vmaxget();

    x = CAR(args); args = CDR(args);
    if (TYPEOF(x) == CLOSXP) envir = CLOENV(x);
    else envir = R_NilValue;
    PROTECT(envir);

    fn = CAR(args); args = CDR(args);
    if (!isString(fn))
	error(_("invalid argument to edit()"));

    if (LENGTH(STRING_ELT(fn, 0)) > 0) {
	const char *ss = translateChar(STRING_ELT(fn, 0));
	filename = R_alloc(strlen(ss), sizeof(char));
	strcpy(filename, ss);
    }
    else filename = DefaultFileName;

    if (x != R_NilValue) {
	if((fp=R_fopen(R_ExpandFileName(filename), "w")) == NULL)
	    errorcall(call, _("unable to open file"));
	if (LENGTH(STRING_ELT(fn, 0)) == 0) EdFileUsed++;
	if (TYPEOF(x) != CLOSXP || isNull(src = getAttrib(x, R_SourceSymbol)))
	    src = deparse1(x, 0, FORSOURCING); /* deparse for sourcing, not for display */
	for (i = 0; i < LENGTH(src); i++)
	    fprintf(fp, "%s\n", translateChar(STRING_ELT(src, i)));
	fclose(fp);
    }
#ifdef Win32
    ti = CAR(args);
#endif
    args = CDR(args);
    ed = CAR(args);
    if (!isString(ed)) errorcall(call, _("argument 'editor' type not valid"));
    cmd = translateChar(STRING_ELT(ed, 0));
    if (strlen(cmd) == 0) errorcall(call, _("argument 'editor' is not set"));
    editcmd = R_alloc(strlen(cmd) + strlen(filename) + 6, sizeof(char));
#ifdef Win32
    if (!strcmp(cmd,"internal")) {
	if (!isString(ti))
	    error(_("'title' must be a string"));
	if (LENGTH(STRING_ELT(ti, 0)) > 0) {
	    title = R_alloc(strlen(CHAR(STRING_ELT(ti, 0)))+1, sizeof(char));
	    strcpy(title, CHAR(STRING_ELT(ti, 0)));
	} else {
	    title = R_alloc(strlen(filename)+1, sizeof(char));
	    strcpy(title, filename);
	}
	Rgui_Edit(filename, CE_NATIVE, title, 1);
    }
    else {
	/* Quote path if necessary */
	if(cmd[0] != '"' && Rf_strchr(cmd, ' '))
	    sprintf(editcmd, "\"%s\" \"%s\"", cmd, filename);
	else
	    sprintf(editcmd, "%s \"%s\"", cmd, filename);
	rc = runcmd(editcmd, CE_NATIVE, 1, 1, NULL, NULL, NULL);
	if (rc == NOLAUNCH)
	    errorcall(call, _("unable to run editor '%s'"), cmd);
	if (rc != 0)
	    warningcall(call, _("editor ran but returned error status"));
    }
#else
    if (ptr_R_EditFile)
	rc = ptr_R_EditFile(filename);
    else {
	sprintf(editcmd, "%s %s", cmd, filename);
	rc = R_system(editcmd);
    }
    if (rc != 0)
	errorcall(call, _("problem with running editor %s"), cmd);
#endif

    if (asLogical(GetOption1(install("keep.source")))) {
	PROTECT(Rfn = findFun(install("readLines"), R_BaseEnv));
	PROTECT(src = lang2(Rfn, ScalarString(mkChar(R_ExpandFileName(filename)))));
	PROTECT(src = eval(src, R_BaseEnv));
	PROTECT(Rfn = findFun(install("srcfilecopy"), R_BaseEnv));
	PROTECT(srcfile = lang3(Rfn, ScalarString(mkChar("<tmp>")), src));
	srcfile = eval(srcfile, R_BaseEnv);
	UNPROTECT(5);
    } else
    	srcfile = R_NilValue;
    PROTECT(srcfile);
//.........这里部分代码省略.........

static void c50(char **namesv,
                char **datav,
                char **costv,
                int *subset,
                int *rules,
                int *utility,
                int *trials,
                int *winnow,
                double *sample,
                int *seed,
                int *noGlobalPruning,
                double *CF,
                int *minCases,
                int *fuzzyThreshold,
                int *earlyStopping,
                char **treev,
                char **rulesv,
                char **outputv)
{
    int val;  /* Used by setjmp/longjmp for implementing rbm_exit */

    // Announce ourselves for testing
    // Rprintf("c50 called\n");

    // Initialize the globals to the values that the c50
    // program would have at the start of execution
    initglobals();

    // Set globals based on the arguments.  This is analogous
    // to parsing the command line in the c50 program.
    setglobals(*subset, *rules, *utility, *trials, *winnow, *sample,
               *seed, *noGlobalPruning, *CF, *minCases, *fuzzyThreshold,
               *earlyStopping, *costv);

    // Handles the strbufv data structure
    rbm_removeall();

    // Deallocates memory allocated by NewCase.
    // Not necessary since it's also called at the end of this function,
    // but it doesn't hurt, and I'm feeling paranoid.
    FreeCases();

    // XXX Should this be controlled via an option?
    // Rprintf("Calling setOf\n");
    setOf();

    // Create a strbuf using *namesv as the buffer.
    // Note that this is a readonly strbuf since we can't
    // extend *namesv.
    STRBUF *sb_names = strbuf_create_full(*namesv, strlen(*namesv));

    // Register this strbuf using the name "undefined.names"
	if (rbm_register(sb_names, "undefined.names", 0) < 0) {
		error("undefined.names already exists");
	}

    // Create a strbuf using *datav and register it as "undefined.data"
    STRBUF *sb_datav = strbuf_create_full(*datav, strlen(*datav));
    // XXX why is sb_datav copied? was that part of my debugging?
    // XXX or is this the cause of the leak?
	if (rbm_register(strbuf_copy(sb_datav), "undefined.data", 0) < 0) {
		error("undefined data already exists");
	}

    // Create a strbuf using *costv and register it as "undefined.costs"
    if (strlen(*costv) > 0) {
        // Rprintf("registering cost matrix: %s", *costv);
        STRBUF *sb_costv = strbuf_create_full(*costv, strlen(*costv));
        // XXX should sb_costv be copied?
	    if (rbm_register(sb_costv, "undefined.costs", 0) < 0) {
		    error("undefined.cost already exists");
	    }
    } else {
        // Rprintf("no cost matrix to register\n");
    }

    /*
     * We need to initialize rbm_buf before calling any code that
     * might call exit/rbm_exit.
     */
    if ((val = setjmp(rbm_buf)) == 0) {

        // Real work is done here
        // Rprintf("Calling c50main\n");
        c50main();

        // Rprintf("c50main finished\n");

        if (*rules == 0) {
            // Get the contents of the the tree file
            STRBUF *treebuf = rbm_lookup("undefined.tree");
            if (treebuf != NULL) {
                char *treeString = strbuf_getall(treebuf);
                char *treeObj = R_alloc(strlen(treeString) + 1, 1);
                strcpy(treeObj, treeString);

                // I think the previous value of *treev will be garbage collected
                *treev = treeObj;
            } else {
                // XXX Should *treev be assigned something in this case?
//.........这里部分代码省略.........

SEXP Dens_bw(const SEXP data_sxp, 
			 const SEXP cluster_sxp, 
			 const SEXP cluster_center_sxp, 
			 const SEXP stdev_sxp, 
			 const SEXP choosen_metric_sxp)
{
	// temporary variables 
	int i, j, k, pos, protect_num=0;
	// some constants 
	int clust_num, dim_num, obj_num;

	double *cluster_center;
	int *cluster_tab;
	cluster_center = REAL(cluster_center_sxp);
	cluster_tab = INTEGER(cluster_sxp);
	

	SEXP data_dim;
	PROTECT( data_dim = getAttrib(data_sxp, R_DimSymbol) );
	protect_num++;

	obj_num = INTEGER(data_dim)[0]; 

	SEXP dim;
	PROTECT( dim = getAttrib(cluster_center_sxp, R_DimSymbol) );
	protect_num++;

	clust_num = INTEGER(dim)[0];
	dim_num = INTEGER(dim)[1];

	SEXP density_middle_point_sxp;
	PROTECT( density_middle_point_sxp = allocMatrix(VECSXP, clust_num, clust_num) );
	protect_num++;
	
	SEXP* midp_coords_sxp;
	midp_coords_sxp = (SEXP*) R_alloc(clust_num*clust_num, sizeof(SEXP) );	

	// initilaize matix of vectors, each vector represent density_middle point between two cluter centers
	double *tmp_tab;
	for(i=0; i<clust_num; i++)
		for(j=0; j<clust_num; j++)
		{
			pos = j + i*clust_num;
			if(i<j) 
			{
				PROTECT( midp_coords_sxp[pos] = allocVector(REALSXP, dim_num) );
				tmp_tab = REAL(midp_coords_sxp[pos]);

				// compute density_middle point between center of cluster i and j
				for(k=0; k<dim_num; k++) 
					tmp_tab[k] = ( cluster_center[i + k*clust_num] + cluster_center[j + k*clust_num] )/2;

			}
			else PROTECT( midp_coords_sxp[pos] = R_NilValue );
			
			protect_num++;
			SET_VECTOR_ELT( density_middle_point_sxp, pos, midp_coords_sxp[pos] );
		}

	// time to get the pointer to the distance function
	double *mean = NULL;
	pMetricFunct metric;
	metric = getMetricFunct(INTEGER(choosen_metric_sxp)[0]);
	if( INTEGER(choosen_metric_sxp)[0] == CORRELATION )
	{
		SEXP mean_sxp;
		PROTECT( mean_sxp = clv_mean(data_sxp, obj_num, dim_num) );
		protect_num++;
		mean = REAL(mean_sxp);
	}

	// vector and matrix create to store information about density of clusters 
	// and "density" of space between centers of clusters 
	SEXP center_density_sxp, midp_density_sxp;
	int *center_density, *midp_density;
	double stdev;
	PROTECT( center_density_sxp = allocVector(INTSXP, clust_num) );
	protect_num++;
	PROTECT( midp_density_sxp = allocMatrix(INTSXP, clust_num, clust_num) );
	protect_num++;
	center_density = INTEGER(center_density_sxp);
	midp_density = INTEGER(midp_density_sxp);
	stdev = REAL(stdev_sxp)[0];

	for(i=0; i<clust_num; i++)
	{
		center_density[i] = 0;
		pos = i*clust_num;
		for(j=0; j<clust_num; j++) midp_density[ j + pos ] = 0;
	}

	// declare return value
	SEXP result_sxp;
	PROTECT( result_sxp = allocVector(REALSXP, 1) );
	protect_num++;

	// compute density over clusters centers and density_middle points between cluster centers 
	int density_center1, density_center2, density_middle;
	double dist;
	double tmp_sum = 0;
//.........这里部分代码省略.........

SEXP ffApply(SEXP result, SEXP data, SEXP margin, SEXP function,
             SEXP nrows, SEXP ncols, int worldRank,
             SEXP out_filename, int worldSize) {

    SEXP ans;

    int my_start, my_end, N, function_nlines;
    int local_check = 0, global_check = 0;

    char *filename, *file_out;
    int  *filesize;
    double *mapped_data_matrix;

    filesize = (int *) R_alloc(1, sizeof(int));

    if(worldRank == MASTER_PROCESS) {
        /* data argument is actually a path to a binary file where data is stored */
        filename = (char *)CHAR((STRING_ELT(data,0)));
        file_out = (char *)CHAR(STRING_ELT(out_filename,0));

        /* function SEXP object is a vector of strings, each element contains
           a single line of the function definition */
        function_nlines = length(function);

    } else {
        filename = (char *) R_alloc(FILENAME_LENGTH, sizeof(char));
        file_out = (char *) R_alloc(FILENAME_LENGTH, sizeof(char));

        PROTECT(nrows = allocVector(INTSXP, 1));
        PROTECT(ncols = allocVector(INTSXP, 1));
        PROTECT(margin = allocVector(INTSXP, 1));
    }

    MPI_Bcast(filename, FILENAME_LENGTH, MPI_CHAR, 0, MPI_COMM_WORLD);
    MPI_Bcast(file_out, FILENAME_LENGTH, MPI_CHAR, 0, MPI_COMM_WORLD);
    MPI_Bcast(INTEGER(nrows), 1, MPI_INT, 0, MPI_COMM_WORLD);
    MPI_Bcast(INTEGER(ncols), 1, MPI_INT, 0, MPI_COMM_WORLD);
    MPI_Bcast(INTEGER(margin), 1, MPI_INT, 0, MPI_COMM_WORLD);
    MPI_Bcast(&function_nlines, 1, MPI_INT, 0, MPI_COMM_WORLD);

    if(worldRank != MASTER_PROCESS) {

        PROTECT(function = allocVector(STRSXP, function_nlines));
    }

    if((mapped_data_matrix = map_file(filename, filesize)) == NULL) {
        local_check = -1;
    }

    /* Check if all processes have successfully mapped the file to memory */
    MPI_Allreduce(&local_check, &global_check, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
    if ( global_check != 0 ) return ScalarInteger(-1);

    /* Matrix dimensions in R are interpreted differently than in C.
       We will refer to R rows and columns ordering, so rows are not alligned
       in memory */

    if (INTEGER(margin)[0] == 1) {
        N = INTEGER(nrows)[0];
    } else if (INTEGER(margin)[0] == 2) {
        N = INTEGER(ncols)[0];
    } else {
        DEBUG("Do not know how to distribute margin number %d\n",
              INTEGER(margin)[0]);
        return R_NilValue;
    }

    loopDistribute(worldRank, worldSize, N, &my_start, &my_end);

    /* Bcast function name or definition, cover case when definition is split into
       several lines and stored as a SEXP string vector */
    bcastRFunction(function, function_nlines, worldRank);

    /* Response container, Vector of SEXPs, margin determines vector length */
    PROTECT(ans = allocVector(VECSXP, 1));

    do_ffApply(ans, mapped_data_matrix, margin, function, my_start, my_end,
               INTEGER(nrows)[0], INTEGER(ncols)[0], worldRank, file_out);

    if(worldRank != MASTER_PROCESS) {
        UNPROTECT(4);
    } else {
        UNPROTECT(1);
    }

    return result;
}

/*  start of AS 182 */
void
forkal(Starma G, int d, int il, double *delta, double *y, double *amse,
       int *ifault)
{
    int p = G->p, q = G->q, r = G->r, n = G->n, np = G->np;
    double *phi = G->phi, *V = G->V, *w = G->w, *xrow = G->xrow;
    double *a, *P, *store;
    int rd = r + d, rz = rd*(rd + 1)/2;
    double phii, phij, sigma2, a1, aa, dt, phijdt, ams, tmp;
    int i, j, k, l, nu = 0;
    int k1;
    int i45, jj, kk, lk, ll;
    int nt;
    int kk1, lk1;
    int ind, jkl, kkk;
    int ind1, ind2;

/*  Finite sample prediction from ARIMA processes. */

/*  This routine will calculate the finite sample predictions
    and their conditional mean square errors for any ARIMA process. */

/*     invoking this routine will calculate the finite sample predictions */
/*     and their conditional mean square errors for any arima process. */

    store = (double *) R_alloc(rd, sizeof(double));
    Free(G->a); G->a = a = Calloc(rd, double);
    Free(G->P); G->P = P = Calloc(rz, double);

/*     check for input faults. */
    *ifault = 0;
    if (p < 0) *ifault = 1;
    if (q < 0) *ifault += 2;
    if (p * p + q * q == 0) *ifault = 4;
    if (r != max(p, q + 1)) *ifault = 5;
    if (np != r * (r + 1) / 2) *ifault = 6;
    if (d < 0) *ifault = 8;
    if (il < 1) *ifault = 11;
    if (*ifault != 0) return;

/*     Find initial likelihood conditions. */

    if (r == 1) {
	a[0] = 0.0;
	V[0] = 1.0;
	P[0] = 1.0 / (1.0 - phi[0] * phi[0]);
    } else starma(G, ifault);

/*     Calculate data transformations */

    nt = n - d;
    if (d > 0) {
	for (j = 0; j < d; j++) {
	    store[j] = w[n - j - 2];
	    if(ISNAN(store[j]))
		error(_("missing value in last %d observations"), d);
	}
	for (i = 0; i < nt; i++) {
	    aa = 0.0;
	    for (k = 0; k < d; ++k) aa -= delta[k] * w[d + i - k - 1];
	    w[i] = w[i + d] + aa;
	}
    }

/*     Evaluate likelihood to obtain final Kalman filter conditions */

    {
	double sumlog = 0.0, ssq = 0.0;
	int nit = 0;
	G->n = nt;
	karma(G, &sumlog, &ssq, 1, &nit);
    }


/*     Calculate m.l.e. of sigma squared */

    sigma2 = 0.0;
    for (j = 0; j < nt; j++) {
	/* MacOS X/gcc 3.5 does/didn't have isnan defined properly */
	tmp = G->resid[j];
	if(!ISNAN(tmp)) { nu++; sigma2 += tmp * tmp; }
    }

    sigma2 /= nu;

/*     reset the initial a and P when differencing occurs */

    if (d > 0) {
	for (i = 0; i < np; i++) xrow[i] = P[i];
	for (i = 0; i < rz; i++) P[i] = 0.0;
	ind = 0;
	for (j = 0; j < r; j++) {
	    k = j * (rd + 1) - j * (j + 1) / 2;
	    for (i = j; i < r; i++) P[k++] = xrow[ind++];
	}
	for (j = 0; j < d; j++) a[r + j] = store[j];
    }

    i45 = 2*rd + 1;
//.........这里部分代码省略.........