C++ NumericVector类代码示例

// tmvrGaussian_density evaluates the PDF of the truncated multivariate Gaussian distribution with specified  bounds
// please note this function returns the un-logged PDF to avoid infinite values
// [[Rcpp::export]]
double tmvrGaussian_pdf(NumericVector x, NumericVector mu, NumericMatrix sigma, NumericVector lower, NumericVector upper){
  
  int n = x.size();
  double out;
  
  //generate boolean vector indicating if points in x fall inside support region
  LogicalVector support(n);
  for(int i=0; i < n; i++){
    if((x(i) >= lower[i]) && (x(i) <= upper[i])){
      support[i] = true;
    } else {
      support[i] = false;
    }
  }
  
  //check if any points fell outside the support region
  bool outside_support = false;
  for(int i=0; i < support.size(); i++){
    if(support[i] == false){
      outside_support = true;
    }
  }
  
  //if points fall outside support region immediately break and return 0
  if(outside_support){
    out = 0;
    return(out);
  }
  
  //if points fall inside support region calculate the density
  Environment mvtnorm("package:mvtnorm");
  Function pmvnorm = mvtnorm["pmvnorm"];
  double pdf;
  SEXP cdf_r;
  pdf = mvrGaussian_pdf(as<arma::colvec>(x),as<arma::colvec>(mu),as<arma::mat>(sigma));
  cdf_r = pmvnorm(lower,upper,mu,R_NilValue,sigma);
  double cdf = as<double>(cdf_r);
  out = pdf / cdf;
  
  return(out);
}

NumericVector drawHiddenVarsSample(const NumericVector &O_vec,  double time, const NumericVector &lambda,  const NumericVector &topo_path, 
                                   const vector< vector<int> > &parents, double &dens, double lambda_s, bool sampling_time_available) {
  int p = lambda.length();
  NumericVector T(p), T_sum(p);
  dens = 0.0;
  
  if(sampling_time_available == false) {
    time = my_rexp(1, lambda_s);
  }
  
  for(int k = 0; k < p; k++) {
    int e = topo_path[k];  
    
    double max_parents_t = 0.0;
    
    for(unsigned int j = 0; j < parents[e].size(); j++) {
      if(T_sum[ parents[e][j] ] > max_parents_t) {
        max_parents_t = T_sum[parents[e][j]];
      }
    }
    
    if( O_vec[e] == 1 ) {
      double tmp = rtexp(lambda[e],  time-max_parents_t);
      
      dens = dens +   my_dexp( tmp, lambda[e], true) - my_pexp(time-max_parents_t, lambda[e], true);
      
      T_sum[e] =  max_parents_t + (tmp);
    } else {
      
      double  tmp = my_rexp(1, lambda[e]) ;
      T_sum[e] =  std::max(time, max_parents_t) + tmp;    
      
      double tmp2 = my_dexp( tmp, lambda[e], true) ;
      dens +=  tmp2;
    }
    
    T[e] = T_sum[e] - max_parents_t ;      
  }
  return(T);
  //  return Rcpp::List::create(Rcpp::Named("T") = T, Rcpp::Named("density")=dens);
}

//data is assumed to be sorted by time
// [[Rcpp::export]]
NumericVector FindIntervalCalibCPPvec(NumericVector w, NumericVector wres) {
  int npoints = w.size();
  NumericVector out(2);
  bool CondLocation;
      int j=0;
      CondLocation = true;
      while(CondLocation == true)
      {
      if (wres[j] == 1)
      {
        if(j==0)
        {
          out[0] = 0;
          out[1] = w(0);
          CondLocation = false;
        } else
          {     
          out[0] = w[j-1];
          out[1] = w[j];
          CondLocation = false;
          }
      } else if (wres[j]==INFINITY)
      {
        if (j==0)
        {
          out[0] = 0;
          out[1] = INFINITY;
        } else {
        out[0] = w[j-1];
        out[1] = INFINITY;
        }
        CondLocation = false;
      } else if (j==npoints-1) {
        out[0] = w[j];
        out[1] = INFINITY;
        CondLocation = false;}
      j += 1;
      }
        
    return out;
  }

double corRcpp(NumericVector x, const NumericVector y)
{
    size_t n = x.size();
    double ex(0), ey(0), sxx(0), sxy(0), syy(0), xt(0), yt(0);
    double tiny = 1e-20;

  for (size_t i = 0; i < n; i++) { // Find the means.
   ex += x[i];
   ey += y[i];
  }
  ex /= n;
  ey /= n;
  for (size_t i = 0; i < n; i++) { // Compute the correlation coefficient.
    xt = x[i] - ex;
    yt = y[i] - ey;
    sxx += xt * xt;
    syy += yt * yt;
    sxy += xt * yt;
  }
return sxy/(sqrt(sxx*syy)+tiny);
}

// [[Rcpp::export]]
List maxcpp(NumericVector tau, int x, int y ,int j) {
    arma::cube mycube(tau.begin(), x, y, j);
    arma::vec sub(j);
    arma::mat max(x,y), max_idx(x,y);
    for (int row = 0; row < x; row++) {
        for (int col = 0; col < y; col++) {
            for (int drug = 0; drug < j; drug++) {
                sub(drug) = mycube(row, col, drug);
            }
            // get the index for max and min value
            arma::uword maxidx;
            // get max and min by slice
            max(row, col) = sub.max(maxidx);
            max_idx(row, col) = maxidx+1;
        }
    }
    return Rcpp::List::create(
               Rcpp::Named("max") = max,
               Rcpp::Named("max.idx") = max_idx
           ) ;
}

// [[Rcpp::export]]
List rwmhUpdate(NumericVector x, NumericVector eps, F<double> f){

  int n = x.size();
  NumericVector y = x + eps;
  double p1= f(x);
  double p2 = f(y);
  double temp = exp(p2-p1);
  double accept = std::min(1.0,temp);
  double u = R::runif(0,1);
  if (u < accept) {
    return List::create(
      _["chain"] = y,
      _["acc_rate"] = accept
    );
    }else{
      return List::create(
        _["chain"] = x,
        _["acc_rate"] = accept
      );
      }
}

float closeness(NumericVector x1, NumericVector x2) {
	float y = 0;
	int y_n = 0;
	int n = x1.size();

	for(int i = 0; i < n; i ++) {
		for(int j = 0; j < n; j ++) {
			if(std::abs(x1[i] - 0) < 1e-6 | std::abs(x2[j] - 0) < 1e-6) {
				continue;
			}
			y += std::abs(i - j);
			y_n += 1;
		}
	}

	if(y_n == 0) {
		return 0;
	} else {
		return y/y_n;
	}
}

// [[Rcpp::export("pcevKernel")]]
SEXP pcevKernelRcpp(NumericMatrix Yr, NumericMatrix Zr, NumericVector eiValuer) {
   int n = Yr.nrow(), p = Yr.ncol(), q = Zr.ncol();

   arma::mat Y(Yr.begin(), n, p, false);       // reuses memory and avoids extra copy
   arma::mat Z(Zr.begin(), n, q, false);
   arma::colvec eiValue(eiValuer.begin(), n, false);

   // Initialize parameters
   arma::mat F(n, p, fill::zeros);
   arma::mat E(n, p, fill::zeros);
   theta param_new(fastLm(Y, Z), F, E);

   // Make copy and update
   theta param_old = param_new;
   param_new.update(Y, Z, eiValue);
   
   return Rcpp::List::create(
   Rcpp::Named("LogLik") = param_new.logLL(eiValue),
   Rcpp::Named("tau") = param_new.getTau(),
   Rcpp::Named("Sigma") = param_new.getSigma());
}

void
Assembly::addCachedResidual(NumericVector<Number> & residual, Moose::KernelType type)
{
  std::vector<Real> & cached_residual_values = _cached_residual_values[type];
  std::vector<unsigned int> & cached_residual_rows = _cached_residual_rows[type];

  mooseAssert(cached_residual_values.size() == cached_residual_rows.size(), "Number of cached residuals and number of rows must match!");

  residual.add_vector(cached_residual_values, cached_residual_rows);

  if (_max_cached_residuals < cached_residual_values.size())
    _max_cached_residuals = cached_residual_values.size();

  // Try to be more efficient from now on
  // The 2 is just a fudge factor to keep us from having to grow the vector during assembly
  cached_residual_values.clear();
  cached_residual_values.reserve(_max_cached_residuals*2);

  cached_residual_rows.clear();
  cached_residual_rows.reserve(_max_cached_residuals*2);
}

// [[Rcpp::export]]
List PredictiveRecursion_DifferentSigma(NumericVector z, double mu0, NumericVector sig0, IntegerVector sweeporder,
    NumericVector grid_x, NumericVector theta_guess,
    double nullprob=0.95, double decay = -0.67) {
  // z: vector of observed data, vector of length N
  // sig0: vector of standard errors se(z_i) under the null, of same length as z
  // sweeporder: a vector of indices in [0...N-1] representing the order in which the z are processed
  //     length(sweeporder) = Npasses*length(z)
  // grid_x: a grid of points at which the alternative density will be approximated
  // theta_guess: an initial guess for the sub-density under the alternative hypothesis
  // nullprob: an initial guess for the fraction of null cases
  // decay: the stochastic-approximation decay parameter, should be in (-1, -2/3)

  // Set-up
  int n = sweeporder.size();
  int k, gridsize = grid_x.size();
  NumericVector theta_subdens(clone(theta_guess));
  double pi0 = nullprob;
  NumericVector joint1(gridsize);
  NumericVector ftheta1(gridsize);
  double m0, m1, mmix, cc;

  // Begin sweep through the data
  for(int i=0; i<n; i++) {
    if(i % 200 == 0) Rcpp::checkUserInterrupt();  
    k = sweeporder[i];
    cc = pow(3.0+(double)i, decay);
    joint1 = dnorm(grid_x, z[k] - mu0, sig0[k]) * theta_subdens;
    m0 = pi0*R::dnorm(z[k] - mu0, 0.0, sig0[k], 0);
    m1 = trapezoid(grid_x, joint1);
    mmix = m0 + m1;
    pi0 = (1.0-cc)*pi0 + cc*(m0/mmix);
    ftheta1 = joint1/mmix;
    theta_subdens = (1.0-cc)*theta_subdens + cc*ftheta1;
  }

  return Rcpp::List::create(Rcpp::Named("grid_x")=grid_x,
          Rcpp::Named("theta_subdens")=theta_subdens,
          Rcpp::Named("pi0")=pi0
          );
}

//' Computes the convex minorant of a vector.
//' @param x,y Vector of x and y values
//' @return x.knots, y.knots, y.slopes and the left derivative at all x values
//' @export
//[[Rcpp::export]]
List GreatestConvexMinorant(NumericVector x, NumericVector y) {
  
  int ny = y.length();
  
  NumericVector XX = x; 
  NumericVector XY = y; 
  NumericVector leftDerivative(ny - 1);
  
  vector<Point> P(ny);
  for (int i = 0; i < ny; i++) {
    P[i].x = XX[i];
    P[i].y = XY[i];
  }
  
  vector<Point> convHull = convex_hull(P);
  
  int            nP = convHull.size();
  vector<double> convHullX(nP); 
  vector<double> convHullY(nP); 
  for (int i = 0; i < nP; i++) {
    convHullX[i] = convHull.at(i).x;
    convHullY[i] = convHull.at(i).y;
  }
  
  NumericVector XXX     = Rcpp::wrap(convHullX); // correct
  NumericVector XYY     = Rcpp::wrap(convHullY); // correct
  NumericVector Slopes  = diff(XYY) / diff(XXX); // has to be corrected
  //return slopes;
  
  for (int i = 0; i < ny-1; i++) {
    for (int j = 0; j < nP; j++) {
      if (XXX[j] < XX[i+1]) leftDerivative[i] = Slopes[j]; 
    }
  }
  
  return List::create(Named("y.slopes") = Slopes,
                      Named("x.knots") = XXX,
                      Named("y.knots") = XYY,
                      Named("left.derivative") = leftDerivative);
}

// [[Rcpp::export]]
NumericVector eloDstC(NumericVector rating, IntegerVector raceCount, NumericVector time, double K, double P, int provisionalN) {
  int n = rating.size();
  NumericVector out = clone(rating);
  double W_e, e, pct, R, delta_a, delta_b;
  
  for (int j1 = 0; j1 < n-1; ++j1){
    for (int j2 = j1 + 1; j2 < n; ++j2){
      e = -1 * (rating[j1] - rating[j2]) / 400;
      W_e = 1 / (pow(10,e) + 1);
      
      pct = 100 * (time[j2] - time[j1]) / time[j1];
      if (pct <= 1){
        R = pct;
      }
      else{
        R = pow(pct,0.25);
      }
      
      delta_a = (K / (n-1)) * R * (1 - W_e);
      delta_b = -1 * delta_a;
      
      if (raceCount[j1] < provisionalN && raceCount[j2] < provisionalN){
        delta_a = P * delta_a;
        delta_b = P * delta_b;
      }
      if (raceCount[j1] < provisionalN && raceCount[j2] >= provisionalN){
        delta_a = P * delta_a;
        delta_b = 0;
      }
      if (raceCount[j1] >= provisionalN && raceCount[j2] < provisionalN){
        delta_a = 0;
        delta_b = P * delta_b;
      }
      
      out[j1] += delta_a;
      out[j2] += delta_b;
    }
  }
  return out;
}

// [[Rcpp::export]]
NumericMatrix Kest_anin_border_c(NumericMatrix coord, NumericVector lambda, 
                                 NumericMatrix bbox, NumericVector bdist, 
                                 NumericVector r, NumericMatrix directions,
                                 double epsilon, int border=1) {
  
  Pp pp(coord, lambda, bbox); 
  
  int nr = r.size();
  
  int dim = directions.ncol();
  int ndir = directions.nrow();
  NumericMatrix out(nr, ndir);
  
  int i,j,l, ui, ri;
  double d, w, dot, ang;
  double rmax = r(nr-1);
  
  for(i=0; i < pp.size(); i++) {
    for(j= 0; j < pp.size(); j++) {
      if(i!=j){
        d = pp.getDist(&i, &j); 
        if(d < rmax) {
          w = 1 / (pp.getMark(&i) * pp.getMark(&j) );
          for(ui=0; ui < ndir; ui++){
            dot = 0;
            for(l=0; l < dim; l++)  dot += (pp.getCoord(&j,&l)-pp.getCoord(&i,&l)) * directions(ui, l);
            ang = acos(dot/d);
            ang = fmin(ang, PI-ang);
            if(ang < epsilon){
              for(ri=0; ri < nr; ri++){
                if(d < r(ri) && (bdist(i) > r(ri) | border==0) ) out(ri, ui) += w;
              }
            }
          }
        }
      }
    }
  }
  return out;
}

// [[Rcpp::export]]
IntegerMatrix raceSimC(NumericVector ratings,int nsim) {
   int nr = ratings.size();
   double p;
   NumericVector r;
   IntegerMatrix out(nr,nsim);
   std::fill( out.begin(), out.end(), 1);
   
   for (int i = 0; i < nsim; ++i){
     for (int j = 0; j < nr-1; ++j){
       for (int k = j+1; k < nr; ++k){
         p = 1 / (pow(10,-(ratings[j] - ratings[k]) / 400) + 1);
         r = runif(1);
         if (r[0] <= p){
           out(k,i) += 1;
         }else{
           out(j,i) += 1;
         }
       }
     }
   }
   return out;
}

// **********************************************************//
//                    Cache time intervals                   //
// **********************************************************//
// [[Rcpp::export]]
List timefinder (NumericVector timestamps, IntegerVector edgetrim, double timeunit) {
    int D = timestamps.size();
    List out(D);
    for (int d = min(edgetrim)-1; d < D; d++) {
		double time0 = timestamps[d-1];
		double time1 = time0-24*timeunit;
		double time2 = time0-96*timeunit;
		double time3 = time0-384*timeunit;
		int id1 = which_num(time1, timestamps);
		int id2 = which_num(time2, timestamps);
		int id3 = which_num(time3, timestamps);
		arma::mat intervals(3, 2);
		intervals(0,0) = id1;
		intervals(0,1) = d;
		intervals(1,0) = id2;
		intervals(1,1) = id1-1;
		intervals(2,0) = id3;
		intervals(2,1) = id2-1;
		out[d] = intervals;		
	}
    return out;
}