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

// [[Rcpp::export]]
Rcpp::List calib( arma::mat Y, 
                  arma::vec C,
                  arma::mat Z,  
                  NumericVector mu_input,
                  IntegerVector mu_dim,
                  NumericVector mu0_input,
                  IntegerVector mu0_dim
                )
{
  
  arma::cube mu(mu_input.begin(), mu_dim[0], mu_dim[1], mu_dim[2]);
  arma::cube mu0(mu0_input.begin(), mu0_dim[0], mu0_dim[1], mu0_dim[2]);
  
  int n = Y.n_rows;
  int p = Y.n_cols;
  int niter = Z.n_rows;
  cube calibration(niter,p,n); calibration.fill(0); 
  mat calibrationMedian(n,p); calibrationMedian.fill(0);
  
  for(int it=0; it<niter; it++)
  {
    for(int i=0; i<n; i++)
      calibration.slice(i).row(it) =   mu.slice(it).cols(Z(it,i)*p,Z(it,i)*p+p-1).row(C(i)) - mu0.slice(it).col(Z(it,i)).t();
  }  
  
  for( int i = 0; i < n; i++)
    calibrationMedian.row(i) = median(calibration.slice(i),0);
  
  return Rcpp::List::create(  
    Rcpp::Named( "Y_cal" ) = Y - calibrationMedian,
    Rcpp::Named( "calibration_distribution" ) = calibration,
    Rcpp::Named( "calibration_median" ) = calibrationMedian 
  ) ;    
}

// [[Rcpp::export]]
SEXP compPvals2(NumericVector nullvec, NumericVector vec) { //Try above except with iterators. Runs about as fast.
  
  int n = nullvec.size();
  int m = vec.size();
  Rcpp::NumericVector pvalVec(m);
  
  typedef Rcpp::NumericVector::iterator vec_iterator;
  vec_iterator it_nv = nullvec.begin();
  vec_iterator it_v = vec.begin();
  vec_iterator it_pv = pvalVec.begin();
  
  for(int i=0; i<m; i++) {
    int count = 0;
    for(int j=0; j<n; j++) {
      count = count + (it_nv[j]>=it_v[i]);
      //cout<<count<<endl;
    }
    double val = ((double)count)/((double)n);
    //cout<<val<<endl;
    it_pv[i] = val;
  }
  
  return pvalVec;
  
}

// [[Rcpp::export]]
S4 CPP_scale_margins_sparse(S4 M, NumericVector rows, NumericVector cols, bool duplicate = true) {
  if (!M.is("dgCMatrix"))
    stop("internal error -- not a canonical sparse matrix");
  IntegerVector dims = M.slot("Dim");
  int nr = dims[0], nc = dims[1];
  if (nr != rows.size() || nc != cols.size())
    stop("internal error -- row/column weights not conformable with matrix");

  if (duplicate)
    M = clone(M);

  IntegerVector p = M.slot("p");
  IntegerVector::iterator _p = p.begin();
  IntegerVector row_of = M.slot("i");
  IntegerVector::iterator _row_of = row_of.begin();
  NumericVector x = M.slot("x");
  NumericVector::iterator _x = x.begin();
  NumericVector::iterator _rows = rows.begin();

  for (int col = 0; col < nc; col++) {
    double col_weight = cols[col];
    for (int i = _p[col]; i < _p[col+1]; i++) {
      _x[i] *= _rows[_row_of[i]] * col_weight;
    }
  }

  return M;
}

// [[Rcpp::export]]
NumericVector CPP_dsm_score_sparse(int nr, int nc, IntegerVector p, IntegerVector row_of, NumericVector f, NumericVector f1, NumericVector f2, double N, int am_code, int sparse, int transform_code) {
  if (am_code < 0 || am_code >= am_table_entries)
    stop("internal error -- invalid AM code");
  am_func AM = am_table[am_code]; /* selected association measure */

  // -- don't check whether sparse=TRUE, so power users can compute non-sparse AMs for nonzero entries of the sparse matrix
  //  if (!sparse) stop("only sparse association scores can be used with sparse matrix representation");
  
  int n_items = f.size();
  NumericVector scores(n_items);
  if (am_code != 0 && (nr != f1.size() || nc != f2.size()))
    stop("internal error -- marginal vectors f1 and f2 not conformable with matrix f");

  IntegerVector::iterator _p = p.begin();
  IntegerVector::iterator _row_of = row_of.begin();
  NumericVector::iterator _f = f.begin();
  NumericVector::iterator _f1 = f1.begin();
  NumericVector::iterator _f2 = f2.begin();
  NumericVector::iterator _scores = scores.begin();
  
  for (int col = 0; col < nc; col++) {
    for (int i = _p[col]; i < _p[col+1]; i++) {
      /* frequeny measure (*am_code == 0) is a special case, since marginals may not be available ("reweight" mode) */
      double score = (am_code == 0) ? _f[i] : AM(_f[i], _f1[_row_of[i]], _f2[col], N, sparse);
      _scores[i] = (transform_code) ? transform(score, transform_code) : score;
    }
  }

  return scores;
}

// [[Rcpp::export]]
NumericVector ptsThresh(NumericVector values,NumericVector endDates,int window){
  int valuesLen = values.length();
  NumericVector out = clone(values);
  NumericVector startDates = endDates - window;

  for (int i = 0; i < valuesLen; ++i){
    LogicalVector idx = (endDates <= endDates[i]) & (endDates >= startDates[i]);
    NumericVector valuesWindow = values[idx];
    int lenCur = valuesWindow.length();

    if (lenCur == 1){
      out[i] = (1.4 / 1.3) * valuesWindow[i];
    }
    if (lenCur == 2){
      out[i] = (1.5 / 2.4) * sum(valuesWindow);
    }
    if (lenCur == 3){
      out[i] = (1.5 / 3.3) * sum(valuesWindow);
    }
    if (lenCur == 4){
      out[i] = (1.5 / 4.0) * sum(valuesWindow);
    }
    if (lenCur >= 5){
      std::nth_element(valuesWindow.begin(),valuesWindow.begin() + 5,valuesWindow.end());
      out[i] = valuesWindow[4];
    }
  }

  return out;
}

// [[Rcpp::export]]
SEXP compPvals3(NumericVector nullvec, NumericVector vec) { //A fancyer version using iterators. Actually runs much slower!!!!
  
  int n = nullvec.size();
  int m = vec.size();
  Rcpp::NumericVector pvalVec(m);
  
  typedef Rcpp::NumericVector::iterator vec_iterator;
  //vec_iterator it_nv = nullvec.begin();
  //vec_iterator it_v = vec.begin();
  //vec_iterator it_pv = pvalVec.begin();
  
  for(vec_iterator it_v = vec.begin(), it_pv = pvalVec.begin(); it_v != vec.end(); ++it_v, ++it_pv) {
    int count = 0;
    for(vec_iterator it_nv = nullvec.begin(); it_nv != nullvec.end(); ++it_nv) {
      count = count + ( *it_nv >= *it_v );
      //cout<<count<<endl;
    }
    double val = ((double)count)/((double)n);
    //cout<<val<<endl;
    *it_pv = val;
  }
  
  return pvalVec;
  
}

// [[Rcpp::export]]
NumericVector ddirichlet(NumericVector xx , NumericVector alphaa){
  
  if (is_true(any(alphaa<=0))){
    return wrap(-1e20);
  }

  if (is_true(any( (xx <=0) ))){
    return wrap(-1e20);
  }

  double alpha_sum = std::accumulate(alphaa.begin(), alphaa.end(), 0.0); 
  
  // this will return alpha_sum = 0 0 0
  //return (alpha_sum);


  NumericVector log_gamma_alpha = lgamma(alphaa);
  NumericVector log_alpha_sum = lgamma(wrap(alpha_sum));


  double sum_log_gamma_alpha = std::accumulate(log_gamma_alpha.begin(), log_gamma_alpha.end(), 0.0); 
  
  double logD = sum_log_gamma_alpha - as<double>(log_alpha_sum);


  NumericVector a = ( alphaa - 1.0 ) * log(xx);
  return wrap( std::accumulate(a.begin(), a.end(), 0.0) - logD ); 
}

// [[Rcpp::export]]
NumericMatrix CPP_dsm_score_dense(NumericMatrix f, NumericVector f1, NumericVector f2, double N, int am_code, int sparse, int transform_code) {
  if (am_code < 0 || am_code >= am_table_entries)
    stop("internal error -- invalid AM code");
  am_func AM = am_table[am_code]; /* selected association measure */

  int nr = f.nrow(), nc = f.ncol();
  if (am_code != 0 && (nr != f1.size() || nc != f2.size()))
    stop("internal error -- marginal vectors f1 and f2 not conformable with matrix f");

  NumericMatrix scores(nr, nc);

  NumericMatrix::iterator _f = f.begin();
  NumericVector::iterator _f1 = f1.begin();
  NumericVector::iterator _f2 = f2.begin();
  NumericMatrix::iterator _scores = scores.begin();

  int i = 0;
  for (int col = 0; col < nc; col++) {
    for (int row = 0; row < nr; row++) {
      /* frequeny measure (am_code == 0) is a special case, since marginals may not be available (in "reweight" mode) */
      double score = (am_code == 0) ? _f[i] : AM(_f[i], _f1[row], _f2[col], N, sparse);
      _scores[i] = (transform_code) ? transform(score, transform_code) : score;
      i++;
    }
  }

  return scores;
}

NumericVector logLikMixHMM(NumericVector transitionMatrix, NumericVector emissionArray, 
  NumericVector initialProbs, IntegerVector obsArray, NumericMatrix coefs, 
  NumericMatrix X_, IntegerVector numberOfStates) {  
  
  
  IntegerVector eDims = emissionArray.attr("dim"); //m,p,r
  IntegerVector oDims = obsArray.attr("dim"); //k,n,r
  
  int q = coefs.nrow();
  arma::mat coef(coefs.begin(),q,coefs.ncol());
  coef.col(0).zeros();
  arma::mat X(X_.begin(),oDims[0],q);
  arma::mat lweights = exp(X*coef).t();
  if(!lweights.is_finite()){
    return wrap(-std::numeric_limits<double>::max());
    
  }
  lweights.each_row() /= sum(lweights,0);
  arma::colvec init(initialProbs.begin(),eDims[0], true);
  arma::mat transition(transitionMatrix.begin(),eDims[0],eDims[0], true);
  arma::cube emission(emissionArray.begin(), eDims[0], eDims[1], eDims[2], true);
  arma::icube obs(obsArray.begin(), oDims[0], oDims[1], oDims[2], false);
  
  arma::vec alpha(eDims[0]);
  NumericVector ll(oDims[0]);  
  double tmp;
  arma::vec initk(eDims[0]);
  
  for(int k = 0; k < oDims[0]; k++){    
    initk = init % reparma(lweights.col(k),numberOfStates);
    
    for(int i=0; i < eDims[0]; i++){      
      alpha(i) = initk(i);
      for(int r = 0; r < oDims[2]; r++){
        alpha(i) *= emission(i,obs(k,0,r),r);
      }
    }    
    
    tmp = sum(alpha);
    ll(k) = log(tmp);
    alpha /= tmp;
    
    arma::vec alphatmp(eDims[0]);
    
    for(int t = 1; t < oDims[1]; t++){  
      for(int i = 0; i < eDims[0]; i++){
        alphatmp(i) = arma::dot(transition.col(i), alpha);
        for(int r = 0; r < oDims[2]; r++){
          alphatmp(i) *= emission(i,obs(k,t,r),r);
        }
      }
      tmp = sum(alphatmp);
      ll(k) += log(tmp);
      alpha = alphatmp/tmp;
    }
  } 
  
  return ll;
}

// [[Rcpp::export]]
NumericVector  f1_gamma(NumericMatrix b,NumericVector y,NumericMatrix x,NumericVector alpha,NumericVector wt)
{
 
    // Get dimensions of x - Note: should match dimensions of
    //  y, b, alpha, and wt (may add error checking)
    
    // May want to add method for dealing with alpha and wt when 
    // constants instead of vectors
    
    int l1 = x.nrow(), l2 = x.ncol();
    int m1 = b.ncol();
    
//    int lalpha=alpha.nrow();
//    int lwt=wt.nrow();

    Rcpp::NumericMatrix b2temp(l2,1);

    arma::mat x2(x.begin(), l1, l2, false); 
    arma::mat alpha2(alpha.begin(), l1, 1, false); 
    arma::mat wt2(wt.begin(), l1, 1, false);
    
    Rcpp::NumericVector xb(l1);
    arma::colvec xb2(xb.begin(),l1,false); // Reuse memory - update both below
     
    // Moving Loop inside the function is key for speed

    NumericVector yy(l1);
    NumericVector res(m1);


    for(int i=0;i<m1;i++){
    b2temp=b(Range(0,l2-1),Range(i,i));
    arma::mat b2(b2temp.begin(), l2, 1, false); 
  

//  mu<-t(exp(alpha+x%*%b))
//  disp2<-1/wt

//  -sum(dgamma(y,shape=1/disp2,scale=mu*disp2,log=TRUE))


    xb2=exp(alpha2+ x2 * b2);
    
    for(int j=0;j<l1;j++){
      
    xb[j]=xb[j]/wt[j];  
    }

    yy=-dgamma_glmb(y,wt,xb,true);
    

    res(i) =std::accumulate(yy.begin(), yy.end(), 0.0);

    }
    
    return res;      
}

// [[Rcpp::export]]
NumericVector row_kth(NumericMatrix toSort, int k) {
  int n = toSort.rows();
  NumericVector meds = NumericVector(n);
  for (int i = 0; i < n; i++) {
    NumericVector curRow = toSort.row(i);
    std::nth_element(curRow.begin(), curRow.begin() + k, curRow.end());
    meds[i] = curRow[k];
  }

  return meds;
}

//' Empirical sample quantile.
//' Calculate empirical sample quantile.
//'
//' @param x A numeric vector, specifying the sample for which the quantile is to be calculated.
//' @param q A real number between 0 and 1 inclusive, specifying the desired quantile.
//'
//' @return The empirical quantile of the provided sample.
//' @examples
//' x<-rnorm(100)
//' qsamp(x, 0.5)
//' @export
// [[Rcpp::export]]
double qsamp(NumericVector x, double q)
{
int n = x.size();
int Q = floor(n*q);
double g = n*q-double(Q);
if (g == 0) Q -= 1; //Q = Q, but c++ uses 0-based indexing
else Q = Q; //Q -= 1 0-based indexing

std::nth_element(x.begin(), x.begin()+Q, x.end());

return x[Q];
}

// Quick and dirty implementation of lowerBound, the complement to R's
// findInterval
// [[Rcpp::export]]
IntegerVector lowerBound(const NumericVector& x, const NumericVector& breaks) {
  int n = x.size();
  IntegerVector out(n);

  for (int i = 0; i < n; i++) {
    NumericVector::const_iterator it =
      std::lower_bound(breaks.begin(), breaks.end(), x[i]);
    if (it == breaks.end()) --it;
    out[i] = it - breaks.begin() + 1;
  }
  return out;
}

// [[Rcpp::export]]
SEXP C_MedianSpec(SEXP x) {
   NumericMatrix VV(x);
   int n_specs = VV.nrow();
   int count_max = VV.ncol();
   int position = n_specs / 2; // Euclidian division
   NumericVector out(count_max);
   for (int j = 0; j < count_max; j++) { 
        NumericVector y = VV(_,j); // Copy column -- original will not be mod
        std::nth_element(y.begin(), y.begin() + position, y.end()); 
        out[j] = y[position];  
   }
   return out;
}

//' Compute complex visibility based on ABCD principle
//'
//' @param ii ABCD
//' @param kx P2VM. By default, c(0,2,2,0)
//' @param sx P2VM. By default, c(4,4,4,4)
//' @param fx Total flux in ABCD
//' 
// [[Rcpp::export]]
ComplexVector gv_abcd2vis(NumericVector ii,
  NumericVector kx = NumericVector::create(4),
  NumericVector sx = NumericVector::create(4),
  NumericVector fx = NumericVector::create(1))
{
  double vv[2] = {0,0};
  ComplexVector vis = ComplexVector(1);
  abcd2vis(ii.begin(), kx.begin(), sx.begin(), &vv[0], as<double>(fx));
  vis[0].r = vv[0];
  vis[0].i = vv[1];

  return vis;
}