自学教程：C++ unif_rand函数代码示例

SEXP vlmc_sim(SEXP vlmc_R, SEXP nsim_){/* Given a fitted VLMC (with |alphabet| = m), simulate a new series y[1..N],   where N = data_len */    int N = asInteger(nsim_), next_ind = 0, nprot = 1;    if (!isInteger(vlmc_R)) {	vlmc_R = PROTECT(coerceVector(vlmc_R, INTSXP)); nprot++;    }    SEXP res = PROTECT(allocVector(INTSXP, N)); // the result    int *y = INTEGER(res),	m = INTEGER(vlmc_R)[0]; // =  |alphabet|    node_t *top = load_tree(INTEGER(vlmc_R), &next_ind, LENGTH(vlmc_R),		    /*level*/ 0, /*Debug*/ 0);    GetRNGstate();    for (int i = 0; i < N; i++) {	/* Find the context, descending the tree, given y[i-1], y[i-2],... : */	node_t *this, *temp;	int j;	for (j = 1, this = top;	     j <= i && (temp = this->child[y[i - j]]) != NULL;	     j++, this = temp) ;	int count = 0;	double r = (double) this->total * unif_rand();	for (j = 0; j < m; j++) {	    count += this->count[j];	    if (r <= count) {		y[i] = j;	break;	    }	}    }    PutRNGstate();    free_node(top);/* do not leak ! */    UNPROTECT(nprot);    return res;}

void permuteOOB(int m, double *x, int *in, int nsample, int mdim) {/* Permute the OOB part of a variable in x.  * Argument: *   m: the variable to be permuted *   x: the data matrix (variables in rows) *   in: vector indicating which case is OOB *   nsample: number of cases in the data *   mdim: number of variables in the data */    double *tp, tmp;    int i, last, k, nOOB = 0;        tp = (double *)  calloc(nsample, sizeof(double));    for (i = 0; i < nsample; ++i) {	/* make a copy of the OOB part of the data into tp (for permuting) */	if (in[i] == 0) {            tp[nOOB] = x[m + i*mdim];            nOOB++;        }    }    /* Permute tp */    last = nOOB;    for (i = 0; i < nOOB; ++i) {	k = (int) (last * unif_rand());	tmp = tp[last - 1];	tp[last - 1] = tp[k];	tp[k] = tmp;	last--;    }    /* Copy the permuted OOB data back into x. */    nOOB = 0;    for (i = 0; i < nsample; ++i) {	if (in[i] == 0) {            x[m + i*mdim] = tp[nOOB];            nOOB++;	}    }    free(tp);}

/* Written by William Constantine     */mutil_errcode mutil_rand_uniform( void *rand_ptr, double *num_out ){  MUTIL_TRACE("Start mutil_rand_uniform()");  /* avoid lint warning */  (void) rand_ptr;  if( !num_out ) {    MUTIL_ERROR( "NULL pointer for output" );    return MUTIL_ERR_NULL_POINTER;  }  /* do not allow function to return endpoints */  do {    *num_out = unif_rand();  } while( *num_out == 0 || *num_out == 1 );  MUTIL_TRACE("Done with mutil_rand_uniform()");  return MUTIL_ERR_OK;}

double rsignrank(double n){    int i, k;    double r;#ifdef IEEE_754    /* NaNs propagated correctly */    if (ISNAN(n)) return(n);#endif    n = floor(n + 0.5);    if (n < 0) ML_ERR_return_NAN;    if (n == 0)	return(0);    r = 0.0;    k = (int) n;    for (i = 0; i < k; ) {	r += (++i) * floor(unif_rand() + 0.5);    }    return(r);}

int *SampleWithoutReplace(int n, int k, int *result, int *buffer){    if (result == NULL)        result = (int *) R_alloc(k, sizeof(int));    if (buffer == NULL)        buffer = (int *) R_alloc(n, sizeof(int));    for (int i = 0; i < n; i++)        buffer[i] = i;    GetRNGstate();    for (int i = 0; i < k; i++)    {        int j = n * unif_rand();        result[i] = buffer[j];        buffer[j] = buffer[--n];    }    PutRNGstate();    return result;}

/* Grid method samping from tomography line*/void rGrid(	   double *Sample,         /* W_i sampled from each tomography line */                 	   double *W1gi,           /* The grid lines of W1[i] */	   double *W2gi,           /* The grid lines of W2[i] */	   int ni_grid,            /* number of grids for observation i*/	   double *mu,             /* mean vector for normal */ 	   double **InvSigma,      /* Inverse covariance matrix for normal */	   int n_dim)              /* dimension of parameters */{  int j;  double dtemp;  double *vtemp=doubleArray(n_dim);  double *prob_grid=doubleArray(ni_grid);     /* density by grid */  double *prob_grid_cum=doubleArray(ni_grid); /* cumulative density by grid */      dtemp=0;  for (j=0;j<ni_grid;j++){    vtemp[0]=log(W1gi[j])-log(1-W1gi[j]);    vtemp[1]=log(W2gi[j])-log(1-W2gi[j]);    prob_grid[j]=dMVN(vtemp, mu, InvSigma, n_dim, 1) -      log(W1gi[j])-log(W2gi[j])-log(1-W1gi[j])-log(1-W2gi[j]);    prob_grid[j]=exp(prob_grid[j]);    dtemp+=prob_grid[j];    prob_grid_cum[j]=dtemp;  }  for (j=0;j<ni_grid;j++)    prob_grid_cum[j]/=dtemp; /*standardize prob.grid */  /*2 sample W_i on the ith tomo line */  j=0;  dtemp=unif_rand();  while (dtemp > prob_grid_cum[j]) j++;  Sample[0]=W1gi[j];  Sample[1]=W2gi[j];  free(vtemp);  free(prob_grid);  free(prob_grid_cum);}

void TestSetError(double *countts, int *jts, int *clts, int *jet, int ntest,        int nclass, int nvote, double *errts,        int labelts, int *nclts, double *cutoff) {    int j, n, ntie;    double cmax, crit;        for (n = 0; n < ntest; ++n) countts[jts[n]-1 + n*nclass] += 1.0;        /*  Prediction is the class with the maximum votes */    for (n = 0; n < ntest; ++n) {        cmax=0.0;        ntie = 1;        for (j = 0; j < nclass; ++j) {            crit = (countts[j + n*nclass] / nvote) / cutoff[j];            if (crit > cmax) {                jet[n] = j+1;                cmax = crit;            }            /*  Break ties at random: */            if (crit == cmax) {                ntie++;                if (unif_rand() > 1.0 / ntie) {                    jet[n] = j+1;                    cmax = crit;                }            }        }    }    if (labelts) {        zeroDouble(errts, nclass + 1);        for (n = 0; n < ntest; ++n) {            if (jet[n] != clts[n]) {                errts[0] += 1.0;                errts[clts[n]] += 1.0;            }        }        errts[0] /= ntest;        for (n = 1; n <= nclass; ++n) errts[n] /= nclts[n-1];    }}

/* sample W via MH for 2x2 table */void rMH(	 double *W,              /* previous draws */	 double *XY,             /* X_i and Y_i */	 double W1min,           /* lower bound for W1 */	 double W1max,           /* upper bound for W1 */	 double *mu,            /* mean vector for normal */ 	 double **InvSigma,     /* Inverse covariance matrix for normal */	 int n_dim)              /* dimension of parameters */{  int j;  double dens1, dens2, ratio;  double *Sample = doubleArray(n_dim);  double *vtemp = doubleArray(n_dim);  double *vtemp1 = doubleArray(n_dim);    /* sample W_1 from unif(W1min, W1max) */  Sample[0] = runif(W1min, W1max);  Sample[1] = XY[1]/(1-XY[0])-Sample[0]*XY[0]/(1-XY[0]);  for (j = 0; j < n_dim; j++) {    vtemp[j] = log(Sample[j])-log(1-Sample[j]);    vtemp1[j] = log(W[j])-log(1-W[j]);  }  /* acceptance ratio */  dens1 = dMVN(vtemp, mu, InvSigma, n_dim, 1) -    log(Sample[0])-log(Sample[1])-log(1-Sample[0])-log(1-Sample[1]);  dens2 = dMVN(vtemp1, mu, InvSigma, n_dim, 1) -    log(W[0])-log(W[1])-log(1-W[0])-log(1-W[1]);  ratio = fmin2(1, exp(dens1-dens2));    /* accept */  if (unif_rand() < ratio)     for (j=0; j<n_dim; j++)       W[j]=Sample[j];    free(Sample);  free(vtemp);  free(vtemp1);}

// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - |// copula for data with missing values // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - |void copula_NA( double Z[], double K[], int R[], int not_continuous[], int *n, int *p ){    int number = *n, dim = *p, nxp = number * dim, dimp1 = dim + 1;        #pragma omp parallel    {	        double sigma, sd_j, mu_ij, lb, ub, runif_value, pnorm_lb, pnorm_ub;        int i, j;                #pragma omp for        for( int counter = 0; counter < nxp; counter++ )        {               j = counter / number;            i = counter % number;                        if( not_continuous[ j ] )            {                sigma = 1.0 / K[ j * dimp1 ]; // 1.0 / K[ j * dim + j ];                sd_j  = sqrt( sigma );                                get_mean( Z, K, &mu_ij, &sigma, &i, &j, &number, &dim );                                if( R[ counter ] != 0 )                {                    get_bounds_NA( Z, R, &lb, &ub, &i, &j, &number );                                        pnorm_lb     = Rf_pnorm5( lb, mu_ij, sd_j, TRUE, FALSE );                    pnorm_ub     = Rf_pnorm5( ub, mu_ij, sd_j, TRUE, FALSE );                    //runif_value = runif( pnorm_lb, pnorm_ub );                    runif_value  = pnorm_lb + unif_rand() * ( pnorm_ub - pnorm_lb );                    Z[ counter ] = Rf_qnorm5( runif_value, mu_ij, sd_j, TRUE, FALSE );                }else                    Z[ counter ] = mu_ij + norm_rand() * sd_j;  // rnorm( mu_ij, sd_j );            }        }    }}

void R_max_col(double *matrix, int *nr, int *nc, int *maxes){    int	  r, c, m, ntie, n_r = *nr;    double a, b, tol, large;    Rboolean isna, used_random=FALSE;    for (r = 0; r < n_r; r++) {	/* first check row for any NAs and find the largest entry */	large = 0.0;	isna = FALSE;	for (c = 0; c < *nc; c++) {	    a = matrix[r + c * n_r];	    if (ISNAN(a)) { isna = TRUE; break; }	    large = fmax2(large, fabs(a));	}	if (isna) { maxes[r] = NA_INTEGER; continue; }	tol = RELTOL * large;	m = 0;	ntie = 1;	a = matrix[r];	for (c = 1; c < *nc; c++) {	    b = matrix[r + c * n_r];	    if (b >= a + tol) {		ntie = 1;		a = b;		m = c;	    } else if (b >= a - tol) {		ntie++;    		if (!used_random) { GetRNGstate(); used_random = TRUE; }		if (ntie * unif_rand() < 1.) m = c;	    }	}	maxes[r] = m + 1;    }    if(used_random) PutRNGstate();}

int sampleTransaction(int T, int *weights) {	int randomWeight;	int imax = T - 1;	int imin = 0;	int imid;	/* We select a weight value randomly */	randomWeight = ((int) (unif_rand() * (weights[ T - 1 ]))) + 1;	/* We identify the first transaction whose cumulative	 * weight is larger or equal to randomWeight */	/* We look for the element using binary search */	while (imax >= imin) {		imid = imin + ((imax - imin) / 2);		if (weights[ imid ] < randomWeight)			imin = imid + 1;		else if (weights[ imid ] > randomWeight)			imax = imid - 1;		else {			while (imid > 0 && weights[ imid - 1 ] == randomWeight)				imid--;			return imid;		}	}	if (weights[ imid ] > randomWeight)		return imid;	return imid + 1;}

double r_random_walk(double alpha, double tau, double beta, double delta){  double dt=0.0001;  double t,sigma=1;  double p = .5 * (1+((delta*sqrt(dt))/sigma));  double a;  //double q = .5 * (1-((mu*sqrt(dt))/sigma));  int i = 0;  double y = beta*alpha;  while(y < alpha && y > 0)  {    GetRNGstate();    a = unif_rand();    PutRNGstate();    if(a <= p) y = y + sigma*sqrt(dt);    else y = y - sigma*sqrt(dt);    i++;  }  if(y >= alpha) t = (i*dt+tau);  else t = -(i*dt+tau);  return t;}

/** * Generate the direction for "running Shake-and-Bake" according to 1.3.3 of * Boender et al. (1991) */void hitandrun_rsabDir(double *d, Matrix *constr, int index) {	const int inc1 = 1; // for BLAS	int n = constr->nCol - 1;	double c[n]; // the constraint vector	for (int i = 0; i < n; ++i) {		c[i] = *get(constr, index, i);	}  if (n == 1) {    d[0] = -c[0];    return;  }	double r = root(unif_rand(), n - 1);	hitandrun_randDir(d, n); // /~{u} in the paper	double cd = F77_CALL(ddot)(&n, c, &inc1, d, &inc1);	double fd = r / sqrt(1 - cd * cd);	double fc = -(r * cd / sqrt(1 - cd * cd) + sqrt(1 - r * r));	F77_CALL(dscal)(&n, &fd, d, &inc1); // d := fd * d	F77_CALL(daxpy)(&n, &fc, c, &inc1, d, &inc1); // d := fc * c + d}

void rrand_(double *r) { *r = unif_rand(); }

//.........这里部分代码省略.........        /* Print header for running output. */        Rprintf("ntree      OOB");        for (n = 1; n <= nclass; ++n) Rprintf("%7i", n);        if (labelts) {            Rprintf("|    Test");            for (n = 1; n <= nclass; ++n) Rprintf("%7i", n);        }        Rprintf("/n");    }    idxByNnode = 0;    idxByNsample = 0;        //Rprintf("addclass %d, ntree %d, cl[300]=%d", addClass,Ntree,cl[299]);    for(jb = 0; jb < Ntree; jb++) {		//Rprintf("addclass %d, ntree %d, cl[300]=%d", addClass,Ntree,cl[299]);        //printf("jb=%d,/n",jb);        /* Do we need to simulate data for the second class? */        if (addClass) createClass(x, nsample0, nsample, mdim);        do {            zeroInt(nodestatus + idxByNnode, *nrnodes);            zeroInt(treemap + 2*idxByNnode, 2 * *nrnodes);            zeroDouble(xbestsplit + idxByNnode, *nrnodes);            zeroInt(nodeclass + idxByNnode, *nrnodes);            zeroInt(varUsed, mdim);            /* TODO: Put all sampling code into a function. */            /* drawSample(sampsize, nsample, ); */            if (stratify) {  /* stratified sampling */                zeroInt(jin, nsample);                zeroDouble(tclasspop, nclass);                zeroDouble(win, nsample);                if (replace) {  /* with replacement */                    for (n = 0; n < nstrata; ++n) {                        for (j = 0; j < sampsize[n]; ++j) {                            ktmp = (int) (unif_rand() * strata_size[n]);                            k = strata_idx[n][ktmp];                            tclasspop[cl[k] - 1] += classwt[cl[k] - 1];                            win[k] += classwt[cl[k] - 1];                            jin[k] = 1;                        }                    }                } else { /* stratified sampling w/o replacement */                    /* re-initialize the index array */                    zeroInt(strata_size, nstrata);                    for (j = 0; j < nsample; ++j) {                        strata_size[strata[j] - 1] ++;                        strata_idx[strata[j] - 1][strata_size[strata[j] - 1] - 1] = j;                    }                    /* sampling without replacement */                    for (n = 0; n < nstrata; ++n) {                        last = strata_size[n] - 1;                        for (j = 0; j < sampsize[n]; ++j) {                            ktmp = (int) (unif_rand() * (last+1));                            k = strata_idx[n][ktmp];                            swapInt(strata_idx[n][last], strata_idx[n][ktmp]);                            last--;                            tclasspop[cl[k] - 1] += classwt[cl[k]-1];                            win[k] += classwt[cl[k]-1];                            jin[k] = 1;                        }                    }                }            } else {  /* unstratified sampling */                anyEmpty = 0;                ntry = 0;                do {                    zeroInt(jin, nsample);

/* *  REFERENCE * *    Ahrens, J.H. and Dieter, U. *    Extensions of Forsythe's method for random sampling from *    the normal distribution. *    Math. Comput. 27, 927-937. * *    The definitions of the constants a[k], d[k], t[k] and *    h[k] are according to the abovementioned article */double norm_rand(void){    const static double a[32] =    {	0.0000000, 0.03917609, 0.07841241, 0.1177699,	0.1573107, 0.19709910, 0.23720210, 0.2776904,	0.3186394, 0.36012990, 0.40225010, 0.4450965,	0.4887764, 0.53340970, 0.57913220, 0.6260990,	0.6744898, 0.72451440, 0.77642180, 0.8305109,	0.8871466, 0.94678180, 1.00999000, 1.0775160,	1.1503490, 1.22985900, 1.31801100, 1.4177970,	1.5341210, 1.67594000, 1.86273200, 2.1538750    };    const static double d[31] =    {	0.0000000, 0.0000000, 0.0000000, 0.0000000,	0.0000000, 0.2636843, 0.2425085, 0.2255674,	0.2116342, 0.1999243, 0.1899108, 0.1812252,	0.1736014, 0.1668419, 0.1607967, 0.1553497,	0.1504094, 0.1459026, 0.1417700, 0.1379632,	0.1344418, 0.1311722, 0.1281260, 0.1252791,	0.1226109, 0.1201036, 0.1177417, 0.1155119,	0.1134023, 0.1114027, 0.1095039    };    const static double t[31] =    {	7.673828e-4, 0.002306870, 0.003860618, 0.005438454,	0.007050699, 0.008708396, 0.010423570, 0.012209530,	0.014081250, 0.016055790, 0.018152900, 0.020395730,	0.022811770, 0.025434070, 0.028302960, 0.031468220,	0.034992330, 0.038954830, 0.043458780, 0.048640350,	0.054683340, 0.061842220, 0.070479830, 0.081131950,	0.094624440, 0.112300100, 0.136498000, 0.171688600,	0.227624100, 0.330498000, 0.584703100    };    const static double h[31] =    {	0.03920617, 0.03932705, 0.03950999, 0.03975703,	0.04007093, 0.04045533, 0.04091481, 0.04145507,	0.04208311, 0.04280748, 0.04363863, 0.04458932,	0.04567523, 0.04691571, 0.04833487, 0.04996298,	0.05183859, 0.05401138, 0.05654656, 0.05953130,	0.06308489, 0.06737503, 0.07264544, 0.07926471,	0.08781922, 0.09930398, 0.11555990, 0.14043440,	0.18361420, 0.27900160, 0.70104740    };    /*----------- Constants and definitions for  Kinderman - Ramage --- */    /*     *  REFERENCE     *     *    Kinderman A. J. and Ramage J. G. (1976).     *    Computer generation of normal random variables.     *    JASA 71, 893-896.     */#define C1		0.398942280401433#define C2		0.180025191068563#define g(x)		(C1*exp(-x*x/2.0)-C2*(A-x))    const static double A =  2.216035867166471;    double s, u1, w, y, u2, u3, aa, tt, theta, R;    int i;    switch(N01_kind) {    case  AHRENS_DIETER: /* see Reference above */	u1 = unif_rand();	s = 0.0;	if (u1 > 0.5)	    s = 1.0;	u1 = u1 + u1 - s;	u1 *= 32.0;	i = (int) u1;	if (i == 32)	    i = 31;	if (i != 0) {	    u2 = u1 - i;	    aa = a[i - 1];	    while (u2 <= t[i - 1]) {		u1 = unif_rand();		w = u1 * (a[i] - aa);		tt = (w * 0.5 + aa) * w;//.........这里部分代码省略.........

/********************************************************************** *  * random_int *    * Generates a random int integer between "low" and "high", inclusive. * *  Input: *  *    low * *    high * **********************************************************************/int random_int(int low, int high){  return((int)(unif_rand()*(double)(high - low + 1)) + low);}

//     rhyper(NR, NB, n) -- NR 'red', NB 'blue', n drawn, how many are 'red'double rhyper(double nn1in, double nn2in, double kkin){    /* extern double afc(int); */    int nn1, nn2, kk;    int ix; // return value (coerced to double at the very end)    Rboolean setup1, setup2;    /* These should become 'thread_local globals' : */    static int ks = -1, n1s = -1, n2s = -1;    static int m, minjx, maxjx;    static int k, n1, n2; // <- not allowing larger integer par    static double tn;    // II :    static double w;    // III:    static double a, d, s, xl, xr, kl, kr, lamdl, lamdr, p1, p2, p3;    /* check parameter validity */    if(!R_FINITE(nn1in) || !R_FINITE(nn2in) || !R_FINITE(kkin))	ML_ERR_return_NAN;    nn1in = R_forceint(nn1in);    nn2in = R_forceint(nn2in);    kkin  = R_forceint(kkin);    if (nn1in < 0 || nn2in < 0 || kkin < 0 || kkin > nn1in + nn2in)	ML_ERR_return_NAN;    if (nn1in >= INT_MAX || nn2in >= INT_MAX || kkin >= INT_MAX) {	/* large n -- evade integer overflow (and inappropriate algorithms)	   -------- */        // FIXME: Much faster to give rbinom() approx when appropriate; -> see Kuensch(1989)	// Johnson, Kotz,.. p.258 (top) mention the *four* different binomial approximations	if(kkin == 1.) { // Bernoulli	    return rbinom(kkin, nn1in / (nn1in + nn2in));	}	// Slow, but safe: return  F^{-1}(U)  where F(.) = phyper(.) and  U ~ U[0,1]	return qhyper(unif_rand(), nn1in, nn2in, kkin, FALSE, FALSE);    }    nn1 = (int)nn1in;    nn2 = (int)nn2in;    kk  = (int)kkin;    /* if new parameter values, initialize */    if (nn1 != n1s || nn2 != n2s) {	setup1 = TRUE;	setup2 = TRUE;    } else if (kk != ks) {	setup1 = FALSE;	setup2 = TRUE;    } else {	setup1 = FALSE;	setup2 = FALSE;    }    if (setup1) {	n1s = nn1;	n2s = nn2;	tn = nn1 + nn2;	if (nn1 <= nn2) {	    n1 = nn1;	    n2 = nn2;	} else {	    n1 = nn2;	    n2 = nn1;	}    }    if (setup2) {	ks = kk;	if (kk + kk >= tn) {	    k = (int)(tn - kk);	} else {	    k = kk;	}    }    if (setup1 || setup2) {	m = (int) ((k + 1.) * (n1 + 1.) / (tn + 2.));	minjx = imax2(0, k - n2);	maxjx = imin2(n1, k);#ifdef DEBUG_rhyper	REprintf("rhyper(nn1=%d, nn2=%d, kk=%d), setup: floor(mean)= m=%d, jx in (%d..%d)/n",		 nn1, nn2, kk, m, minjx, maxjx);#endif    }    /* generate random variate --- Three basic cases */    if (minjx == maxjx) { /* I: degenerate distribution ---------------- */#ifdef DEBUG_rhyper	REprintf("rhyper(), branch I (degenerate)/n");#endif	ix = maxjx;	goto L_finis; // return appropriate variate    } else if (m - minjx < 10) { // II: (Scaled) algorithm HIN (inverse transformation) ----	const static double scale = 1e25; // scaling factor against (early) underflow	const static double con = 57.5646273248511421;					  // 25*log(10) = log(scale) { <==> exp(con) == scale }	if (setup1 || setup2) {	    double lw; // log(w);  w = exp(lw) * scale = exp(lw + log(scale)) = exp(lw + con)	    if (k < n2) {		lw = afc(n2) + afc(n1 + n2 - k) - afc(n2 - k) - afc(n1 + n2);//.........这里部分代码省略.........

/********************* void MetropolisHastings In this function, theta is a m->n_stats-vector just as in MCMCSample, but now networkstatistics is merely another m->n_stats-vector because this function merely iterates nsteps times through the Markov chain, keeping track of the cumulative change statistics along the way, then returns, leaving the updated change statistics in the networkstatistics vector.  In other words, this function  essentially generates a sample of size one*********************/MCMCStatus MetropolisHastings(MHproposal *MHp,			      double *theta, double *networkstatistics,			      int nsteps, int *staken,			      int fVerbose,			      Network *nwp,			      Model *m) {  unsigned int taken=0, unsuccessful=0;/*  if (fVerbose)    Rprintf("Now proposing %d MH steps... ", nsteps); */  for(unsigned int step=0; step < nsteps; step++) {    MHp->logratio = 0;    (*(MHp->func))(MHp, nwp); /* Call MH function to propose toggles */    if(MHp->toggletail[0]==MH_FAILED){      if(MHp->togglehead[0]==MH_UNRECOVERABLE)	error("Something very bad happened during proposal. Memory has not been deallocated, so restart R soon.");      if(MHp->togglehead[0]==MH_IMPOSSIBLE){	Rprintf("MH Proposal function encountered a configuration from which no toggle(s) can be proposed./n");	return MCMC_MH_FAILED;      }      if(MHp->togglehead[0]==MH_UNSUCCESSFUL){	warning("MH Proposal function failed to find a valid proposal.");	unsuccessful++;	if(unsuccessful>taken*MH_QUIT_UNSUCCESSFUL){	  Rprintf("Too many MH Proposal function failures./n");	  return MCMC_MH_FAILED;	}       	continue;      }    }        if(fVerbose>=5){      Rprintf("Proposal: ");      for(unsigned int i=0; i<MHp->ntoggles; i++)	Rprintf(" (%d, %d)", MHp->toggletail[i], MHp->togglehead[i]);      Rprintf("/n");    }    /* Calculate change statistics,       remembering that tail -> head */    ChangeStats(MHp->ntoggles, MHp->toggletail, MHp->togglehead, nwp, m);    if(fVerbose>=5){      Rprintf("Changes: (");      for(unsigned int i=0; i<m->n_stats; i++)	Rprintf(" %f ", m->workspace[i]);      Rprintf(")/n");    }        /* Calculate inner product */    double ip=0;    for (unsigned int i=0; i<m->n_stats; i++){      ip += theta[i] * m->workspace[i];    }    /* The logic is to set cutoff = ip+logratio ,       then let the MH probability equal min{exp(cutoff), 1.0}.       But we'll do it in log space instead.  */    double cutoff = ip + MHp->logratio;    if(fVerbose>=5){      Rprintf("log acceptance probability: %f + %f = %f/n", ip, MHp->logratio, cutoff);    }        /* if we accept the proposed network */    if (cutoff >= 0.0 || log(unif_rand()) < cutoff) {       if(fVerbose>=5){	Rprintf("Accepted./n");      }      /* Make proposed toggles (updating timestamps--i.e., for real this time) */      for(unsigned int i=0; i < MHp->ntoggles; i++){	ToggleEdge(MHp->toggletail[i], MHp->togglehead[i], nwp);		if(MHp->discord)	  for(Network **nwd=MHp->discord; *nwd!=NULL; nwd++){	    ToggleEdge(MHp->toggletail[i],  MHp->togglehead[i], *nwd);	  }      }      /* record network statistics for posterity */      for (unsigned int i = 0; i < m->n_stats; i++){	networkstatistics[i] += m->workspace[i];      }      taken++;    }else{      if(fVerbose>=5){	Rprintf("Rejected./n");      }//.........这里部分代码省略.........

//.........这里部分代码省略.........				if (branch->one != NULL) branch = branch->one;				else newmodel = 1;			} else {				if (branch->zero != NULL)  branch = branch->zero;				else newmodel = 1;			} 			pmodel  += bit;		}		if (pmodel  == n_sure || pmodel == n + n_sure) {			MH = 1.0/(1.0 - problocal);		}		if (newmodel == 1) {		  new_loc = nUnique;		  PROTECT(Rmodel_m = allocVector(INTSXP,pmodel));		  GetModel_m(Rmodel_m, model, p);		  		  glm_fit = PROTECT(glm_FitModel(X, Y, Rmodel_m, Roffset, Rweights,						 glmfamily, Rcontrol, Ra, Rb, Rs, Rlaplace));			  prior_m = compute_prior_probs(model,pmodel,p, modelprior);					  logmargy = REAL(getListElement(getListElement(glm_fit, "lpy"),"lpY"))[0];		  shrinkage_m = REAL(getListElement(getListElement(glm_fit, "lpy"),							  "shrinkage"))[0];		  postnew = logmargy + log(prior_m);		} else {		  new_loc = branch->where;		  postnew =  REAL(logmarg)[new_loc] + log(REAL(priorprobs)[new_loc]);      		} 		MH *= exp(postnew - postold);		//    Rprintf("MH new %lf old %lf/n", postnew, postold);		if (unif_rand() < MH) {			if (newmodel == 1)  {			  new_loc = nUnique;			  insert_model_tree(tree, vars, n, model, nUnique);			  INTEGER(modeldim)[nUnique] = pmodel;				//Rprintf("model %d: %d variables/n", m, pmodel);				  SetModel2(logmargy, shrinkage_m, prior_m, sampleprobs, logmarg, shrinkage, priorprobs, nUnique);			  SetModel1(glm_fit, Rmodel_m, beta, se, modelspace, deviance, R2, Q, Rintercept, nUnique);			  			  UNPROTECT(2);				  ++nUnique; 			}			old_loc = new_loc;			postold = postnew;			pmodel_old = pmodel;			memcpy(modelold, model, sizeof(int)*p);		} else  {			if (newmodel == 1) UNPROTECT(2);		}		INTEGER(counts)[old_loc] += 1;		for (i = 0; i < n; i++) {			// store in opposite order so nth variable is first 			real_model[n-1-i] = (double) modelold[vars[i].index];			REAL(MCMCprobs)[vars[i].index] += (double) modelold[vars[i].index];		}		m++;	}	for (i = 0; i < n; i++) {		REAL(MCMCprobs)[vars[i].index] /= (double) m;	}		// Compute marginal probabilities  

//.........这里部分代码省略.........	      // Update position of previous last transition in group	      piTransitionPositionInGroup[ppiGroup[iOldGroup][iOldPositionInGroup]] = 		iOldPositionInGroup;	    }	    piTransitionInGroup[iTransition] = -1;	  }	  pdTransitionHazard[iTransition] = dNewHazard;	}		// Get Time to transition	dTime += exp_rand() / pnodeRoot->dPartialAcumHazard;	iUsedRandomNumbers++;		while (dTime >= dTarget) {	  ++iStep;	  // Update the state for the fixed incremented time.	  for(iPlace = 0; iPlace < iPlaces; iPlace++)	    ppdMarking[iPlace][iStep] = pdCrntMarking[iPlace];	  if (iStep == iSectionSteps - 1)	    goto EXIT_LOOP;	  dTarget += dDelta;	  // Force check if user interrupted	  iInterruptCnt = 1;	}	if (! --iInterruptCnt) {	  // Allow user interruption	  R_CheckUserInterrupt();	  iInterruptCnt = 10000000;	}	do {	  // Find group containing firing transition	  double dRnd = unif_rand() * pnodeRoot->dPartialAcumHazard;	  iUsedRandomNumbers++;	  pnode = pnodeRoot;	  do {	    if (dRnd < pnode->left->dPartialAcumHazard) {	      pnode = pnode->left;	    } else {	      dRnd -= pnode->left->dPartialAcumHazard;	      pnode = pnode->right;	    }	      	  } while (pnode->left);	  // Next check is because	  // once in a while it is generated a number that goes past	  // the last group or selects a group with zero elements	  // due to accumulated truncation errors.	  // Discard this random number and try again.	} while (piGroupElm[iGroup = pnode->iGroup] == 0);	double dMaxInGroup = dMinHazard * pow(2, iGroup + 1);	// Find transition in group	while (1) {	  if (! --iInterruptCnt) {	    // Allow user interruption	    R_CheckUserInterrupt();	    iInterruptCnt = 10000000;	  }	  iTransitionPtr = (int) (unif_rand() * piGroupElm[iGroup]);	  iUsedRandomNumbers++;	  iTransition = ppiGroup[iGroup][iTransitionPtr];	  iUsedRandomNumbers++;	  if (pdTransitionHazard[iTransition] > unif_rand() * dMaxInGroup) {	    piTotTransitions[iLastTransition = iTransition]++;	    for(iPlacePtr = 0; iPlacePtr < piSNZxRowTot[iTransition]; iPlacePtr++) {

GBMRESULT CGBM::iterate(    double *adF,    double &dTrainError,    double &dValidError,    double &dOOBagImprove,    int &cNodes,    int cNumClasses,    int cClassIdx){    GBMRESULT hr = GBM_OK;    unsigned long i = 0;    unsigned long cBagged = 0;    int cIdxOff = cClassIdx * (cTrain + cValid); //   for(i=0; i < cTrain + cIdxOff; i++){ adF[i] = 0;}    if(!fInitialized)    {        hr = GBM_FAIL;        goto Error;    }    dTrainError = 0.0;    dValidError = 0.0;    dOOBagImprove = 0.0;    vecpTermNodes.assign(2*cDepth+1,NULL);    // randomly assign observations to the Bag    if (cClassIdx == 0)    {        if (!IsPairwise())        {            // regular instance based training            for(i=0; i<cTrain; i++) /* && (cBagged < cTotalInBag); i++) */            {                if(unif_rand()*(cTrain-i) < cTotalInBag-cBagged)                {                    afInBag[i] = true;                    cBagged++;                }                else                {                    afInBag[i] = false;                }/*                if (cBagged >= cTotalInBag){                    break; 		                } */            }	    std::fill(afInBag + i, afInBag + cTrain, false);        }        else        {            // for pairwise training, sampling is per group            // therefore, we will not have exactly cTotalInBag instances            double dLastGroup = -1;            bool chosen = false;            unsigned int cBaggedGroups = 0;            unsigned int cSeenGroups   = 0;            unsigned int cTotalGroupsInBag = (unsigned long)(dBagFraction * cGroups);            if (cTotalGroupsInBag <= 0)            {                cTotalGroupsInBag = 1;            }            for(i=0; i<cTrain; i++)            {                const double dGroup = pData->adMisc[i];                if (dGroup != dLastGroup)                {                    if (cBaggedGroups >= cTotalGroupsInBag)                    {                        break;                    }                    // Group changed, make a new decision                    chosen = (unif_rand()*(cGroups - cSeenGroups) < cTotalGroupsInBag - cBaggedGroups);                    if (chosen)                    {                        cBaggedGroups++;                    }                    dLastGroup = dGroup;                    cSeenGroups++;                }                if (chosen)                {                    afInBag[i] = true;                    cBagged++;                }                else                {                    afInBag[i] = false;                }            }            // the remainder is not in the bag	    std::fill(afInBag + i, afInBag + cTrain, false);        }    }//.........这里部分代码省略.........

/*  Define the R RNG for use from Fortran. */void F77_SUB(rrand)(double *r) { *r = unif_rand(); }

voidrcont2(int *nrow, int *ncol,       /* vectors of row and column totals, and their sum ntotal: */       int *nrowt, int *ncolt, int *ntotal,       double *fact, int *jwork, int *matrix){    int j, l, m, ia, ib, ic, jc, id, ie, ii, nll, nlm, nr_1, nc_1;    double x, y, dummy, sumprb;    Rboolean lsm, lsp;    nr_1 = *nrow - 1;    nc_1 = *ncol - 1;    ib = 0; /* -Wall */    /* Construct random matrix */    for (j = 0; j < nc_1; ++j)	jwork[j] = ncolt[j];    jc = *ntotal;    for (l = 0; l < nr_1; ++l) { /* -----  matrix[ l, * ] ----- */	ia = nrowt[l];	ic = jc;	jc -= ia;/* = n_tot - sum(nr[0:l]) */	for (m = 0; m < nc_1; ++m) {	    id = jwork[m];	    ie = ic;	    ic -= id;	    ib = ie - ia;	    ii = ib - id;	    if (ie == 0) { /* Row [l,] is full, fill rest with zero entries */		for (j = m; j < nc_1; ++j)		    matrix[l + j * *nrow] = 0;		ia = 0;		break;	    }	    /* Generate pseudo-random number */	    dummy = unif_rand();	    do {/* Outer Loop */		/* Compute conditional expected value of MATRIX(L, M) */		nlm = (int)(ia * (id / (double) ie) + 0.5);		x = exp(fact[ia] + fact[ib] + fact[ic] + fact[id]			- fact[ie] - fact[nlm]			- fact[id - nlm] - fact[ia - nlm] - fact[ii + nlm]);		if (x >= dummy)		    break;		if (x == 0.)/* MM: I haven't seen this anymore */		    error(_("rcont2 [%d,%d]: exp underflow to 0; algorithm failure"), l, m);		sumprb = x;		y = x;		nll = nlm;		do {		    /* Increment entry in row L, column M */		    j = (int)((id - nlm) * (double)(ia - nlm));		    lsp = (j == 0);		    if (!lsp) {			++nlm;			x = x * j / ((double) nlm * (ii + nlm));			sumprb += x;			if (sumprb >= dummy)			    goto L160;		    }		    do {			R_CheckUserInterrupt();			/* Decrement entry in row L, column M */			j = (int)(nll * (double)(ii + nll));			lsm = (j == 0);			if (!lsm) {			    --nll;			    y = y * j / ((double) (id - nll) * (ia - nll));			    sumprb += y;			    if (sumprb >= dummy) {				nlm = nll;				goto L160;			    }			    /* else */			    if (!lsp)				break;/* to while (!lsp) */			}		    } while (!lsm);		} while (!lsp);		dummy = sumprb * unif_rand();	    } while (1);L160:	    matrix[l + m * *nrow] = nlm;//.........这里部分代码省略.........

/* Cross-Validation-routine from svm-train */void do_cross_validation(struct svm_problem *prob,			 struct svm_parameter *param,			 int nr_fold,			 double* cresults,			 double* ctotal1,			 double* ctotal2){	int i;	int total_correct = 0;	double total_error = 0;	double sumv = 0, sumy = 0, sumvv = 0, sumyy = 0, sumvy = 0;	/* random shuffle */	GetRNGstate();	for(i=0; i<prob->l; i++)	{	        int j = i+((int) (unif_rand() * (prob->l-i)))%(prob->l-i);		struct svm_node *tx;		double ty;					tx = prob->x[i];		prob->x[i] = prob->x[j];		prob->x[j] = tx;		ty = prob->y[i];		prob->y[i] = prob->y[j];		prob->y[j] = ty;	}	PutRNGstate();	for(i=0; i<nr_fold; i++)	{		int begin = i*prob->l/nr_fold;		int end = (i+1)*prob->l/nr_fold;		int j,k;		struct svm_problem subprob;		subprob.l = prob->l-(end-begin);		subprob.x = Malloc(struct svm_node*,subprob.l);		subprob.y = Malloc(double,subprob.l);					k=0;		for(j = 0; j < begin; j++)		{			subprob.x[k] = prob->x[j];			subprob.y[k] = prob->y[j];			++k;		}		for(j = end; j<prob->l; j++)		{			subprob.x[k] = prob->x[j];			subprob.y[k] = prob->y[j];			++k;		}		if(param->svm_type == EPSILON_SVR ||		   param->svm_type == NU_SVR)		{			struct svm_model *submodel = svm_train(&subprob,param);			double error = 0;			for(j=begin;j<end;j++)			{				double v = svm_predict(submodel,prob->x[j]);				double y = prob->y[j];				error += (v-y)*(v-y);				sumv += v;				sumy += y;				sumvv += v*v;				sumyy += y*y;				sumvy += v*y;			}			svm_free_and_destroy_model(&submodel);			/* printf("Mean squared error = %g/n",			   error/(end-begin)); */			cresults[i] = error/(end-begin);			total_error += error;					}		else		{			struct svm_model *submodel = svm_train(&subprob,param);			int correct = 0;			for(j=begin;j<end;j++)			{				double v = svm_predict(submodel,prob->x[j]);				if(v == prob->y[j])					++correct;			}			svm_free_and_destroy_model(&submodel);			/* printf("Accuracy = %g%% (%d/%d)/n", */			/* 100.0*correct/(end-begin),correct,(end-begin)); */			cresults[i] = 100.0*correct/(end-begin);			total_correct += correct;		}		free(subprob.x);		free(subprob.y);	}		if(param->svm_type == EPSILON_SVR || param->svm_type == NU_SVR)//.........这里部分代码省略.........

size_t rewire_ex(unsigned short *incidence,size_t ncol, size_t nrow,size_t max_iter,size_t verbose,size_t MAXITER,unsigned int seed){	set_rand(seed); 	size_t i,j,kk,n,rand1,rand2,t=0;	size_t e=0;    //copy of the original incidence matrix	size_t *from;	size_t *to;	size_t a,b,c,d;	for(i=0;i<nrow;i++)		for(j=0;j<ncol;j++)			e+=incidence[i*ncol+j];	e/=2;	//initialization of score vector overwriting the original TO CHECK	do from=(size_t*)calloc(e,sizeof(size_t)); while(from==NULL);	do to=(size_t*)calloc(e,sizeof(size_t));   while(to==NULL);	kk=0;	for(i=0;i<nrow;++i)		for(j=0;j<i;++j)			if(incidence[i*ncol+j]==1)				{					from[kk]=i;					to[kk]=j;          kk++;        }	time_t  tin,tfin;  tin = time (NULL); 	//GetRNGstate();	for(n=0;n<max_iter;t++)		{			if(verbose==1)				loadBar( n,  max_iter, 100,  50);		//random rewiring		  rand1=(size_t) (unif_rand()*e);      do rand2=(size_t) (unif_rand()*e); while (rand1==rand2);      a=from[rand1];      c=from[rand2];      b=to[rand1];      d=to[rand2];	if(t>MAXITER)			{				tfin = time (NULL); 				if(verbose==1)    		printf("DONE in %d seconds /n",-(tin-tfin));				//PutRNGstate();				return (-1);			}		//  printf("%d %d %d %d %d %d %d %d %d %d/n ",rand1,rand2,a+1,b+1,c+1,d+1,incidence[a*ncol+d],incidence[c*ncol+b],incidence[a*ncol+c],incidence[d*ncol+b]);     if(a!=c && b!=d&&  a!=d && c!=b &&           (	(incidence[a*ncol+d]==0 && incidence[c*ncol+b]==0  ) ||            (incidence[a*ncol+c]==0 && incidence[d*ncol+b]==0  ) ))        {            if(incidence[a*ncol+d]==0 && incidence[c*ncol+b]==0 && incidence[a*ncol+c]==0 && incidence[d*ncol+b]==0 )            {                if(unif_rand()>=0)                {                    incidence[a*ncol+d]=1;incidence[d*ncol+a]=1;                    incidence[c*ncol+b]=1;incidence[b*ncol+c]=1;                    incidence[a*ncol+b]=0;incidence[b*ncol+a]=0;                    incidence[c*ncol+d]=0;incidence[d*ncol+c]=0;                    to[rand1]=d;                    to[rand2]=b;                    n++;                }                else                {                    incidence[a*ncol+c]=1;incidence[c*ncol+a]=1;                    incidence[d*ncol+b]=1;incidence[b*ncol+d]=1;                    incidence[a*ncol+b]=0;incidence[b*ncol+a]=0;                    incidence[c*ncol+d]=0;incidence[d*ncol+c]=0;                    //	from[rand1]=d;                    to[rand1]=c;                    from[rand2]=b;                    to[rand2]=d;                    n++;                }            }            else                if(incidence[a*ncol+d]==0 && incidence[c*ncol+b]==0)                {                    incidence[a*ncol+d]=1;incidence[d*ncol+a]=1;                    incidence[c*ncol+b]=1;incidence[b*ncol+c]=1;                    incidence[a*ncol+b]=0;incidence[b*ncol+a]=0;                    incidence[c*ncol+d]=0;incidence[d*ncol+c]=0;                    to[rand1]=d;                    to[rand2]=b;                    n++;                }//.........这里部分代码省略.........

void classForest(int *mdim, int *ntest, int *nclass, int *maxcat,        int *nrnodes, int *ntree, double *x, double *xbestsplit,        double *pid, double *cutoff, double *countts, int *treemap,        int *nodestatus, int *cat, int *nodeclass, int *jts,        int *jet, int *bestvar, int *node, int *treeSize,        int *keepPred, int *prox, double *proxMat, int *nodes) {    int j, n, n1, n2, idxNodes, offset1, offset2, *junk, ntie;    double crit, cmax;        zeroDouble(countts, *nclass * *ntest);    idxNodes = 0;    offset1 = 0;    offset2 = 0;    junk = NULL;    // Rprintf("nclass %d/n", *nclass);    for (j = 0; j < *ntree; ++j) {        // Rprintf("pCT nclass %d /n", *nclass);        /* predict by the j-th tree */		// Rprintf( "#ntree: %d, idxNodes: %d/n", j, idxNodes );        predictClassTree(x, *ntest, *mdim, nrnodes, treemap + 2*idxNodes,                nodestatus + idxNodes, xbestsplit + idxNodes,                bestvar + idxNodes, nodeclass + idxNodes,                treeSize[j], cat, *nclass,                jts + offset1, node + offset2, *maxcat);		//// original code		//predictClassTree(x, *ntest, *mdim, treemap + 2*idxNodes,		//	nodestatus + idxNodes, xbestsplit + idxNodes,		//	bestvar + idxNodes, nodeclass + idxNodes,		//	treeSize[j], cat, *nclass,		//	jts + offset1, node + offset2, *maxcat);                /* accumulate votes: */        for (n = 0; n < *ntest; ++n) {            countts[jts[n + offset1] - 1 + n * *nclass] += 1.0;        }                /* if desired, do proximities for this round */        if (*prox) computeProximity(proxMat, 0, node + offset2, junk, junk,                *ntest);        idxNodes += *nrnodes;        if (*keepPred) offset1 += *ntest;        if (*nodes)    offset2 += *ntest;    }        //Rprintf("ntest %d/n", *ntest);    /* Aggregated prediction is the class with the maximum votes/cutoff */    for (n = 0; n < *ntest; ++n) {        //Rprintf("Ap: ntest %d/n", *ntest);        cmax = 0.0;        ntie = 1;        for (j = 0; j < *nclass; ++j) {            crit = (countts[j + n * *nclass] / *ntree) / cutoff[j];            if (crit > cmax) {                jet[n] = j + 1;                cmax = crit;            }            /* Break ties at random: */            if (crit == cmax) {                ntie++;                if (unif_rand() > 1.0 / ntie) jet[n] = j + 1;            }        }    }        //Rprintf("ntest %d/n", *ntest);    /* if proximities requested, do the final adjustment     * (division by number of trees) */        //Rprintf("prox %d",*prox);    if (*prox) {        //Rprintf("prox: ntest %d/n", *ntest);        for (n1 = 0; n1 < *ntest; ++n1) {            for (n2 = n1 + 1; n2 < *ntest; ++n2) {                proxMat[n1 + n2 * *ntest] /= *ntree;                proxMat[n2 + n1 * *ntest] = proxMat[n1 + n2 * *ntest];            }            proxMat[n1 + n1 * *ntest] = 1.0;        }    }    //Rprintf("END ntest %d/n", *ntest);    }

//.........这里部分代码省略.........    if (*labelts) zeroDouble(yTestPred, ntest);    /* print header for running output */    if (*jprint <= *nTree) {	Rprintf("     |      Out-of-bag   ");	if (*testdat) Rprintf("|       Test set    ");	Rprintf("|/n");	Rprintf("Tree |      MSE  %%Var(y) ");	if (*testdat) Rprintf("|      MSE  %%Var(y) ");	Rprintf("|/n");    }    GetRNGstate();    /*************************************     * Start the loop over trees.     *************************************/    for (j = 0; j < *nTree; ++j) {    /* multinomial */    /*unsigned int coeffs[*sampsize];*/    /* for loop implementation */    /*double probs[*sampsize];*/    for (k = 0; k < *sampsize; ++k) {        probs[k] = 1/(*sampsize);    }    ran_multinomial(*sampsize,100,probs,coeffs);		idx = keepF ? j * *nrnodes : 0;		zeroInt(in, nsample);        zeroInt(varUsed, mdim);        /* Draw a random sample for growing a tree. */		if (*replace) { /* sampling with replacement */			for (n = 0; n < *sampsize; ++n) {				xrand = unif_rand();				k = xrand * nsample;				in[k] = 1;				yb[n] = y[k];				for(m = 0; m < mdim; ++m) {					xb[m + n * mdim] = x[m + k * mdim];				}			}		} else { /* sampling w/o replacement */			for (n = 0; n < nsample; ++n) nind[n] = n;			last = nsample - 1;			for (n = 0; n < *sampsize; ++n) {				ktmp = (int) (unif_rand() * (last+1));                k = nind[ktmp];                swapInt(nind[ktmp], nind[last]);				last--;				in[k] = 1;				yb[n] = y[k];				for(m = 0; m < mdim; ++m) {					xb[m + n * mdim] = x[m + k * mdim];				}			}		}		if (keepInbag) {			for (n = 0; n < nsample; ++n) inbag[n + j * nsample] = in[n];		}        /* grow the regression tree */		regTree(xb, yb, mdim, *sampsize, lDaughter + idx, rDaughter + idx,                upper + idx, avnode + idx, nodestatus + idx, *nrnodes,                treeSize + j, *nthsize, *mtry, mbest + idx, cat, tgini,                varUsed, coeffs);        /* predict the OOB data with the current tree */		/* ytr is the prediction on OOB data by the current tree */

//.........这里部分代码省略.........    REAL(sampleprobs)[m] = 1.0;    REAL(logmarg)[m] = logmargy;    REAL(shrinkage)[m] = shrinkage_m;    prior_m  = compute_prior_probs(model,pmodel,p, modelprior);    REAL(priorprobs)[m] = prior_m;    UNPROTECT(3);    old_loc = 0;    pmodel_old = pmodel;    nUnique=1;    INTEGER(counts)[0] = 0;    postold =  REAL(logmarg)[m] + log(REAL(priorprobs)[m]);    memcpy(modelold, model, sizeof(int)*p);  /*   Rprintf("model %d max logmarg %lf/n", m, REAL(logmarg)[m]); */    /*  Rprintf("Now Sample the Rest of the Models /n");  */  m = 0;  while (nUnique < k && m < INTEGER(BURNIN_Iterations)[0]) {    memcpy(model, modelold, sizeof(int)*p);    pmodel =  n_sure;    MH = 1.0;    if (pmodel_old == n_sure || pmodel_old == n_sure + n){	MH =  random_walk(model, vars,  n);	MH =  1.0 - problocal;    }    else {      if (unif_rand() < problocal) {      // random	MH =  random_switch(model, vars, n, pmodel_old, varin, varout );      }      else {      // Randomw walk proposal flip bit//	MH =  random_walk(model, vars,  n);      }    }    branch = tree;    newmodel= 0;    for (i = 0; i< n; i++) {      bit =  model[vars[i].index];      if (bit == 1) {	if (branch->one != NULL) branch = branch->one;	else newmodel = 1;	}      else {	if (branch->zero != NULL)  branch = branch->zero;	else newmodel = 1.0;      }      pmodel  += bit;    }    if (pmodel  == n_sure || pmodel == n + n_sure)  MH = 1.0/(1.0 - problocal);    if (newmodel == 1) {      new_loc = nUnique;      PROTECT(Rmodel_m = allocVector(INTSXP,pmodel));      model_m = INTEGER(Rmodel_m);

size_t analysis(unsigned short *incidence,size_t ncol,size_t nrow,double *scores,size_t step,size_t max_iter,size_t verbose,unsigned int seed){	set_rand(seed);	size_t i,j,kk,n,rand1,rand2;	size_t dim=max_iter+1;	size_t *from;  	size_t *to;  	unsigned short *matrix;	size_t a,b,c,d,e=0;	size_t index=1;  do matrix=(unsigned short *)calloc(nrow*ncol,sizeof(unsigned short)); while(matrix==NULL);	for(i=0;i<nrow;++i)	{		    for(j=0;j<ncol;++j)			{				matrix[i*ncol+j]=incidence[i*ncol+j];				e+=incidence[i*ncol+j];			}  }	//initialization of score vector overwriting the original	//do	scores=(double*)calloc(dim,sizeof(double));	while(scores==NULL);	for(i=0;i<dim;scores[i++]=0.0);	scores[0]=1.0; 	do	from=(size_t*)calloc(e,sizeof(size_t));	while(from==NULL);	do	to=(size_t*)calloc(e,sizeof(size_t));	  while(to==NULL); 	kk=0;	for(i=0;i<nrow;++i)		for(j=0;j<ncol;++j)			if(matrix[i*ncol+j]==1)				{					from[kk]=i;					to[kk]=j;					kk++;				}	time_t  tin,tfin;	tin = time (NULL);	//GetRNGstate();	for(n=0;n<max_iter;n++)	{		//random rewiring    if(verbose==1)  		loadBar( n,  max_iter, 100,  50);		rand1=(size_t) (unif_rand()*e);		do	rand2=(size_t) (unif_rand()*e);	while (rand1==rand2);		a=from[rand1];		c=from[rand2];		b=to[rand1];		d=to[rand2];		//printf("%d %d %d %d %d %d %d %d %d/n",e,a,b,c,d,rand1,rand2,t,n);		if(a!=c && d!=b && incidence[a*ncol+d]==0 && incidence[c*ncol+b]==0)			{				to[rand1]=d;				to[rand2]=b;				incidence[a*ncol+d]=incidence[c*ncol+b]=1;				incidence[a*ncol+b]=incidence[c*ncol+d]=0;			}			if(n%step==0)			(scores)[index++]=similarity(incidence,matrix, ncol,nrow, e);	}	tfin = time (NULL);	//PutRNGstate();	if(verbose==1)    printf("DONE in %d seconds /n",-(tin-tfin));	return (index-1);}