Search code examples
c++rrcpp

Rcpp: one function works on NumericVector, the other needs to treat it as scalar


I have two functions written with Rcpp.

The first function f takes a NumericVector as input and returns a NumericVector as output. But often, f is meant to be treated as a scalar function. When I call it from R, it doesn't matter, since R will treat any scalar as a vector anyways, so I don't have to bother about it. I can both use f as a scalar or a vector function, e.g. f(1) and f(c(1,2,3,4) both work.

  1. Question 1: Is this a bad idea? Should I make TWO separate functions, one scalar, and one vectorized? For example, is there some unwanted overhead when I call f from R using scalar input rather than a vector?

The second function g needs to use f purely as a scalar function. This causes a problem because I need to call f from Rcpp when I am defining g.... now I can't use it as a scalar function by calling f(x) for some double x.

I tried doing something like this from inside the definition of g:

NumericVector x = NumericVector::create(value); //create vector of length 1
f(x) //use f(x) as scalar 

but timing the code took absurdly long (since the function g needs to perform the above operation many times).

  1. Question 2: Is there a way to more quickly treat f as a scalar function just like I would if I was in R?

My current solution is to do what I asked in Queston 1: create another f function which is scalar, and then call THAT one from g, so that inside g we simply say

alternative_f(x) //use alternative f as scalar

Solution

  • You could factor out the scalar code into a separate (inline) function and call that from both f (with std::transform) and g (+ additional logic). Simplified example:

    #include <Rcpp.h>
    
    inline double f_impl (double x) {
        return x * 2;
    }
    
    
    // [[Rcpp::export]]
    double g (double x) {
        return f_impl(x) + 2;
    }
    
    // [[Rcpp::export]]
    Rcpp::NumericVector f (Rcpp::NumericVector x) {
        Rcpp::NumericVector y(Rcpp::no_init(x.length()));
        std::transform(x.begin(), x.end(), y.begin(), f_impl);
        return y;
    }