r - 将 RcppArmadillo 向量转换为 Rcpp 向量

Question

我正在尝试将 RcppArmadillo 向量(例如 arma::colvec )转换为 Rcpp 向量( NumericVector )。我知道我可以先转换arma::colvec至SEXP然后转换SEXP至NumericVector (例如 as<NumericVector>(wrap(temp) )，假设 temp 是 arma::colvec目的)。但是什么是做到这一点的好方法呢？

我想这样做只是因为我不确定是否可以通过 arma::colvec对象作为 Rcpp::Function 的参数对象。

Answer 1

我试图用参数 Rcpp::Function 评估 arma::vec ，它似乎采用四种形式的参数，没有编译错误。也就是说，如果 f 是 Rcpp::Function 并且 a 是 arma::vec ，那么

1. f(a)
2. f(wrap(a))
3. f(as<NumericVector>(wrap(a)))
4. f(NumericVector(a.begin(),a.end()))

至少表面上不会产生编译和运行时错误。

为此，我对四个版本的论证进行了一点测试。由于我怀疑垃圾收集中会出现问题，因此我再次测试它们 gctorture 。

gctorture(on=FALSE)
Rcpp::sourceCpp(code = '
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

using namespace Rcpp;

// [[Rcpp::export]]
double foo1(arma::vec a, arma::vec b, Function f){
    double sum = 0.0;
    for(int i=0;i<100;i++){
        sum += as<double>(f(a, b));
    }
    return sum;
}

// [[Rcpp::export]]
double foo2(arma::vec a, arma::vec b, Function f){
    double sum = 0.0;
    for(int i=0;i<100;i++){
        sum += as<double>(f(wrap(a),wrap(b)));
    }
    return sum;
}

// [[Rcpp::export]]
double foo3(arma::vec a, arma::vec b, Function f){
    double sum = 0.0;
    for(int i=0;i<100;i++){
        sum += as<double>(f(as<NumericVector>(wrap(a)),as<NumericVector>(wrap(b))));
    }
    return sum;
}

// [[Rcpp::export]]
double foo4(arma::vec a, arma::vec b, Function f){
    double sum = 0.0;
    for(int i=0;i<100;i++){
        sum += as<double>(f(NumericVector(a.begin(),a.end()),NumericVector(b.begin(),b.end())));
    }
    return sum;
}
')
# note that when gctorture is on, the program will be very slow as it
# tries to perfrom GC for every allocation.
# gctorture(on=TRUE)
f = function(x,y) {
    mean(x) + mean(y)
}
# all three functions should return 700
foo1(c(1,2,3), c(4,5,6), f) # error
foo2(c(1,2,3), c(4,5,6), f) # wrong answer (occasionally)!
foo3(c(1,2,3), c(4,5,6), f) # correct answer
foo4(c(1,2,3), c(4,5,6), f) # correct answer

因此，第一个方法会产生错误，第二个方法会产生错误答案，只有第三个和第四个方法会返回正确答案。

> # they should return 700
> foo1(c(1,2,3), c(4,5,6), f) # error
Error: invalid multibyte string at '<80><a1><e2>'
> foo2(c(1,2,3), c(4,5,6), f) # wrong answer (occasionally)!
[1] 712
> foo3(c(1,2,3), c(4,5,6), f) # correct answer
[1] 700
> foo4(c(1,2,3), c(4,5,6), f) # correct answer
[1] 700

请注意，如果 gctorture 设置为 FALSE ，则所有函数都会返回正确的结果。

> foo1(c(1,2,3), c(4,5,6), f) # error
[1] 700
> foo2(c(1,2,3), c(4,5,6), f) # wrong answer (occasionally)!
[1] 700
> foo3(c(1,2,3), c(4,5,6), f) # correct answer
[1] 700
> foo4(c(1,2,3), c(4,5,6), f) # correct answer
[1] 700

这意味着方法1和方法2在运行时垃圾收集时会受到破坏，但我们不知道它何时发生。因此，不正确包装参数是危险的。

编辑:截至 2017 年 12 月 5 日，所有四次转换均产生正确的结果。

1. f(a)
2. f(wrap(a))
3. f(as<NumericVector>(wrap(a)))
4. f(NumericVector(a.begin(),a.end()))

这是基准

> microbenchmark(foo1(c(1,2,3), c(4,5,6), f), foo2(c(1,2,3), c(4,5,6), f), foo
3(c(1,2,3), c(4,5,6), f), foo4(c(1,2,3), c(4,5,6), f))
Unit: milliseconds
                            expr      min       lq     mean   median       uq
 foo1(c(1, 2, 3), c(4, 5, 6), f) 2.575459 2.694297 2.905398 2.734009 2.921552
 foo2(c(1, 2, 3), c(4, 5, 6), f) 2.574565 2.677380 2.880511 2.731615 2.847573
 foo3(c(1, 2, 3), c(4, 5, 6), f) 2.582574 2.701779 2.862598 2.753256 2.875745
 foo4(c(1, 2, 3), c(4, 5, 6), f) 2.378309 2.469361 2.675188 2.538140 2.695720
      max neval
 4.186352   100
 5.336418   100
 4.611379   100
 3.734019   100

f(NumericVector(a.begin(),a.end())) 比其他方法稍微快一些。