Rcpp函数用于添加向量的元素

For instance,if I had paras = [1,5] and indices = [3,6]
then I would want to find the cumulative sum of the third value (3)
twice,the first value (1) and the sixth (5),to get 12. There is
additionally the option of warping the parameter values according to
their location.

因为我最清楚.

您的C代码存在一些问题.首先,不要这样做 – NumericVector warp = NA_REAL – 使用Rcpp :: Nullable<>模板(如下所示).这将解决一些问题：

>它更具可读性.如果你不熟悉Nullable类,它几乎就是它听起来的样子 – 一个可能是也可能不为null的对象.
>您不必进行任何尴尬的初始化,例如NumericVector warp = NA_REAL.坦率地说,我很惊讶编译器接受了这一点.
>您不必担心意外忘记C使用从零开始的索引,与R不同,如此行：if(NumericVector :: is_na(warp [1])){.这有不明确的行为写在它上面.

这是一个修订版本,取消了您对上述问题的引用说明：

#include <Rcpp.h>

typedef Rcpp::Nullable<Rcpp::NumericVector> nullable_t;
// [[Rcpp::export]]
double DotProd(Rcpp::NumericVector indices,Rcpp::NumericVector params,nullable_t warp_ = R_NilValue) {
  R_xlen_t i = 0,n = indices.size();
  double result = 0.0;

  if (warp_.isNull()) {
    for ( ; i < n; i++) {
      result += params[indices[i]];
    }    
  } else {
    Rcpp::NumericVector warp(warp_);
    for ( ; i < n; i++) {
      result += params[indices[i]] * warp[i];
    }  
  }

  return result;
}

您有一些精心设计的代码来生成示例数据.我没有花时间来完成这个,因为没有必要,基准测试也没有.你说自己C版本没有产生正确的结果.您的首要任务应该是让您的代码处理简单数据.然后给它提供一些更复杂的数据.然后基准.上面的修订版本适用于简单数据：

args <- list(
  indices = c(3,6),params = c(1,5),warp = c(.25,.75,1.25,1.75)
)

all.equal(
  DotProd(args[[1]],args[[2]]),dot_prod(args[[1]],args[[2]]))
#[1] TRUE

all.equal(
  DotProd(args[[1]],args[[2]],args[[3]]),args[[3]]))
#[1] TRUE

它也比此样本数据上的R版本更快.我没有理由相信它不适用于更大,更复杂的数据 – * apply函数没有什么神奇或特别的效率;它们只是更惯用/可读的R.

microbenchmark::microbenchmark(
  "Rcpp" = DotProd(args[[1]],"R" = dot_prod(args[[1]],args[[2]]))
#Unit: microseconds
#expr    min      lq     mean  median      uq    max neval
#Rcpp  2.463  2.8815  3.52907  3.3265  3.8445 18.823   100
#R    18.869 20.0285 21.60490 20.4400 21.0745 66.531   100
#
microbenchmark::microbenchmark(
  "Rcpp" = DotProd(args[[1]],args[[3]]))
#Unit: microseconds
#expr    min      lq     mean median      uq    max neval
#Rcpp  2.680  3.0430  3.84796  3.701  4.1360 12.304   100
#R    21.587 22.6855 23.79194 23.342 23.8565 68.473   100

我从上面的例子中省略了NA检查,但是也可以通过使用一点Rcpp糖将其修改成更惯用的东西.以前,你这样做：

LogicalVector indices_ok;
for (int i = 0; i < len; i++){
  indices_ok.push_back(R_IsNA(indices[i]));
}
if(is_true(any(indices_ok))){
  return NA_REAL;
}

它有点激进 – 您正在测试整个值向量(使用R_IsNA),然后应用is_true(any(indices_ok)) – 当您可能过早地中断并在R_IsNA的第一个实例(indices [i])上返回NA_REAL时结果是真的.另外,使用push_back会使你的函数变慢一些 – 你最好将indices_ok初始化为已知大小并通过循环中的索引访问来填充它.不过,这是压缩操作的一种方法：

if (Rcpp::na_omit(indices).size() != indices.size()) return NA_REAL;

为了完整起见,这里有一个完全糖化的版本,可以让你完全避免循环：

#include <Rcpp.h> 

typedef Rcpp::Nullable<Rcpp::NumericVector> nullable_t;
// [[Rcpp::export]]
double DotProd3(Rcpp::NumericVector indices,nullable_t warp_ = R_NilValue) {
  if (Rcpp::na_omit(indices).size() != indices.size()) return NA_REAL; 

  if (warp_.isNull()) {
    Rcpp::NumericVector tmp = params[indices];
    return Rcpp::sum(tmp);    
  } else {
    Rcpp::NumericVector warp(warp_),tmp = params[indices];
    return Rcpp::sum(tmp * warp); 
  }
}

/*** R

all.equal(
  DotProd3(args[[1]],args[[2]]))
#[1] TRUE

all.equal(
  DotProd3(args[[1]],args[[3]]))
#[1] TRUE

*/