将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用参数评估a arma::vec,似乎它以三种形式获取参数而没有编译错误.也就是说,如果f是a Rcpp::Function并且a是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

结果,第一种方法产生错误,第二种方法产生错误答案,只有第三种方法返回正确答案.

> # all three functions 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-05,所有四次转换都会产生正确的结果.

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()))比其他方法略快.

Answer 2

这应该适用于arma::vec,arma::rowvec和arma::colvec：

template <typename T>
Rcpp::NumericVector arma2vec(const T& x) {
    return Rcpp::NumericVector(x.begin(), x.end());
}

Answer 3

我也有同样的问题。我用wrap在几层for循环的核心进行转换，速度非常慢。我认为换行功能是拖慢速度的罪魁祸首，所以我想知道是否有一种优雅的方法来做到这一点。

至于雷蒙德的问题，您可能想尝试包含名称空间，例如：而不是在代码开头Rcpp::as<Rcpp::NumericVector>(wrap(A))包含一行。using namespace Rcpp;