14

我正在尝试将 RcppArmadillo 向量(例如arma::colvec)转换为 Rcpp 向量(NumericVector)。我知道我可以先转换arma::colvecSEXP,然后再转换SEXPNumericVector(例如as<NumericVector>(wrap(temp)),假设 temp 是一个arma::colvec对象)。但是有什么好的方法呢?

我想这样做只是因为我不确定是否可以将arma::colvec对象作为参数传递给Rcpp::Function对象。

4

3 回答 3

9

I was trying to Evaluate a Rcpp::Function with argument arma::vec, it seems that it takes the argument in four forms without compilation errors. That is, if f is a Rcpp::Function and a is a arma::vec, then

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

produce no compilation and runtime errors, at least apparently.

For this reason, I have conducted a little test for the four versions of arguments. Since I suspect that somethings will go wrong in garbage collection, I test them again 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

As a result, the first method produces an error, the second method produces a wrong answer and only the third and the fourth method return the 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

Note that, if gctorture is set FALSE, then all functions will return a correct result.

> 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

It means that method 1 and method 2 are subjected to break when garbage is collected during runtime and we don't know when it happens. Thus, it is dangerous to not wrap the parameter properly.

Edit: as of 2017-12-05, all four conversions produce the correct result.

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

and this is the benchmark

> 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

And f(NumericVector(a.begin(),a.end())) is marginally faster than other methods.

于 2014-03-12T07:08:26.340 回答
6

这应该与arma::vec,arma::rowvec和一起使用arma::colvec

template <typename T>
Rcpp::NumericVector arma2vec(const T& x) {
    return Rcpp::NumericVector(x.begin(), x.end());
}
于 2016-05-24T14:25:25.413 回答
1

I had the same question. I used wrap to do the conversion at the core of several layers of for loops and it was very slow. I think the wrap function is to blame for dragging the speed down so I wish to know if there is an elegant way to do this.

As for Raymond's question, you might want to try including the namespace like: Rcpp::as<Rcpp::NumericVector>(wrap(A)) instead or include a line using namespace Rcpp; at the beginning of your code.

于 2013-01-12T19:47:18.203 回答