C++ NumericVector::size方法代码示例

NumericVector inclusionProb(const NumericVector& prob, const int& size) {
//	prob ... vector of initial probability weights
//	size ... sample size

	// check initial result
    int i;						// index for loops
    int nneg = 0;				// number of negative values
    int N = prob.size();		// number of observations
    double sum = 0.0;			// sum of probability weights
    NumericVector result(N);	// inclusion probabilities
    for(i = 0; i < N; i++) {
        if(prob[i] < 0.0) {
            result[i] = 0.0;	// set negative values to 0
            nneg++;  			// count negative values
        } else {
        	result[i] = prob[i];
        	if(prob[i] > 0.0) {
        		sum += prob[i];  // add current element to sum
        	}
        }
    }
    if(nneg > 0) warning("negative probability weights are set to 0");

    // compute inclusion probabilities
    if(sum > 0.0) {
        int ngeq1 = 0;	// number of values >= 1
        for(i = 0; i < N; i++) {
            if(result[i] > 0.0) {
                result[i] = result[i] * size / sum;	// initial adjustment
                if(result[i] >= 1.0) ngeq1++;  // count values >= 1
            }
        }
        // values >= 1 are set to 1 and the others are readjusted
        if(ngeq1 > 0) {
            int nset = 0;	// ???
        	// I think this results in an infinite loop
            while(ngeq1 != nset) {
                // compute sum of values less than 1
                sum = 0.0;	// reset sum
                for(i = 0; i < N; i++) {
                    if(result[i] > 0.0 && result[i] < 1.0) {
                        sum += result[i];  // add current element to sum
                    }
                }
                // readjust values
                if(sum > 0.0) {
                    for(i = 0; i < N; i++) {
                        if(result[i] > 0.0 && result[i] < 1.0) {
                            result[i] = result[i] * (size-ngeq1) / sum;
                        }
                    }
                }
                nset = ngeq1;
                // recount values >= 1 and set them to 1
                ngeq1 = 0;
                for(i = 0; i < N; i++) {
                    if(result[i] >= 1.0) {
                        result[i] = 1.0;
                        ngeq1++;  // count values >= 1
                    }
                }
            }
        }
    }

    return result;
}

// [[Rcpp::export]]
double sample_beta_scalar(NumericVector y, List x, NumericMatrix groups,
                   NumericMatrix beta, List gamma, NumericVector delta,
                   NumericMatrix z, double alpha, double sigma2,
                   List bs_beta, double s2_beta, int p_id, int k_id) {
  double out;
  NumericVector params(2);
  
  double sum_term = 0.0;
  double sum_term2 = 0.0;
  double sum_term2a = 0.0;
  double sum_term2b = 0.0;
  double sum_term3 = 0.0;
  double sum_term4 = 0.0;
  NumericMatrix x_tmp = x(k_id);
  NumericMatrix bs_beta_tmp = bs_beta(k_id);
  for (int i = 0; i < y.size(); i++) {
    double term = 0.0;
    for (int j = 0; j < x_tmp.ncol(); j++) {
      term += bs_beta_tmp(j, p_id) * x_tmp(i, j);
    }
    sum_term += term * term;
    sum_term2 += y[i] * term;
    sum_term2a += alpha * term;
    for (int k = 0; k < z.ncol(); k++) {
      sum_term2b += term * delta[k] * z(i, k);
    }
    for (int q = 0; q < gamma.size(); q++) {
      NumericVector gamma_q = gamma[q];
      int group_set = groups(i, q);
      sum_term3 += gamma_q[group_set] * term;
    }
    for (int k2 = 0; k2 < beta.nrow(); k2++) {
      NumericMatrix bs_beta_k = bs_beta[k2];
      NumericMatrix x_k = x[k2];
      if (k2 != k_id) {
        for (int s = 0; s < beta.ncol(); s++) {
          for (int j = 0; j < x_k.ncol(); j++) {
            sum_term4 += term * beta(k2, s) * bs_beta_k(j, s) * x_k(i, j);
          }
        }
      } else {
        for (int s = 0; s < beta.ncol(); s++) {
          if (s != p_id) {
            for (int j = 0; j < x_k.ncol(); j++) {
              sum_term4 += term * beta(k2, s) * bs_beta_k(j, s) * x_k(i, j);
            }
          }
        }
      }
    }
  }
  
  //Calculate params[1]
  params[1] = (sigma2 * s2_beta) / (s2_beta * sum_term + sigma2);
  
  //Calculate params[0]
  params[0] = params[1] * (-1.0 / (2.0 * sigma2)) * (-2.0 * sum_term2 + 2.0 * sum_term2a +
                                                     2.0 * sum_term2b + 2.0 * sum_term4 +
                                                     2.0 * sum_term3);

  //Calculate loglik
  out = rnorm(1, params[0], sqrt(params[1]))[0];
  return out;
}

// [[Rcpp::export]]
List naive_cluster(NumericVector D, IntegerVector K, double threshold) {
    int M = K.size();
    std::vector<std::vector<int> > result(M);
    std::vector<std::vector<double> > result_distances(M);
    std::vector<pair_dist> dst(D.size());
    std::map<int, std::set<topic_index> > clusters;

    // initialization
    int c = 0;
    for (int m = 0; m < M; ++m) {
        result[m].resize(K[m]);
        result_distances[m].resize(K[m]);
        for (int k = 0; k < K[m]; ++k) {
            result[m][k] = c;
            clusters[c].insert(topic_index(m, k));
            ++c;
        }
    }

    // brute force: we'll just write down which index in D
    // corresponds to which pair of topics, copying over D like space
    // is cheap or something. D is assumed to go block-by-block row-wise
    // along a block-upper-triangular matrix whose (I, J) block is the
    // K_I x K_J matrix of distances between topics in models I and J.
    // Within each block D runs rowwise along the block.
    // Zero blocks are omitted.
    int d = 0; // runs along D
    for (int m1 = 0; m1 < M - 1; ++m1) {
        for (int m2 = m1 + 1; m2 < M; ++m2) {
            for (int k1 = 0; k1 < K[m1]; ++k1) {
                for (int k2 = 0; k2 < K[m2]; ++k2) {
                    if (d >= dst.size()) {
                        stop("The length of D is not consistent with K");
                    }
                    dst[d].d = D[d];
                    dst[d].t1.first = m1;
                    dst[d].t1.second = k1;
                    dst[d].t2.first = m2;
                    dst[d].t2.second = k2;
                    d += 1;
                }
            }
        }
    }
    if (d != D.size()) {
        stop("The length of D is not consistent with K");
    }

    pair_dist_cmp pdc;
    std::sort(dst.begin(), dst.end(), pdc);

    int r1, r2;
    std::set<int> sect;
    bool allow;
    for (std::vector<pair_dist>::iterator d = dst.begin();
            d != dst.end(); ++d) {
        if (d->d > threshold) { // then we're done
            break;
        }

#ifdef __NAIVE_CLUSTER_LOG__
        Rcout << d->t1.first << " ";
        Rcout << d->t1.second << " | ";
        Rcout << d->t2.first << " ";
        Rcout << d->t2.second << " [";
        Rcout << d->d << "]";
#endif

        // get current cluster assignments
        r1 = result[d->t1.first][d->t1.second];
        r2 = result[d->t2.first][d->t2.second];
        if (r1 == r2) {
#ifdef __NAIVE_CLUSTER_LOG__
            Rcout << "...already merged." << std::endl;
#endif
            continue; // they're already merged
        }
        if (r1 > r2) {
            std::swap(r1, r2);  // guarantee r1 < r2
        }
        std::set<topic_index> &c1 = clusters[r1];
        std::set<topic_index> &c2 = clusters[r2];

        if (c1.size() == M || c2.size() == M) {
            // if either cluster is full, merge is impossible
#ifdef __NAIVE_CLUSTER_LOG__
            Rcout << "...cluster full." << std::endl;
#endif
            continue;
        }
        // Check whether the two clusters share topics from the same model
        sect.clear();
        allow = true;
#ifdef __NAIVE_CLUSTER_LOG__
        Rcout << " c1 " << "(" << r1 << "):";
#endif
        for (std::set<topic_index>::iterator t = c1.begin();
                t != c1.end(); ++t) {
#ifdef __NAIVE_CLUSTER_LOG__
//.........这里部分代码省略.........

// [[Rcpp::export]]
List SDP_single_foreq(int tmax, NumericVector res_bs, int cons_bs, int xc, NumericVector rep_gain, 
NumericMatrix f_m, NumericVector mort, List dec_ls, NumericVector rho_vec, NumericMatrix c_learn, 
NumericVector g_forage, NumericVector c_forage, List W_nr, List istar_nr, int N_init, int tsim, 
double alpha, double beta, double comp) {
   
   
   //Note:
   //Here, state variables are in standard format
   double xc_state = (double) xc; //state
   
   //To use a state variable X as an INDEX, it must be an INTEGER, and X=X-1
   xc = xc - 1; //Index
   
   //Book-keeping
   //Establish iniital variables required for the SDP
   int num_res = res_bs.size();
   //How many decisions? Count columns on first decision matrix
   NumericMatrix dec_ex = dec_ls(0);
   //int num_dec = dec_ex.ncol();
   int max_enc = f_m.nrow(); //Rows = encounters; Columns = resources
   
   //This is a state, not an index, so should equal cons_bs
   double xmax = (double) cons_bs;
   
   //Vectors for output
   IntegerVector pop_size(tsim);
   IntegerVector recruit_size(tsim);
   //OTHERS
   
   //Number of starting individuals
   int N = N_init;
   
   //Initialize the record-keeping of states
   //Record the state vector for each individual
   
   //Which of these individuals are alive (they all start out alive)
   IntegerVector alive(N);
   for (int n=0; n<N; n++) {
     alive(n) = 1;
   }
   int num_alive = sum(alive); //so num_alive = N at first
   
   //This records the energetic state of each individual... values will change over time
   IntegerVector state_v(N);
   for (int i=0; i<N; i++) {
     state_v(i) = xmax; //The initial state is full health
   }
   
   //Record the time vector for each individual
   IntegerVector time_v(N);
   for (int i=0; i<N; i++) {
     time_v(i) = 1; //tmax-1;
   }
   
   //Record the current resource for each individual
   IntegerVector res_v(N);
   for (int i=0; i<N; i++) {
    NumericVector rdraw_v = runif(1);
    double rdraw = as<double>(rdraw_v);
    //The initial resource is randomly drawn for each individual
    //Random draw is between 0 and num_res-1... so already 'indexed'
    res_v(i) = (int) floor(rdraw*(num_res)); 
   }
   
   
   //What are the optimal decisions for the current individual-level states?
   IntegerVector dec_v(N);
   for (int i=0; i<N; i++) {
     //res_v is already indexed
     //state_v and time_v are states, not indices, so alter by -1
     IntegerMatrix istar_t = istar_nr(res_v(i)); //Grab the istar for the focal resource
     dec_v(i) = istar_t(state_v(i)-1,time_v(i)-1); //Grab the decision for a given state and time t=0;
   }
   
   //TO EXPORT
   //Record initial pop size
   pop_size(0) = num_alive;
   recruit_size(0) = 0;
   
   List trophic_int(tsim);
   List energetic_state(tsim);
   List temporal_state(tsim);
   List decision_state(tsim);
   List fitness(tsim);
   
   trophic_int(0) = res_v;
   energetic_state(0) = state_v;
   temporal_state(0) = time_v;
   decision_state(0) = dec_v;
   
   NumericVector ind_fit(N);
   for (int i=0; i<N; i++) {
     NumericMatrix W_matrix = W_nr(res_v(i));
     ind_fit(i) = W_matrix(state_v(i)-1,time_v(i)-1);
   }
   fitness(0) = ind_fit;
   
   Rcpp::Rcout << "I got here " << num_alive << std::endl;
   
//.........这里部分代码省略.........

/* Function to run breadth-first search, returning only sets of connected 
   components that reside on the boundary of the districts */
List bsearch_boundary(List aList,
		      arma::vec boundary)
{

  /* Inputs to function:
     aList: adjacency list

     boundary: vector of boundary element indicators (as arma)
   */

  // Get indices of boundary units
  arma::uvec boundary_indices = find(boundary == 1);

  // Container - outputted of breadth search, a list
  List bsearch;

  // Container - partition vector, gets added to bsearch when queue is empty
  NumericVector partition;

  // Set mark vector - ledger of which indices have been reached
  NumericVector mark(aList.size());

  // Set queue vector
  NumericVector q;

  // Initialize breadth search with first element in boundary_indices
  mark(boundary_indices(0)) = boundary_indices(0);
  partition.push_back(boundary_indices(0));
  q = aList(boundary_indices(0));

  // Initialize objects inside loop
  int u; bool in_part; NumericVector adj_u; int i; int v; 

  // Begin do{} loop - run until number of elements in boundary_indices is 0
  do{

    // Begin while{} loop - run until q is empty
    while(q.size() > 0){
      
      // Dequeue first element in queue
      u = q(0);

      // Mark that element in ledger
      mark(u) = u;

      // Check if element is in the partition - add to partition if false
      in_part = is_true(any(partition == u));
      if(in_part == false){
	partition.push_back(u);
      }
      
      // Get adjacency vector for unit u
      adj_u = aList(u);

      // Loop through elements of adj_u, add to queue and mark if not reached
      if(adj_u.size() > 0){
	
	// Start loop
	for(i = 0; i < adj_u.size(); i++){
	  
	  // Reach element v
	  v = adj_u(i);

	  /* Check if already reached - if false, mark, add to partition, and
	     add to queue */
	  if(is_true(any(mark == v)) == FALSE){
	    mark(v) = v;
	    partition.push_back(v);
	    q.push_back(v);
	  }

	}

      }

      // Erase dequeued element from queue when done searching
      q.erase(q.begin());

    }

    // Handling an empty queue
    if(q.size() == 0){

      /* First, find boundary units that are in the reached partition and
	 remove them from boundary_units vector */
      for(i = boundary_indices.n_elem - 1; i >= 0; i--){
	if(is_true(any(partition == boundary_indices(i))) == TRUE){
	  boundary_indices.shed_row(i);
	}
      }
      
      // Store the partition, clear partition vector
      bsearch.push_back(partition);
      partition.erase(partition.begin(), partition.end());

      // Re-initialize breadth search from new starting value if nonempty
      if(boundary_indices.n_elem > 0){
	q = aList(boundary_indices(0));
//.........这里部分代码省略.........

// [[Rcpp::export]]
double g(NumericVector x) {
  return x.size();
}

// just an R interface to the coxfit routine, for regression testing purposes.
SEXP
cpp_coxfit(SEXP R_interface)
{
    // ---------------
    // get R objects:

    // extract survival times
    R_interface = CDR(R_interface);
    SEXP R_survTimes = CAR(R_interface);

    // censoring status
    R_interface = CDR(R_interface);
    SEXP R_censInd = CAR(R_interface);

    // offsets
    R_interface = CDR(R_interface);
    SEXP R_offsets = CAR(R_interface);

    // design matrix
    R_interface = CDR(R_interface);
    SEXP R_X= CAR(R_interface);

    // and the tie method
    R_interface = CDR(R_interface);
    SEXP R_method = CAR(R_interface);

    // ---------------
    // unpack R objects:

    // survival times
    const NumericVector n_survTimes = R_survTimes;
    //const AVector survTimes(n_survTimes.begin(), n_survTimes.size(),
    //                      false);
    
    // errors with const and stuff... what if we copy into new memory?
    const AVector survTimes(n_survTimes.begin(), n_survTimes.size());

    // censoring status
    const IntVector censInd = as<IntVector>(R_censInd);

    // offsets
    const AVector offsets = as<NumericVector>(R_offsets);

    // design matrix
    const NumericMatrix n_X = R_X;
    //const AMatrix X(n_X.begin(), n_X.nrow(),
    //                n_X.ncol(), false);
    
    //Same issue as above L717:721
    const AMatrix X(n_X.begin(), n_X.nrow(),
                    n_X.ncol());
    

    // tie method
    const int method = as<int>(R_method);

    // ---------------
    // assign remaining arguments for Coxfit

    const int nObs = survTimes.size();

    const AVector weights = arma::ones<AVector>(nObs);

    // ---------------

    // get new Coxfit object
    Coxfit cox(survTimes,
               censInd,
               X,
               weights,
               offsets,
               method);

    // do the fitting
    const int nIter = cox.fit();

    // get results: need to do that first!
    CoxfitResults fit = cox.finalizeAndGetResults();

    // check results
    cox.checkResults();

    // pack results into R list and return that
    return List::create(_["coef"] = fit.coefs,
                        _["imat"] = fit.imat,
                        _["loglik"] = fit.loglik,
                        _["nIter"] = nIter);

}

// [[Rcpp::export]]
NumericVector cumsum3(NumericVector x) {
  NumericVector res(x.size());
  std::partial_sum(x.begin(), x.end(), res.begin());
  return res;
}

//[[Rcpp::export]]
List  quasiTR(NumericVector start,
	      Function fn,
	      Function gr,
	      const List control) {
  
  
  using Eigen::VectorXi;
  using Eigen::Map;

  typedef MatrixXd optHessType;
  typedef LLT<optHessType> optPrecondType;
  
  int nvars = start.size();
  
  double rad = as<double>(control["start.trust.radius"]);
  const double min_rad = as<double>(control["stop.trust.radius"]);
  const double tol = as<double>(control["cg.tol"]);
  const double prec = as<double>(control["prec"]);
  const int report_freq = as<int>(control["report.freq"]);
  const int report_level = as<int>(control["report.level"]);
  const int report_precision = as<int>(control["report.precision"]);
  const int maxit = as<int>(control["maxit"]);
  const double contract_factor = as<double>(control["contract.factor"]);
  const double expand_factor = as<double>(control["expand.factor"]);
  const double contract_threshold = as<double>(control["contract.threshold"]);
  const double expand_threshold_rad = as<double>(control["expand.threshold.radius"]);
  const double expand_threshold_ap = as<double>(control["expand.threshold.ap"]);
  const double function_scale_factor = as<double>(control["function.scale.factor"]);
  const int precond_refresh_freq = as<int>(control["precond.refresh.freq"]);
  const int precond_ID = as<int>(control["preconditioner"]);
  const int quasi_newton_method = as<int>(control["quasi.newton.method"]);
  const int trust_iter = as<int>(control["trust.iter"]);
  
  std::string method_string;
  
  if (quasi_newton_method==1) {
    method_string = "SR1";
  } else {
    method_string = "BFGS";
  }
  
  Rfunc func(nvars, fn, gr);
  
  Map<VectorXd> startX(start.begin(),nvars); 
  
  // Control parameters for optimizer
  
  Trust_CG_Optimizer<Map<VectorXd>, Rfunc,
		     optHessType, optPrecondType> opt(func,
						      startX, rad, min_rad, tol,
						      prec, report_freq,
						      report_level,
						      report_precision,
						      maxit, contract_factor,
						      expand_factor,
						      contract_threshold,
						      expand_threshold_rad,
						      expand_threshold_ap,
						      function_scale_factor,
						      precond_refresh_freq,
						      precond_ID,
						      quasi_newton_method,
						      trust_iter);
  
  opt.run();
  
  // collect results and return
  
  VectorXd P(nvars);
  VectorXd grad(nvars);
  
  // return sparse hessian information in CSC format
  
  double fval, radius;
  int iterations;
  MB_Status status;
  
  status = opt.get_current_state(P, fval, grad,
				 iterations, radius);
  
  List res;
  res = List::create(Rcpp::Named("fval") = Rcpp::wrap(fval),
		     Rcpp::Named("solution") = Rcpp::wrap(P),
		     Rcpp::Named("gradient") = Rcpp::wrap(grad),	
		     Rcpp::Named("iterations") = Rcpp::wrap(iterations),
		     Rcpp::Named("status") = Rcpp::wrap((std::string) MB_strerror(status)),
		     Rcpp::Named("trust.radius") = Rcpp::wrap(radius),
		     Rcpp::Named("method") = Rcpp::wrap(method_string),
		     Rcpp::Named("hessian.update.method") = Rcpp::wrap(quasi_newton_method)
		     );
  
  return(res);
  
}

NumericVector sqrtCpp(const NumericVector & orig) {
  NumericVector vec(orig.size());     // create a target vector of the same size
  std::transform(orig.begin(), orig.end(), vec.begin(), sqrt_double);

  return(vec);
}

// [[Rcpp::export]]
List wasserstein_auto_(NumericVector p_, NumericVector q_, NumericMatrix cost_matrix_,
                  double epsilon, double desired_alpha){
  // compute distance between p and q
  // p corresponds to the weights of a N-sample
  // each q corresponds to the weights of a M-sample
  // Thus cost_matrix must be a N x M cost matrix
  // epsilon is a regularization parameter, equal to 1/lambda in some references
  int N = p_.size();
  int M = q_.size();
  
  Map<VectorXd> p(as<Map<VectorXd> >(p_));
  Map<VectorXd> q(as<Map<VectorXd> >(q_));
  Map<MatrixXd> cost_matrix(as<Map<MatrixXd> >(cost_matrix_));
  // avoid to take exp(k) when k is less than -500,
  // as K then contains zeros, and then the upcoming computations divide by zero
  MatrixXd K = (cost_matrix.array() * (-1./epsilon)); // K =  exp(- M / epsilon)
  for (int i = 0; i < N; i++){
    for (int j = 0; j < M; j++){
      if (K(i,j) < -500){
        K(i,j) = exp(-500);
      } else {
        K(i,j) = exp(K(i,j));
      }
    }
  }
  MatrixXd K_transpose = K.transpose();
  MatrixXd K_tilde = p.array().inverse().matrix().asDiagonal() * K; // diag(1/p) %*% K
  VectorXd u = VectorXd::Constant(N, 1./N);
  //
  VectorXd marginal1, marginal2;
  MatrixXd transportmatrix;
  VectorXd v;
  double alpha = 0;
  double beta = 0;
  int niterations_max = 1000;
  int iteration = 0;
  // for (int iteration = 0; iteration < niterations; iteration ++){
  while ((iteration < niterations_max) and (alpha < desired_alpha)){
    iteration ++;
    u = 1. / (K_tilde * (q.array() / (K_transpose * u).array()).matrix()).array();
    if (iteration % 10 == 1){
      // check if criterion is satisfied
      v = q.array() / (K_transpose * u).array();
      transportmatrix = u.col(0).asDiagonal() * K * v.col(0).asDiagonal();
      marginal1 = transportmatrix.rowwise().sum();
      marginal2 = transportmatrix.colwise().sum();
      alpha = 10;
      for (int i = 0; i < N; i++){
        beta = std::min(p(i) / marginal1(i), q(i) / marginal2(i));
        alpha = std::min(alpha, beta);
      }
      // cerr << "alpha = " << alpha << endl;
    }
  }
  v = q.array() / (K_transpose * u).array();
  // compute the optimal transport matrix between p and the first q
  transportmatrix = u.col(0).asDiagonal() * K * v.col(0).asDiagonal();
  // MatrixXd uXIv = u.array() * ((K.array() * cost_matrix.array()).matrix() * v).array();
  // NumericVector d = wrap(uXIv.colwise().sum());
  return List::create(Named("transportmatrix") = wrap(transportmatrix),
                      Named("u") = wrap(u),
                      Named("v") = wrap(v),
                      Named("iteration") = iteration);
}

// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
List inseq(mat M, bool adjust=true)
{
  int i, m, trun;
  mat mu=mean(M);
  //center the rows
  M.each_row() -= mu;
  int n=M.n_rows, p=M.n_cols;
  //Dtm is the vector of det(Sig)'s
  NumericVector Dtm;
  //gam_0 and gam_1 are for gam_2m and gam_2m+1, resp.
  mat gam0(p,p), gam1(p,p), Gam(p,p), Sig(p,p), Sig1(p,p), Gamadj(p,p), eigvec(p,p);
  //for adjustment, set initial increment in Gam=0
  Gamadj.zeros();
  //store the eigenvalues and eigenvectors of each Gam
  vec eigval(p),eigvalneg(p);
  double dtm;
  int sn= n/2; 
  for (m=0; m<n/2; m++)
  {
    gam0.zeros(); gam1.zeros();
    //calculate gam_2m (gam0) and gam_2m+1 (gam1)
    for(i=0; i<(n-2*m);i++) gam0+=trans(M.row(i))*M.row(i+2*m);
    gam0=gam0/n;
    for(i=0; i<(n-2*m-1);i++) gam1+=trans(M.row(i))*M.row(i+2*m+1);
    gam1=gam1/n;
    //Gam_m=gam_2m+gam_2m+1, then symmetrize
    Gam=gam0+gam1; Gam=(Gam+Gam.t())/2;
    
    if (m==0) Sig=-gam0+2*Gam;
    else Sig=Sig+2*Gam;
    
    if (eig_sym(Sig)(0)>0)
    {
      sn=m;
      break;
    }
  }
  if (sn>n/2-1) 
  {
    stop("Not enough samples.");
  }
  Dtm=det(Sig);
  for (m=sn+1; m<n/2; m++)
  {
    gam0.zeros(); gam1.zeros();
    //calculate gam_2m (gam0) and gam_2m+1 (gam1)
    for(i=0; i<(n-2*m);i++) gam0+=trans(M.row(i))*M.row(i+2*m);
    gam0=gam0/n;
    for(i=0; i<(n-2*m-1);i++) gam1+=trans(M.row(i))*M.row(i+2*m+1);
    gam1=gam1/n;
    //Gam_m=gam_2m+gam_2m+1, then symmetrize
    Gam=gam0+gam1; Gam=(Gam+Gam.t())/2;
    
    //Sig_m=Sig_m-1+2Gam_m
    Sig1=Sig+2*Gam;
    dtm=det(Sig1);
    //if dtm1>dtm, continue
    if (dtm<=Dtm(m-sn-1)) break;
    //update Sig
    Sig=Sig1;
    //record dtm
    Dtm.push_back(dtm);

    //to adjust the original Sig, subtract the negative part of Gam
    if (adjust) 
    {
      //calculate eigenvalues and eigenvectors of Gam
      eig_sym(eigval,eigvec,Gam);
      eigvalneg=eigval;
      eigvalneg.elem(find(eigvalneg>0)).zeros();
      Gamadj-=eigvec*diagmat(eigvalneg)*eigvec.t();
    }
  }
  trun = Dtm.size()-1+sn;
  List res; res["Sig"]=Sig; res["Dtm"]=Dtm; res["trunc"]= trun; res["sn"]=sn;
  if (adjust) res["Sigadj"]=Sig+2*Gamadj;
  return res;
}

// [[Rcpp::export]]
List rocUtilsPerfsAllC(NumericVector thresholds, NumericVector controls, NumericVector cases, std::string direction) {
  NumericVector se(thresholds.size());
  NumericVector sp(thresholds.size());
  long tp, tn;
  long i; // iterator over cases & controls
  
  if (direction == ">") {
    for (long t = 0; t < thresholds.size(); t++) {
      if (t % 100 == 0) Rcpp::checkUserInterrupt();
      double threshold = thresholds(t);
        tp = 0;
        for (i = 0; i < cases.size(); i++) {
          if (cases(i) <= threshold) {
            tp++;
        }
      }
      se(t) = (double)tp / cases.size();
      tn = 0;
      for (i = 0; i < controls.size(); i++) {
        if (controls(i) > threshold) {
          tn++;
        }
      }
      sp(t) = (double)tn / controls.size();
    }
  }
 else {
    for (long t = 0; t < thresholds.size(); t++) {
      if (t % 100 == 0) Rcpp::checkUserInterrupt();
      double threshold = thresholds(t);
      tp = 0;
      for (i = 0; i < cases.size(); i++) {
        if (cases(i) >= threshold) {
          tp++;
        }
      }
      se(t) = (double)tp / cases.size();
      long tn = 0;
      for (i = 0; i < controls.size(); i++) {
        if (controls(i) < threshold) {
          tn++;
        }
      }
      sp(t) = (double)tn / controls.size();
    }
  }
  List ret;
  ret["se"] = se;
  ret["sp"] = sp;
  return(ret);
}

SEXP
cpp_sampleGlm(SEXP r_interface)
{
    // ----------------------------------------------------------------------------------
    // extract arguments
    // ----------------------------------------------------------------------------------

    r_interface = CDR(r_interface);
    List rcpp_model(CAR(r_interface));

    r_interface = CDR(r_interface);
    List rcpp_data(CAR(r_interface));

    r_interface = CDR(r_interface);
    List rcpp_fpInfos(CAR(r_interface));

    r_interface = CDR(r_interface);
    List rcpp_ucInfos(CAR(r_interface));
    
    r_interface = CDR(r_interface);
    List rcpp_fixInfos(CAR(r_interface));
    
    r_interface = CDR(r_interface);
    List rcpp_distribution(CAR(r_interface));

    r_interface = CDR(r_interface);
    List rcpp_searchConfig(CAR(r_interface));

    r_interface = CDR(r_interface);
    List rcpp_options(CAR(r_interface));

    r_interface = CDR(r_interface);
    List rcpp_marginalz(CAR(r_interface));

    // ----------------------------------------------------------------------------------
    // unpack the R objects
    // ----------------------------------------------------------------------------------

    // data:
    const NumericMatrix n_x = rcpp_data["x"];
    const AMatrix x(n_x.begin(), n_x.nrow(),
                   n_x.ncol());

    const NumericMatrix n_xCentered = rcpp_data["xCentered"];
    const AMatrix xCentered(n_xCentered.begin(), n_xCentered.nrow(),
                           n_xCentered.ncol());

    const NumericVector n_y = rcpp_data["y"];
    const AVector y(n_y.begin(), n_y.size());

    const IntVector censInd = as<IntVector>(rcpp_data["censInd"]);

    // FP configuration:

    // vector of maximum fp degrees
    const PosIntVector fpmaxs = as<PosIntVector>(rcpp_fpInfos["fpmaxs"]);
    // corresponding vector of fp column indices
    const PosIntVector fppos = rcpp_fpInfos["fppos"];
    // corresponding vector of power set cardinalities
    const PosIntVector fpcards = rcpp_fpInfos["fpcards"];
    // names of fp terms
    const StrVector fpnames = rcpp_fpInfos["fpnames"];


    // UC configuration:

    const PosIntVector ucIndices = rcpp_ucInfos["ucIndices"];
    List rcpp_ucColList = rcpp_ucInfos["ucColList"];

    std::vector<PosIntVector> ucColList;
    for (R_len_t i = 0; i != rcpp_ucColList.length(); ++i)
    {
        ucColList.push_back(as<PosIntVector>(rcpp_ucColList[i]));
    }
    
    // fixed covariate configuration:
    
    const PosIntVector fixIndices = rcpp_fixInfos["fixIndices"];
    List rcpp_fixColList = rcpp_fixInfos["fixColList"];
    
    std::vector<PosIntVector> fixColList;
    for (R_len_t i = 0; i != rcpp_fixColList.length(); ++i)
    {
      fixColList.push_back(as<PosIntVector>(rcpp_fixColList[i]));
    }
    
    
    

    // distributions info:

    const double nullModelLogMargLik = as<double>(rcpp_distribution["nullModelLogMargLik"]);
    const double nullModelDeviance = as<double>(rcpp_distribution["nullModelDeviance"]);
    S4 rcpp_gPrior = rcpp_distribution["gPrior"];
    List rcpp_family = rcpp_distribution["family"];
    const bool tbf = as<bool>(rcpp_distribution["tbf"]);
    const bool doGlm = as<bool>(rcpp_distribution["doGlm"]);
    
    const double empiricalMean = as<double>(rcpp_distribution["yMean"]);
    const bool empiricalgPrior = as<bool>(rcpp_distribution["empiricalgPrior"]);
//.........这里部分代码省略.........

// [[Rcpp::export]]
List or_gdp_em(NumericVector ryo, NumericMatrix rxo, SEXP ralpha, SEXP reta){

	//Define Variables//
	int p=rxo.ncol();
	int no=rxo.nrow();
	double eta=Rcpp::as<double >(reta);
	double alpha=Rcpp::as<double >(ralpha);

	//Create Data//
	arma::mat xo(rxo.begin(), no, p, false);
	arma::colvec yo(ryo.begin(), ryo.size(), false);
	yo-=mean(yo);

	//Pre-Processing//
	Col<double> xoyo=xo.t()*yo;
	Col<double> B=xoyo/no;
	Col<double> Babs=abs(B);
	Mat<double> xoxo=xo.t()*xo;
	Mat<double> D=eye(p,p);
	Mat<double> Ip=eye(p,p);
	double yoyo=dot(yo,yo);
	double deltaB;
	double deltaphi;
	double phi=no/dot(yo-xo*B,yo-xo*B);
	double lp;

	//Create Trace Matrices
	Mat<double> B_trace(p,20000);
	Col<double> phi_trace(20000);
	Col<double> lpd_trace(20000);

	//Run EM Algorithm//
	cout << "Beginning EM Algorithm" << endl;
	int t=0;
	B_trace.col(t)=B;
	phi_trace(t)=phi;
	lpd_trace(t)=or_gdp_log_posterior_density(no,p,alpha,eta,yo,xo,B,phi);
	do{
		t=t+1;


		Babs=abs(B);
		D=diagmat(sqrt(((eta+sqrt(phi)*Babs)%(sqrt(phi)*Babs))/(alpha+1)));
		B=D*solve(D*xoxo*D+Ip,D*xoyo);

		phi=(no+p-3)/(yoyo-dot(xoyo,B));

		//Store Values//
		B_trace.col(t)=B;
		phi_trace(t)=phi;
		lpd_trace(t)=or_gdp_log_posterior_density(no,p,alpha,eta,yo,xo,B,phi);

		deltaB=dot(B_trace.col(t)-B_trace.col(t-1),B_trace.col(t)-B_trace.col(t-1));
		deltaphi=phi_trace(t)-phi_trace(t-1);
	} while((deltaB>0.00001 || deltaphi>0.00001) && t<19999);
	cout << "EM Algorithm Converged in " << t << " Iterations" << endl;

	//Resize Trace Matrices//
	B_trace.resize(p,t);
	phi_trace.resize(t);
	lpd_trace.resize(t);

	return Rcpp::List::create(
			Rcpp::Named("B") = B,
			Rcpp::Named("B_trace") = B_trace,
			Rcpp::Named("phi") = phi,
			Rcpp::Named("phi_trace") = phi_trace,
			Rcpp::Named("lpd_trace") = lpd_trace
			) ;

}