Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Convert RcppArmadillo vector to Rcpp vector

Tags:

r

rcpp

armadillo

I am trying to convert RcppArmadillo vector (e.g. arma::colvec) to a Rcpp vector (NumericVector). I know I can first convert arma::colvec to SEXP and then convert SEXP to NumericVector (e.g. as<NumericVector>(wrap(temp)), assuming temp is an arma::colvec object). But what is a good way to do that?

I want to do that simply because I am unsure if it is okay to pass arma::colvec object as a parameter to an Rcpp::Function object.

like image 569
Raymond Wong Avatar asked Jan 10 '13 07:01

Raymond Wong


3 Answers

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.

like image 155
Randy Lai Avatar answered Nov 07 '22 23:11

Randy Lai


This should works with arma::vec, arma::rowvec and arma::colvec:

template <typename T>
Rcpp::NumericVector arma2vec(const T& x) {
    return Rcpp::NumericVector(x.begin(), x.end());
}
like image 43
Artem Klevtsov Avatar answered Nov 07 '22 23:11

Artem Klevtsov


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.

like image 37
Clavichord Avatar answered Nov 08 '22 00:11

Clavichord