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

Question

我有一个很长的参数向量（大约 4^10 个元素）和一个索引向量。我的目标是将索引向量中索引的所有参数值加在一起。

例如，如果我有 paras = [1,2,3,4,5,5,5] 和 indices = [3,3,1,6] 那么我想找到的累计总和第三个值 (3) 两次，第一个值 (1) 和第六个 (5)，得到 12。另外还有根据位置扭曲参数值的选项。

我正在尝试加快 R 实现的速度，因为我已调用它数百万次。

我现在的代码总是returnsNA，我看不出哪里错了

这是 Rcpp 函数：

double dot_prod_c(NumericVector indices, NumericVector paras, 
                   NumericVector warp = NA_REAL) {
int len = indices.size();
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;
}
double counter = 0;
if(NumericVector::is_na(warp[1])){
    for (int i = 0; i < len; i++){
        counter += paras[indices[i]];
    }
} else {
    for (int i = 0; i < len; i++){
        counter += paras[indices[i]] * warp[i];
    }
}
return counter;
}

这是可用的 R 版本：

dot_prod <- function(indices, paras, warp = NA){
    if(is.na(warp[1])){
        return(sum(sapply(indices, function(ind) paras[ind + 1])))
    } else {
        return(sum(sapply(1:length(indices), function(i){
            ind <- indices[i]
            paras[ind + 1] * warp[i]
        })))
    }
}

下面是一些使用 microbenchmark 包进行测试和基准测试的代码：

# testing
library(Rcpp)
library(microbenchmark)

parameters <- list()
indices <- list()
indices_trad <- list()

set.seed(2)
for (i in 4:12){
    size <- 4^i
    window_size <- 100
    parameters[[i-3]] <- runif(size)
    indices[[i-3]] <- floor(runif(window_size)*size)
    temp <- rep(0, size)
    for (j in 1:window_size){
        temp[indices[[i-3]][j] + 1] <- temp[indices[[i-3]][j] + 1] + 1
    }
    indices_trad[[i-3]] <- temp
}

microbenchmark(
    x <- sapply(1:9, function(i) dot_prod(indices[[i]], parameters[[i]])),
    x_c <- sapply(1:9, function(i) dot_prod_c(indices[[i]], parameters[[i]])),
    x_base <- sapply(1:9, function(i) indices_trad[[i]] %*% parameters[[i]])
)
all.equal(x, x_base) # is true, does work
all.equal(x_c, x_base) # not true - C++ version returns only NAs

Answer 1

我在尝试通过您的代码来解释您的总体目标时遇到了一些麻烦，所以我将继续这个解释

For instance, if I had paras = [1,2,3,4,5,5,5] and indices = [3,3,1,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 class，它几乎就是它听起来的样子 - 一个可能为 null 也可能不为 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, 3, 1, 6),
  params = c(1, 2, 3, 4, 5, 5, 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]]), 
  dot_prod(args[[1]], args[[2]], args[[3]]))
#[1] TRUE

在此示例数据上，它也比 R 版本更快。我没有理由相信它也不会用于更大、更复杂的数据——*apply 函数没有什么神奇或特别有效的；它们只是更惯用/可读的 R。

microbenchmark::microbenchmark(
  "Rcpp" = DotProd(args[[1]], args[[2]]), 
  "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[[2]], args[[3]]), 
  "R" = dot_prod(args[[1]], args[[2]], 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)) - 当您可能过早地中断并且 return NA_REAL 在 R_IsNA(indices[i]) 的第一个实例上导致 true。此外，使用 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, Rcpp::NumericVector params, 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]]), 
  dot_prod(args[[1]], args[[2]]))
#[1] TRUE

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

*/

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

Rcpp function for adding elements of a vector

r

vector

rcpp