//.........这里部分代码省略.........
				freqg[j] = 0;
			for (j=0; j<*n; j++) { /* n times randomly select an index and increase its frequency */
				r = (unif_rand()**n);
				freqg[r] += 1;
			}
			/* now take corresponding number of observations from the data */
			r = 0;
			for (j=0; j<*n; j++) {
				for (k=0; k<freqg[j]; k++) { /* put the j-th observation freqg[j]-times to the bootstrap sample */
					eventg[r] = event[j];
					++r;
				}
			}
			for (first_event_g=0; (event[first_event_g]==0)&&(first_event_g<*n) ; first_event_g++)
				;
			for (j=0; j<*n; j++) { /* each observation is randomly assigned to one of the groups; group from input is reused */
				if (unif_rand() < ((double) y1/(y1+y2))) {
					group[j] = 1;
					--y1;
				} else {
					group[j] = 2;
					--y2;
				}
			}
			y1 = y1bak;
			y2 = y2bak;
			ks_cm_ad_process(eventg, group, n, &y1, &y2, rho, gamma, &first_event_g, du_temp, dsigma, ks_b, cm_w, cm_b, ad_w, ad_b);
			/* store this bootstrapped process for plotting */
			if (i<*nboot_plot) {
				for (j=0; j<*n; j++)
					du_boot_plot[j+i**n] = du_temp[j];
			}
			y1 = y1bak;
			y2 = y2bak;
			ks_cm_ad_stat(n,du_temp,g,ks_b,cm_w,cm_b,ad_w,ad_b,&temp1,&temp2,&temp3,&temp4,&temp5,&temp6);
			*pval_ks_w_boot += (temp1>*stat_ks_w);
			*pval_ks_b_boot += (temp2>*stat_ks_b);
			*pval_cm_w_boot += (temp3>*stat_cm_w);
			*pval_cm_b_boot += (temp4>*stat_cm_b);
			*pval_ad_w_boot += (temp5>*stat_ad_w);
			*pval_ad_b_boot += (temp6>*stat_ad_b);
		
		}
		*pval_ks_w_boot /= (double) *nboot;
		*pval_ks_b_boot /= (double) *nboot;
		*pval_cm_w_boot /= (double) *nboot;
		*pval_cm_b_boot /= (double) *nboot;
		*pval_ad_w_boot /= (double) *nboot;
		*pval_ad_b_boot /= (double) *nboot;
	}
	
	/* PERMUTATION tests (g is 1) */
	*pval_ks_w_perm = *pval_ks_b_perm = *pval_cm_w_perm = *pval_cm_b_perm = *pval_ad_w_perm = *pval_ad_b_perm = 0.;
	for (i=0; i<*n; i++)
		g[i] = 1.;
	if (*nperm>0) {
		for (i=0; i<*nperm; i++) {
			for (j=0; j<*n; j++) { /* each observation is randomly assigned to one of the groups; group from input is reused (destroyed) */
				if (unif_rand() < ((double) y1/(y1+y2))) {
					group[j] = 1;
					--y1;
				} else {
					group[j] = 2;
					--y2;
				}
			}
			y1 = y1bak;
			y2 = y2bak;
			ks_cm_ad_process(event, group, n, &y1, &y2, rho, gamma, &first_event, du_temp, dsigma, ks_b, cm_w, cm_b, ad_w, ad_b);
			/* store this permutation process for plotting */
			if (i<*nperm_plot) {
				for (j=0; j<*n; j++)
					du_perm_plot[j+i**n] = du_temp[j];
			}
			y1 = y1bak;
			y2 = y2bak;
			ks_cm_ad_stat(n,du_temp,g,ks_b,cm_w,cm_b,ad_w,ad_b,&temp1,&temp2,&temp3,&temp4,&temp5,&temp6);
			*pval_ks_w_perm += (temp1>*stat_ks_w);
			*pval_ks_b_perm += (temp2>*stat_ks_b);
			*pval_cm_w_perm += (temp3>*stat_cm_w);
			*pval_cm_b_perm += (temp4>*stat_cm_b);
			*pval_ad_w_perm += (temp5>*stat_ad_w);
			*pval_ad_b_perm += (temp6>*stat_ad_b);
		
		}
		*pval_ks_w_perm /= (double) *nperm;
		*pval_ks_b_perm /= (double) *nperm;
		*pval_cm_w_perm /= (double) *nperm;
		*pval_cm_b_perm /= (double) *nperm;
		*pval_ad_w_perm /= (double) *nperm;
		*pval_ad_b_perm /= (double) *nperm;
	}
	
	/* ASYMPTOTIC p-values of KS */
	*pval_ks_w_asympt = 1. - psupw(*stat_ks_w,50);
	*pval_ks_b_asympt = 1. - psupb(*stat_ks_b,.5,50);
	
	PutRNGstate();
	
}

void SimOneMVN_mxIW(double *nu, double *Lbdinvhlf, double *f1f2, int *pd, 
                    int *pnreps, int *pN, double *es, double *YY)
{
  int i, j, k, l, d, d2, N, nreps, mxnreps, Jrand;
  int *lbuff;

  double xd, sm, tstnu, nu_1, nu_2, zz, lambda, f_1, f_2;

  double *df, *pW, *SgmHlf, *Y, *xbuff, *Sigma, *sig, *SigInv;

  N = *pN;
  d = *pd;
  xd = (double) d;
  d2 = d*d;

  mxnreps=0;
  for(l=0;l<N;l++) if(mxnreps < *(pnreps+l)) mxnreps = *(pnreps+l);

  lbuff         = (int   *)S_alloc(        1,sizeof(int));

  df            = (double *)S_alloc(        1, sizeof(double));
  pW            = (double *)S_alloc(       d2, sizeof(double));
  xbuff         = (double *)S_alloc(        d, sizeof(double));
  SgmHlf        = (double *)S_alloc(       d2, sizeof(double));
  Y             = (double *)S_alloc(mxnreps*d, sizeof(double));
  Sigma         = (double *)S_alloc(       d2, sizeof(double));
  SigInv        = (double *)S_alloc(       d2, sizeof(double));
  sig           = (double *)S_alloc(        d, sizeof(double));

  f_1 = *f1f2;
  f_2 = *(f1f2+1);
  lambda = *nu/(2 * xd + 2.0);

  nu_1 = (2.0 * xd + 2.0)*(1.0 + (lambda - 1.0) *(1.0 - f_1)/(1-f_1/f_2));
  nu_2 = (2.0 * xd + 2.0)*(1.0 + (lambda - 1.0) *          f_2          );

  tstnu = f_1/(nu_2 - (2.0*xd+2.0)) + (1.0-f_1)/(nu_1 - (2.0*xd+2.0));
  tstnu = 1.0/tstnu + 2.0*xd + 2.0;
  Rprintf("nu_1=%g, nu_2=%g, nu ?= %g", nu_1, nu_2, tstnu);

  GetRNGstate();

  /* NOTE:                                                                              */
  /* this block computes the average std dev over genes from the model                  */
  /* its diagonal elements, passed to the pointer, sig (of size 3)                      */
  /* are used for the purposes of assigning mean value to Y's under the alternative     */
  /*                                                                                    */
  for(i=0;i<d;i++)                                                   
    for(j=0;j<d;j++){
      sm = 0.0;
      for(k=0;k<d;k++) 
        sm += *(Lbdinvhlf + d*i + k) * (*(Lbdinvhlf + d*j + k));
      *(SigInv + d*j + i) = sm;
    }
  matinv(SigInv, Sigma, pd);
  for(i=0;i<d;i++) *(sig + i) = pow((*(Sigma + d*i + i))/(*nu - 2.0*xd - 2.0), 0.5);

  for(l=0;l<N;l++){  

    /* Pick J = 1 w.p. 1-f_1, J= 2 w.p. f_1.                                            */
    /* Then draw an InvWish_d(nu_J, Lambda) matrix.  This is done                       */
    /* using the result:  if Sigma^(-1) ~ Wish_d(nu-d-1, Lambda^(-1)) then              */
    /* Sigma ~ InvWish_d(nu, Lambda).  I simulate N i.i.d. Wish_d(nu-d-1,Lambda^(-1))   */
    /* matrices and then invert to get Sigma.  One more catch, my rwishart routine      */
    /* uses the cholesky square root of the parameter matrix, Lambda instead of Lambda  */
    /* itself.  Since I want the parameter matrix in the Wishart to be Lambda^(-1) then */
    /* I should pass its cholesky square root which is Lbdinvhlf, e.g. the cholesky     */
    /* square root of Lambda inverse.  That is computed in the calling R script and     */
    /* passed in.  Notice the need to check that Lambda is nonsingular and that         */
    /* nu > 2*d + 2 (required so that the expected value of the inverse wishart         */
    /* is finite.)                                                                      */
    /*                                                                                  */
    zz = unif_rand();
    Jrand = 1*(zz < f_1);
    *df = (1-Jrand)*nu_1 + Jrand*nu_2 - xd - 1.0;
    rwishart1(df, pd, Lbdinvhlf, pW);
    matinv(pW, Sigma, pd);
    /*                                                                                  */
    /* Sigma ~ (1-f_1)*InvWish_d(nu_1, Lambda) + f_1*InvWish_d(nu_2, Lambda)            */
    /*                                                                                  */
    /* Next, use Sigma to simulate i.i.d. N(0_d, Sigma)'s                               */
    nreps = *(pnreps + l);
    *lbuff = nreps*d;
    rnormn(lbuff, Y); 
    chol(Sigma, SgmHlf, pd);

    for(i=0;i<nreps;i++){
      for(j=0;j<d;j++){
        sm = 0.0;
        for(k=0;k<d;k++) sm += (*(SgmHlf +d*j +k))*(*(Y +d*i +k));
        *(xbuff+j) = sm;
      }
      for(j=0;j<d;j++) *(Y + d*i + j) = *(xbuff + j) + *(es + l)*(*(sig + j));
    }
    for(i=0;i<(nreps*d);i++) *(YY + mxnreps*d*l + i) = *(Y+i);
  }
  PutRNGstate();

}

//.........这里部分代码省略.........
      double u, v;

      int card = (i != 0) * (i != halfM) + 2 * (j != 0) * (j != halfM);
      
      switch (card){
      case 3:
	//B(1) = {1}, B^c(1) = {2}
	j1 = (m - i) + m * j;
	j2 = i + m * (m - j);
	u = norm_rand();
	v = norm_rand();
	a[j1] = ia[j1] = M_SQRT1_2 * rho[j1];
	a[j1] *= u;
	ia[j1] *= v;
	a[j2] = ia[j2] = M_SQRT1_2 * rho[j2];
	a[j2] *= u;
	ia[j2] *= -v;
	
	//B(2) = {1,2}, B^c(2) = {0}
	j1 = (m - i) + m * (m - j);
	j2 = i + m * j;
	u = norm_rand();
	v = norm_rand();
	a[j1] = ia[j1] = M_SQRT1_2 * rho[j1];
	a[j1]*= u;
	ia[j1] *= v;
	a[j2] = ia[j2] = M_SQRT1_2 * rho[j2];
	a[j2]*= u;
	ia[j2] *= -v;      
	break;
      case 1:
	//B(1) = 0, B^c(1) = {1}
	j1 = i + m * j;
	j2 = m - i + m * j;
	u = norm_rand();
	v = norm_rand();
	a[j1] = ia[j1] = M_SQRT1_2 * rho[j1];
	a[j1] *= u;
	ia[j1] *= v;
	a[j2] = ia[j2] = M_SQRT1_2 * rho[j2];
	a[j2] *= u;
	ia[j2] *= -v;
	break;
      case 2:
	//B(1) = 0, B^c(1) = {2}
	j1 = i + m * j;
	j2 = i + m * (m - j);
	u = norm_rand();
	v = norm_rand();
	a[j1] = ia[j1] = M_SQRT1_2 * rho[j1];
	a[j1] *= u;
	ia[j1] *= v;
	a[j2] = ia[j2] = M_SQRT1_2 * rho[j2];
	a[j2] *= u;
	ia[j2] *= -v;
	break;
      case 0:
	j1 = i + m * j;
	a[j1] = rho[j1] * norm_rand();
	ia[j1] = 0;
	break;      
      }
    }

    /* ---------- Computation of Q \Lambda^1/2 Q* Z ------------ */
    int maxf, maxp, *iwork;
    double *work;
    
    /* The next lines is only valid for 2d random fields. I need to
       change if m_1 \neq m_2 as here I suppose that m_1 = m_2 = m */
    fft_factor(m, &maxf, &maxp);
    work = (double *)R_alloc(4 * maxf, sizeof(double));
    iwork = (int *)R_alloc(maxp, sizeof(int));
    fft_work(a, ia, m, m, 1, -1, work, iwork);
    
    fft_factor(m, &maxf, &maxp);
    work = (double *)R_alloc(4 * maxf, sizeof(double));
    iwork = (int *)R_alloc(maxp, sizeof(int));
    fft_work(a, ia, 1, m, m, -1, work, iwork);
        
    for (i=nbar;i--;)
      ans[i + k * nbar] = isqrtMbar * a[i % *ngrid + m * (i / *ngrid)];
  }
  PutRNGstate();  
    
  if (*nugget > 0){
    int dummy = *nsim * nbar;
    double sqrtNugget = sqrt(*nugget);
    
    GetRNGstate();
    for (i=dummy;i--;)
      ans[i] += sqrtNugget * norm_rand();

    PutRNGstate();
  }

  free(a); free(ia);

  return;
}

/*****************
 void MCMC_wrapper

 Wrapper for a call from R.

 and don't forget that tail -> head
*****************/
void MCMC_wrapper(int *dnumnets, int *nedges,
		  int *tails, int *heads,
		  int *dn, int *dflag, int *bipartite, 
		  int *nterms, char **funnames,
		  char **sonames, 
		  char **MHproposaltype, char **MHproposalpackage,
		  double *inputs, double *theta0, int *samplesize, 
		  double *sample, int *burnin, int *interval,  
		  int *newnetworktails, 
		  int *newnetworkheads, 
		  int *fVerbose, 
		  int *attribs, int *maxout, int *maxin, int *minout,
		  int *minin, int *condAllDegExact, int *attriblength, 
		  int *maxedges,
		  int *status){
  int directed_flag;
  Vertex n_nodes, nmax, bip;
  /* Edge n_networks; */
  Network nw[1];
  Model *m;
  MHproposal MH;
  
  n_nodes = (Vertex)*dn; 
  /* n_networks = (Edge)*dnumnets;  */
  nmax = (Edge)abs(*maxedges);
  bip = (Vertex)*bipartite; 
  
  GetRNGstate();  /* R function enabling uniform RNG */
  
  directed_flag = *dflag;

  m=ModelInitialize(*funnames, *sonames, &inputs, *nterms);

  /* Form the network */
  nw[0]=NetworkInitialize(tails, heads, nedges[0], 
                          n_nodes, directed_flag, bip, 0, 0, NULL);
  
  MH_init(&MH,
	  *MHproposaltype, *MHproposalpackage,
	  inputs,
	  *fVerbose,
	  nw, attribs, maxout, maxin, minout, minin,
	  *condAllDegExact, *attriblength);

  *status = MCMCSample(&MH,
		       theta0, sample, *samplesize,
		       *burnin, *interval,
		       *fVerbose, nmax, nw, m);
  
  MH_free(&MH);
        
/* Rprintf("Back! %d %d\n",nw[0].nedges, nmax); */

  /* record new generated network to pass back to R */
  if(*status == MCMC_OK && *maxedges>0 && newnetworktails && newnetworkheads)
    newnetworktails[0]=newnetworkheads[0]=EdgeTree2EdgeList(newnetworktails+1,newnetworkheads+1,nw,nmax-1);
  
  ModelDestroy(m);
  NetworkDestroy(nw);
  PutRNGstate();  /* Disable RNG before returning */
}

//.........这里部分代码省略.........
		error("'den' must be of type 'double'.");
    
	
	int *neighbors = INTEGER(sneighbors);
    int nneigh = asInteger(snneigh);
    if (nneigh <= 0)
		error("The number of neighbors must be positive.");
	int k = asInteger(sk);
	if (k <= 0)
		error("The number of components must be positive.");
	int *Z = INTEGER(sZ);
    int ldZ = LENGTH(sZ) / k;
    if (ldZ <= 0)
		error("The leading dimension of 'Z' must be positive.");
	int ldN = LENGTH(sneighbors) / nneigh;
	if (ldN <= 0)
		error("The leading dimension of 'neighbors' must be positive.");
	int ldD = LENGTH(sden) / k;
	if (ldD <= 0)
		error("The leading dimension of 'den' must be positive.");
	
	if (ldZ - 1 != ldN || ldZ - 1 != ldD || ldN != ldD)
		error("The leading dimension of 'Z', 'neighbors' and 'den' do not match.");
	
    int ldC = LENGTH(scheck) / k;
    double *den = REAL(sden);
    double *check = REAL(scheck);
    int nblocks = LENGTH(blocks);
    int b;
	
	
    GetRNGstate();
	
	
    for (b = 0; b < nblocks; b++) {
		SEXP spoints = VECTOR_ELT(blocks, b);
		int n = LENGTH(spoints);
		int *points = INTEGER(spoints);
		double *U = (double *) R_alloc(n, sizeof(double));
		/*double *prob = (double *) R_alloc(k, sizeof(double));*/
		/*  int *Ni = (int *) R_alloc(k, sizeof(int));*/
   
		int i;
             
		for (i = 0; i < n; i++)
			U[i] = unif_rand();
		

#pragma omp parallel for firstprivate(k, ldD, ldN, ldZ, ldC, nneigh, neighbors, \
				      Z, den, check, points, U)
		for (i = 0; i < n; i++) {
			int j, m;
			double s = 0.0;
			int number = 0;
			double prob[10];
			int Ni[10];
			/* compute the posterior weights for the different classes */
			for (j = 0; j < k; j++) {
				Ni[j] = 0;
				for (m = 0; m < nneigh; m++) {
					int mm = neighbors[points[i] - 1 + m * ldN] - 1;
					Ni[j] += Z[mm + j * ldZ];
				}
			}
			
            /*compute the index to use the check table*/            
			for (j = 0; j < k; j++) {
				number += Ni[j] * pow(nneigh+1, j);
			}
			
			
			for (j = 0; j < k; j++) {
				/* use tabled value check[or + j*ldC] = exp(beta * Ni[j]) */
				prob[j] = den[points[i] - 1 + j * ldD] * check[number + j * ldC];
				s += prob[j];
			}
			
			/* fix up the weights and convert to probabilities */
			if (s != 0.0 && R_FINITE(s))
				for (j = 0; j < k; j++)
					prob[j] = prob[j] / s;
			else
				for (j = 0; j < k; j++)
					prob[j] = 1.0 / k;
			
			/* generate new Z[points[i] - 1,] entries */
			for (j = 0; j < k; j++)
				Z[points[i] - 1 + j * ldZ] = 0;
			for (m = 0, s = 0.0, j = 0; j < k - 1; j++) {
				s += prob[j];
				if (U[i] > s) m++;
			}
			Z[points[i] - 1 + m * ldZ] = 1;
		}
    }
	
    PutRNGstate();
	
    return sZ;
}

void rgeomdirect(double *coord, int *nObs, int *nSite, int *dim,
		 int *covmod, int *grid, double *sigma2, double *nugget,
		 double *range, double *smooth, double *uBound,
		 double *ans){
  /* This function generates random fields for the geometric model

     coord: the coordinates of the locations
      nObs: the number of observations to be generated
    nSite: the number of locations
       dim: the random field is generated in R^dim
    covmod: the covariance model
      grid: Does coord specifies a grid?
    sigma2: the variance of the geometric gaussian process
      nugget: the nugget parameter
     range: the range parameter
    smooth: the smooth parameter
       ans: the generated random field */

  int i, j, neffSite, lagi = 1, lagj = 1, oneInt = 1;
  const double loguBound = log(*uBound), halfSigma2 = 0.5 * *sigma2;
  double sigma = sqrt(*sigma2), sill = 1 - *nugget;

  if (*grid){
    neffSite = R_pow_di(*nSite, *dim);
    lagi = neffSite;
  }

  else{
    neffSite = *nSite;
    lagj = *nObs;
  }

  double *covmat = malloc(neffSite * neffSite * sizeof(double)),
    *gp = malloc(neffSite * sizeof(double));

  buildcovmat(nSite, grid, covmod, coord, dim, nugget, &sill, range,
	      smooth, covmat);
  
  /* Compute the Cholesky decomposition of the covariance matrix */
  int info = 0;
  F77_CALL(dpotrf)("U", &neffSite, covmat, &neffSite, &info);

  if (info != 0)
    error("error code %d from Lapack routine '%s'", info, "dpotrf");

  GetRNGstate();
  
  for (i=*nObs;i--;){
    double poisson = 0;
    int nKO = neffSite;
    
    while (nKO) {
      /* The stopping rule is reached when nKO = 0 i.e. when each site
	 satisfies the condition in Eq. (8) of Schlather (2002) */
      poisson += exp_rand();
      double ipoisson = -log(poisson), thresh = loguBound + ipoisson;
	
      /* We simulate one realisation of a gaussian random field with
	 the required covariance function */
      for (j=neffSite;j--;)
	gp[j] = norm_rand();
      
      F77_CALL(dtrmv)("U", "T", "N", &neffSite, covmat, &neffSite, gp, &oneInt);
      
      nKO = neffSite;
      double ipoissonMinusHalfSigma2 = ipoisson - halfSigma2;
      for (j=neffSite;j--;){
	ans[j * lagj + i * lagi] = fmax2(sigma * gp[j] + ipoissonMinusHalfSigma2,
				      ans[j * lagj + i * lagi]);
	
	nKO -= (thresh <= ans[j * lagj + i * lagi]);
	  
      }
    }
  }
  
  PutRNGstate();

  /* So fare we generate a max-stable process with standard Gumbel
     margins. Switch to unit Frechet ones */
  for (i=*nObs * neffSite;i--;)
    ans[i] = exp(ans[i]);

  free(covmat); free(gp);

  return;
}

//.........这里部分代码省略.........
    case 1:
      whittleMatern(dist, mbar, zero, sill, *range, *smooth, rho);
      break;
    case 2:
      cauchy(dist, mbar, zero, sill, *range, *smooth, rho);
      break;
    case 3:
      powerExp(dist, mbar, zero, sill, *range, *smooth, rho);
      break;
    case 4:
      bessel(dist, mbar, *dim, zero, sill, *range, *smooth, rho);
      break;
    }

    /* Compute the eigen values to check if the circulant embbeding
       matrix is positive definite */

    /* Note : The next lines is only valid for 2d random fields. I
       need to change if there are m_1 \neq m_2 as I suppose that m_1
       = m_2 = m */
    int maxf, maxp;

    fft_factor(m, &maxf, &maxp);
    double *work = (double *)R_alloc(4 * maxf, sizeof(double));
    int *iwork = (int *)R_alloc(maxp, sizeof(int));
    fft_work(rho, irho, m, m, 1, -1, work, iwork);

    fft_factor(m, &maxf, &maxp);
    work = (double *)R_alloc(4 * maxf, sizeof(double));
    iwork = (int *)R_alloc(maxp, sizeof(int));
    fft_work(rho, irho, 1, m, m, -1, work, iwork);

    //Check if the eigenvalues are all positive
    for (i=mbar;i--;){
      notPosDef |= (rho[i] <= 0) || (fabs(irho[i]) > 0.001);
    }

    if (notPosDef){
      k++;
      m = HCN[k];
      halfM = m / 2;
      mbar = m * m;
    }

    if (k > 30)
      error("Impossible to embbed the covariance matrix");
    
  } while (notPosDef);
  /* --------- end of the embedding stage --------- */

  /* Computation of the square root of the eigenvalues */
  for (i=mbar;i--;){
    rho[i] = sqrt(rho[i]);
    irho[i] = 0;//No imaginary part
  }

  int mdag = m / 2 + 1, mdagbar = mdag * mdag;
  double isqrtMbar = 1 / sqrt(mbar);

  double *a = (double *)R_alloc(mbar, sizeof(double));
  double *ia = (double *)R_alloc(mbar, sizeof(double));
  
  GetRNGstate();
  for (i=*nObs;i--;){
    int nKO = nbar;
    double poisson = 0;
    
    while (nKO) {
      /* The stopping rule is reached when nKO = 0 i.e. when each site
	 satisfies the condition in Eq. (8) of Schlather (2002) */
      int j;
      double *gp = (double *)R_alloc(nbar, sizeof(double));
      
      poisson += exp_rand();
      double ipoisson = -log(poisson), thresh = loguBound + ipoisson;
      
      /* We simulate one realisation of a gaussian random field with
	 the required covariance function */
      circcore(rho, a, ia, m, halfM, mdag, mdagbar, *ngrid, nbar, isqrtMbar, *nugget, gp);
      
      nKO = nbar;
      double ipoissonMinusHalfSigma2 = ipoisson - halfSigma2;
      for (j=nbar;j--;){
	ans[j + i * nbar] = fmax2(sigma * gp[j] + ipoissonMinusHalfSigma2,
				  ans[j + i * nbar]);
	nKO -= (thresh <= ans[j + i * nbar]);
	
      }
    }
  }
  
  PutRNGstate();

  /* So fare we generate a max-stable process with standard Gumbel
     margins. Switch to unit Frechet ones */
  for (i=*nObs * nbar;i--;)
    ans[i] = exp(ans[i]);
  
  return;
}

 // [[Rcpp::register]]
 unsigned long exitRNGScope() {
     RNGScopeCounter--;
     if (RNGScopeCounter == 0) PutRNGstate();
     return RNGScopeCounter ;
 }

 void count_permus_with_at_least_k_unfixed_points ( int*n, int *k, double * res){
     GetRNGstate();
     Generic gen;
     *res = (double)gen.count_permus_with_at_least_k_unfixed_points(*n, *k);
     PutRNGstate();
 }

//.........这里部分代码省略.........
		// update beta
		for(j = 0; j < *q; j++){
			beta[j] = alphaLong[j]*ksi[j];
		}

		//update eta, eta+offset
		F77_CALL(dgemm)("N", "N", n, &oneInt, q, one, X, n, beta, q, zero, eta, n);
		for(int i=0; i<*n; i++) etaOffset[i] = eta[i] + offset[i];

		//update sigma_eps
		if(*family == 0){
			//resid = y - eta - offset
			F77_CALL(dcopy)(n, y, &nrv, resid, &nrv);  //resid <- y
			F77_CALL(daxpy)(n, minusOne, etaOffset, &nrv, resid, &nrv); //resid <- resid - eta - offset

			//rss = resid'resid
			rss = F77_CALL(ddot)(n, resid, &oneInt, resid, &oneInt);

			//update sigma2
			invSigma2 = rgamma(*n/2 + *b1, 1/(rss/2 + *b2));
			sqrtInvSigma2 = R_pow(invSigma2, 0.5);
			scale[0] = sqrtInvSigma2;
			*sigma2 = 1 / invSigma2;
		}


		if(i >= *burnin){
			/* report progress */
			if(*verbose){
				for(j=0; j<9; j++){
					if(i == pcts[j]){
						Rprintf(".");
						#ifdef Win32
							R_FlushConsole();
						#endif
						break;
					}
				}
			}
			/* save samples*/
			if(i == keep){
				for(j = 0; j < *q; j++){
					(betaMatR)[save + j*nSamp] = beta[j];
					(ksiMatR)[save + j*nSamp] = ksi[j];
				}
				for(j=0; j < *p; j++){
					(alphaMatR)[save + j*nSamp] = alpha[j];
				}
				for(j=0; j < *pPen; j++){
					(tau2MatR)[save + j*nSamp] = tau2[j];
					(gammaMatR)[save + j*nSamp] = gamma[j];
					(probV1MatR)[save + j*nSamp] = p1[j];
				}
				(wMatR)[save] = *w;
				(sigma2MatR)[save] = *sigma2;
				likMatR[save] = logLik(y, etaOffset, *family, scale, *n);
				(logPostMatR)[save] = updateLogPost(y, 	alpha, varAlpha,
						ksi, varKsi, scale, *b1, *b2, gamma, *w, *alphaW, *betaW,
						tau2, *a1, *a2,	*n, *q, *p, *pPen, pIncluded, qKsiNoUpdate, priorMeanKsi, *family, likMatR[save]);
				keep = keep + *thin;
				save ++;
				R_CheckUserInterrupt();
			}
		} else {
			if(*verbose){
				if(i == (*burnin-1)){
					Rprintf("b");
					#ifdef Win32
						R_FlushConsole();
					#endif
				}
			}
		}
	} /* end for i*/

	PutRNGstate();

	if(*verbose) Rprintf(".");
	if(*family > 0) {
		acceptance = 0.0;
		for(j=0; j<*blocksAlpha; j++) acceptance += acceptAlpha[j];
		acceptance = 0.0;
		if(qKsiNoUpdate < *q){
			for(j=0; j<*blocksKsi; j++) acceptance += acceptKsi[j];
		}
	}

	Free(etaOffset); Free(XKsiUpdate); Free(XAlpha);  Free(resid); Free(eta);
	Free(offsetKsi); Free(modeKsi); Free(priorMeanKsi);	Free(ksiUpdate);
	Free(offsetAlpha);
	Free(modeAlpha);
	Free(priorMeanAlpha);
	Free(varAlpha);
	Free(alphaLong);
	Free(penAlphaSq);
	freeXBlockQR(AlphaBlocks, *blocksAlpha);
	if(qKsiNoUpdate < *q) freeXBlockQR(KsiBlocks, *blocksKsi);
	Free(p1);
	return(R_NilValue);
}/* end sampler ()*/

//.........这里部分代码省略.........
  
        move = 70;     */
        
        if(move == 3)
        {
            BweibSurv_updateSC(beta, &alpha, &kappa, V, survTime, survEvent, cluster, survCov, c, d);
        }

  

        
        
        /* updating cluster-specific random effect: V 
         
         move = 130;*/
    
        if(move == 4)
        {
            BweibSurv_updateCP(beta, alpha, kappa, V, zeta, survTime, survEvent, cluster, survCov, n_j, accept_V, mhProp_V_var);
 
        }
        
        
        /* updating precition parameter of the random effect distribution: zeta
        
        move = 160;*/
        
        if(move == 5)
        {
            BweibSurv_updateVP(V, &zeta, rho1, rho2);
        }
        

        /*        */
        
        
        /* Storing posterior samples */
        
        
        if( ( (MM+1) % *thin ) == 0 && (MM+1) > (*numReps * *burninPerc))
        {
            StoreInx = (MM+1)/(*thin)- (*numReps * *burninPerc)/(*thin);
 
            samples_alpha[StoreInx - 1] = alpha;
            samples_kappa[StoreInx - 1] = kappa;

            if(*p >0)
            {
                for(j = 0; j < *p; j++) samples_beta[(StoreInx - 1) * (*p) + j] = gsl_vector_get(beta, j);
            }

            for(j = 0; j < *J; j++) samples_V[(StoreInx - 1) * (*J) + j] = gsl_vector_get(V, j);
            
            samples_zeta[StoreInx - 1] = zeta;


            if(MM == (*numReps - 1))
            {
                            /*                             */
                if(*p >0)
                {
                    for(j = 0; j < *p; j++) samples_misc[j] = (int) gsl_vector_get(accept_beta, j);
                }

                /*                 */
                samples_misc[*p]    = accept_alpha;
                /*    */
                for(i = 0; i < *J; i++) samples_misc[*p + 1 + i] = (int) gsl_vector_get(accept_V, i);

            }
            



        }
        
        if( ( (MM+1) % 10000 ) == 0)
        {
            
            time(&now);
            
            Rprintf("iteration: %d: %s\n", MM+1, ctime(&now));
            
            
            R_FlushConsole();
            R_ProcessEvents();
            
            
        }
        

        
    }

    
    PutRNGstate();
    return;
    
    
}

//.........这里部分代码省略.........
      case 2:			// fixed step
	dt = dtt;
	nstep = num_map_steps(t,time[step],dt);
	break;
      default:
	errorcall(R_NilValue,"unrecognized 'method'"); // # nocov
	break;
      }

      for (k = 0; k < nstep; k++) { // loop over Euler steps

	// interpolate the covar functions for the covariates
	table_lookup(&covariate_table,t,cp);

	for (j = 0, pm = ps, xm = xt; j < nreps; j++, pm += npars, xm += nvars) { // loop over replicates
	  
	  switch (mode) {

	  case Rfun: 		// R function

	    {
	      double *xp = REAL(xvec);
	      double *pp = REAL(pvec);
	      double *tp = REAL(tvec);
	      double *dtp = REAL(dtvec);
	      double *ap;
	      
	      *tp = t;
	      *dtp = dt;
	      memcpy(xp,xm,nvars*sizeof(double));
	      memcpy(pp,pm,npars*sizeof(double));
	      
	      if (first) {

	      	PROTECT(ans = eval(fcall,rho));	nprotect++; // evaluate the call
	      	if (LENGTH(ans) != nvars) {
	      	  errorcall(R_NilValue,"user 'step.fun' returns a vector of %d state variables but %d are expected: compare initial conditions?",
	      		LENGTH(ans),nvars);
	      	}
		
	      	PROTECT(nm = GET_NAMES(ans)); nprotect++;
	      	use_names = !isNull(nm);
	      	if (use_names) {
	      	  posn = INTEGER(PROTECT(matchnames(Snames,nm,"state variables"))); nprotect++;
	      	}

	      	ap = REAL(AS_NUMERIC(ans));
		
	      	first = 0;

	      } else {
	      
		ap = REAL(AS_NUMERIC(eval(fcall,rho)));

	      }
	      
	      if (use_names) {
	      	for (i = 0; i < nvars; i++) xm[posn[i]] = ap[i];
	      } else {
	      	for (i = 0; i < nvars; i++) xm[i] = ap[i];
	      }

	    }

	    break;
	      
	  case native: 		// native code

	    (*ff)(xm,pm,sidx,pidx,cidx,ncovars,cp,t,dt);

	    break;

	  default:

	    errorcall(R_NilValue,"unrecognized 'mode' %d",mode); // # nocov

	    break;

	  }

	}

	t += dt;
	
	if ((meth == 0) && (k == nstep-2)) { // penultimate step
	  dt = time[step]-t;
	  t = time[step]-dt;
	}
      }
    }
  }

  if (mode==1) {
    PutRNGstate();
    unset_pomp_userdata();
  }
  
  UNPROTECT(nprotect);
  return X;
}

//.........这里部分代码省略.........
                                    (nright[n] - nrightimp[n])) / nout[n];
                        }
                    }
                    if (noutall > 0) {
                        imprt[m + nclass*mdim] +=
                                ((double)(nrightall - nrightimpall)) / noutall;
                        impsd[m + nclass*mdim] +=
                                ((double) (nrightall - nrightimpall) *
                                (nrightall - nrightimpall)) / noutall;
                    }
                }
            }
        }
        
        /*  DO PROXIMITIES */
        if (iprox) {
            computeProximity(prox, oobprox, nodex, jin, oobpair, near);
            /* proximity for test data */
            if (*testdat) {
                computeProximity(proxts, 0, nodexts, jin, oobpair, ntest);
                /* Compute proximity between testset and training set. */
                for (n = 0; n < ntest; ++n) {
                    for (k = 0; k < near; ++k) {
                        if (nodexts[n] == nodex[k])
                            proxts[n + ntest * (k+ntest)] += 1.0;
                    }
                }
            }
        }
        
        if (keepf) idxByNnode += *nrnodes;
        if (keepInbag) idxByNsample += nsample0;
    }
    PutRNGstate();
   
    
    /*  Final processing of variable importance. */
    for (m = 0; m < mdim; m++) tgini[m] /= Ntree;
      
    if (imp) {
        for (m = 0; m < mdim; ++m) {
            if (localImp) { /* casewise measures */
                for (n = 0; n < nsample; ++n) impmat[m + n*mdim] /= out[n];
            }
            /* class-specific measures */
            for (k = 0; k < nclass; ++k) {
                av = imprt[m + k*mdim] / Ntree;
                impsd[m + k*mdim] =
                        sqrt(((impsd[m + k*mdim] / Ntree) - av*av) / Ntree);
                imprt[m + k*mdim] = av;
                /* imprt[m + k*mdim] = (se <= 0.0) ? -1000.0 - av : av / se; */
            }
            /* overall measures */
            av = imprt[m + nclass*mdim] / Ntree;
            impsd[m + nclass*mdim] =
                    sqrt(((impsd[m + nclass*mdim] / Ntree) - av*av) / Ntree);
            imprt[m + nclass*mdim] = av;
            imprt[m + (nclass+1)*mdim] = tgini[m];
        }
    } else {
        for (m = 0; m < mdim; ++m) imprt[m] = tgini[m];
    }
   
    /*  PROXIMITY DATA ++++++++++++++++++++++++++++++++*/
    if (iprox) {
        for (n = 0; n < near; ++n) {

void SummStats(Edge n_edges, Vertex n_nodes, Vertex *tails, Vertex *heads,
Network *nwp, Model *m, double *stats, Network *y0){

	Vertex *nodes; /*may need to change*/

  nodes = (Vertex *)malloc( 1 * sizeof(Vertex));

  GetRNGstate();  /* R function enabling uniform RNG */
  
  ShuffleEdges(tails,heads,n_edges); /* Shuffle edgelist. */
  
  for (unsigned int termi=0; termi < m->n_terms; termi++)
    m->termarray[termi].dstats = m->workspace;
  
  /* Doing this one toggle at a time saves a lot of toggles... */
  for(Edge e=0; e<n_edges; e++){
    ModelTerm *mtp = m->termarray;
    double *statspos=stats;
    
    nodes[0] = 0;

    for (unsigned int termi=0; termi < m->n_terms; termi++, mtp++){
      if(!mtp->s_func){
        (*(mtp->d_func))(1, tails+e, heads+e, nodes,
        mtp, nwp, y0);  /* Call d_??? function */
        for (unsigned int i=0; i < mtp->nstats; i++,statspos++)
          *statspos += mtp->dstats[i];
      }else statspos += mtp->nstats;
    }
    ToggleEdge(tails[e],heads[e],nwp);
  }
  

  /* Doing this one toggle at a time saves a lot of toggles... */
  for(int v=0; v<n_nodes; v++){
	  if(m->nodalstatus[v]== 2.0){

		  nodes[0] =  (int) (v+1);

		  ModelTerm *mtp = m->termarray;

		double *statspos=stats;

    for (unsigned int termi=0; termi < m->n_terms; termi++, mtp++){
      if(!mtp->s_func){
        (*(mtp->d_func))(1, tails, heads, nodes,
        mtp, nwp, y0);  /* Call d_??? function */
        for (unsigned int i=0; i < mtp->nstats; i++,statspos++)
          *statspos += mtp->dstats[i];
      }else statspos += mtp->nstats;
    }
    ToggleNode(v+1,nwp);
  }
  }



  ModelTerm *mtp = m->termarray;
  double *dstats = m->workspace;
  double *statspos=stats;
  for (unsigned int termi=0; termi < m->n_terms; termi++, dstats+=mtp->nstats, mtp++ ){
    if(mtp->s_func){
      (*(mtp->s_func))(mtp, nwp);  /* Call s_??? function */
      for (unsigned int i=0; i < mtp->nstats; i++,statspos++)
        *statspos = mtp->dstats[i];
    }else statspos += mtp->nstats;
  }
  
  PutRNGstate();
}

void simStahl(int *n_sim, double *nu, double *p, double *L,
              int *nxo, double *loc, int *max_nxo,
              int *n_bins4start)
{
    double **Loc, scale;
    double curloc=0.0, u;
    double *startprob, step;
    int i, j, n_nixo;

    /* re-organize loc as a doubly index array */
    Loc = (double **)R_alloc(*n_sim, sizeof(double *));
    Loc[0] = loc;
    for(i=1; i < *n_sim; i++)
        Loc[i] = Loc[i-1] + *max_nxo;

    GetRNGstate();

    if(fabs(*nu - 1.0) < 1e-8) { /* looks like a Poisson model */
        for(i=0; i< *n_sim; i++) {
            R_CheckUserInterrupt(); /* check for ^C */

            nxo[i] = rpois(*L);
            if(nxo[i] > *max_nxo)
                error("Exceeded maximum number of crossovers.");

            for(j=0; j < nxo[i]; j++)
                Loc[i][j] = runif(0.0, *L);
        }
    }
    else {
        scale = 1.0 / (2.0 * *nu * (1.0 - *p));

        /* set up starting distribution */
        startprob = (double *)R_alloc(*n_bins4start, sizeof(double));
        step = *L/(double)*n_bins4start;

        startprob[0] = 2.0*(1.0 - *p)*pgamma(0.5*step, *nu, scale, 0, 0)*step;
        for(i=1; i< *n_bins4start; i++) {
            R_CheckUserInterrupt(); /* check for ^C */

            startprob[i] = startprob[i-1] +
                2.0*(1.0 - *p)*pgamma(((double)i+0.5)*step, *nu, scale, 0, 0)*step;
        }

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

            nxo[i] = 0;

            /* locations of chiasmata from the gamma model */
            /* shape = nu, rate = 2*nu*(1-p) [scale = 1/{2*nu*(1-p)}] */

            u = unif_rand();
            if( u > startprob[*n_bins4start-1] )
                curloc = *L+1;
            else {
                for(j=0; j< *n_bins4start; j++) {
                    if(u <= startprob[j]) {
                        curloc = ((double)j+0.5)*step;
                        if(unif_rand() < 0.5) {
                            nxo[i] = 1;
                            Loc[i][0] = curloc;
                        }
                        break;
                    }
                }
            }

            if(curloc < *L) {
                while(curloc < *L) {
                    curloc += rgamma(*nu, scale);
                    if(curloc < *L && unif_rand() < 0.5) {
                        if(nxo[i] > *max_nxo)
                            error("Exceeded maximum number of crossovers.");

                        Loc[i][nxo[i]] = curloc;
                        (nxo[i])++;
                    }
                }
            }

            /* locations of crossovers from the no interference mechanism */
            if(*p > 0) {
                n_nixo = rpois(*L * *p);
                if(n_nixo + nxo[i] > *max_nxo)
                    error("Exceeded maximum number of crossovers.");

                for(j=0; j < n_nixo; j++)
                    Loc[i][nxo[i]+j] = runif(0.0, *L);
                nxo[i] += n_nixo;
            }
        }
    }

    /* sort the results */
    for(i=0; i< *n_sim; i++)
        R_rsort(Loc[i], nxo[i]);

    PutRNGstate();
}

 void random_permutation(int *n, int *sigma){
   GetRNGstate();
     Generic gen;
     gen.generate_random_permutation(*n, FIRST_ITEM, sigma);
     PutRNGstate();
 }

void dirichlet_wrapper(int *n, double *epsilon)
{
  GetRNGstate(); 
  rdirich(*n, epsilon);
  PutRNGstate();
}

//.........这里部分代码省略.........
  th = create_matrix(*d+1,*d+1);
  nuu = create_matrix(*d+1,*d+1);
  fam = create_intmatrix(*d+1,*d+1);
  // param in matrices:
  k = 0;
  for(i=0;i<((*d)-1);i++)
    {
      for(j=0;j<((*d)-i-1);j++)
	{
	  fam[i][j] = family[k];
	  nuu[i][j] = nu[k];
	  th[i][j] = par[k];
	  k++;
	  //printf("%d \t",fam[i][j]);
	}
      //printf("\n");
    }
  // Simulation
  GetRNGstate();

	/*
	Declare variable to hold seconds on clock.
*/
//time_t seconds;
/*
Get value from system clock and
place in seconds variable.
*/
//time(&seconds);
/*
Convert seconds to a unsigned
integer.
*/
//srand((unsigned int) seconds);


  // for i = 0
  uf[0][0] = u1[0];
  ub[0][0] = u1[0];
  // for i = 1,... d1-1
  // compute uf and ub
  for (int i = 1; i < (*d1); ++i)
    {
      uf[i][i] = u1[i];
      ub[i][i] = u1[i];
      for (int j = (i-1); j >= 0; --j)
	{
	  Hfunc(&fam[i-j-1][j],n,&ub[i][j+1], &uf[i-1][j],&th[i-j-1][j],&nuu[i-j-1][j],&ub[i][j]); //backward
  	  //printf("ub: %d,%d : %d, %5.2f : %10.8f   \n",i,j,fam[i-j-1][j], th[i-j-1][j], ub[i][j]);
	} 
      //printf("\n");
      for (int j = 0; j <= i-1; ++j)
	{
	  Hfunc(&fam[i-j-1][j],n, &uf[i-1][j], &ub[i][j+1],&th[i-j-1][j],&nuu[i-j-1][j],&uf[i][j]); //forward
  	  //printf("uf: %d,%d : %d, %5.2f : %10.8f   \n",i,j,fam[i-j-1][j], th[i-j-1][j], uf[i][j]);
	}
      //printf("\n");
    }
  // for  i= d1,.. d-1
  for (int i = (*d1); i < (*d); ++i)
    {
      // (a) Simulate uniform
      //out[i-(*d1)] =  rand()/(RAND_MAX+1.0);
	  out[i-(*d1)]=runif(0,1);
      // (b) inverse transformation:
      for (int j = 0; j < i; ++j)
	{
  	  //printf("inv: %d,%d : %d, %5.2f : %10.8f   \t",i-j-1,j,fam[i-j-1][j], th[i-j-1][j], uf[i-1][j]);
	  Hinv(&fam[i-j-1][j], n, &out[i-*d1] , &uf[i-1][j], &th[i-j-1][j], &nuu[i-j-1][j],&aux );
	  out[i-(*d1)]  = aux;
	  //printf("%10.8f   \n ", aux);
	}
      //printf("\n");
      if (i <((*d)-1))
	{
	  // forward and backward steps:
	  uf[i][i] = out[i-(*d1)];
	  ub[i][i] = out[i-(*d1)];
	  for (int j = i-1; j >= 0; --j)
	    {
	      Hfunc(&fam[i-j-1][j],n,&ub[i][j+1], &uf[i-1][j],&th[i-j-1][j],&nuu[i-j-1][j],&ub[i][j]); //backward
	      //printf("ub: %d,%d : %d, %5.2f : %10.8f   \n",i-j-1,j,fam[i-j-1][j], th[i-j-1][j], ub[i][j]);
	    } 
	  //printf("\n");
	  for (int j = 0; j <= i-1; ++j)
	    {
	      Hfunc(&fam[i-j-1][j],n, &uf[i-1][j], &ub[i][j+1],&th[i-j-1][j],&nuu[i-j-1][j],&uf[i][j]); //forward
	      //printf("uf: %d,%d : %d, %5.2f : %10.8f   \n",i-j-1,j,fam[i-j-1][j], th[i-j-1][j], uf[i][j]);
	    } 
	  //printf("\n");
	}
    }
  // free memory
  free_matrix(th,*d);    
  free_matrix(ub,*d);    
  free_matrix