SC2

Statistical Computing 2

3. Rcpp basics

Rcpp via R CMD SHLIB: the old but instructive way

To illustrate how Rcpp works, let us consider again our C function for exponential smoothing, that is:

#include <R.h>
#include <Rinternals.h>

SEXP expSmooth2(SEXP y, SEXP a)
{
  int ni;
  double *xy, *xys;
  double ai;
  SEXP ys;
  
  y = PROTECT(coerceVector(y, REALSXP));
  ni = length(y);
  ys = PROTECT(allocVector(REALSXP, ni));
  
  ai = REAL(a)[0];
  xy = REAL(y); 
  xys = REAL(ys);
  
  xys[0] = xy[0];
  for(int i = 1; i < ni; i++){
    xys[i] = ai * xys[i-1] + (1 - ai) * xy[i];
  }
  
  UNPROTECT(2);
  
  return ys;
}

Now, the same function can be implemented in “raw” Rcpp as follows (as we will explain later, this is NOT the recommended way of using Rcpp):

#include <Rcpp.h>
using namespace Rcpp;

RcppExport SEXP expSmoothRcpp_manual(SEXP ySEXP, SEXP aSEXP)
{
  const NumericVector y = as<const NumericVector>(ySEXP);
  const double a = as<const double>(aSEXP);
  
  int ni = y.size();
  NumericVector ys(ni);
  
  ys[0] = y[0];
  for(int i = 1; i < ni; i++){
    ys[i] = a * ys[i-1] + (1 - a) * y[i];
  }
  
  return ys;
}

Don’t worry if you don’t understand what is going on in the code above, we’ll explain it below. As you can see the Rcpp function for exponential smoothing is much more compact. On the Linux command line it can be compiled by doing:

library(Rcpp)
system("export PKG_CXXFLAGS=`Rscript -e \"Rcpp:::CxxFlags()\"` && 
        R CMD SHLIB expSmoothRcpp_manual.cpp")

Note that the call to R CMD SHLIB is preceded by some extra code, which is needed to let SHLIB know where Rcpp.h can be found. We then load the compiled code as usual:

dyn.load("expSmoothRcpp_manual.so")

Let’s see whether the new function actually works:

system("R CMD SHLIB expSmooth2.c")
dyn.load("expSmooth2.so")

y <- rnorm(100)
max( abs(.Call("expSmoothRcpp_manual", y, 0.9) - .Call("expSmooth2", y, 0.9)) )
## [1] 0

So the C and Rcpp functions produce exactly the same output.

Now, let us try to understand how the Rcpp function for exponential smoothing works:

It is worth repeating that compiling Rcpp code via R CMD SHLIB is not recommended, and the next section will illustrate some better methods. The main advantage of Rcpp illustrated so far is that we didn’t explicitly create any object of class SEXP, which means that we didn’t have to call PROTECT and UNPROTECT to do memory management. In fact, this was done automatically for us via the NumericVector class, which wraps an R object of class SEXP and protects it from R’s garbage collector while it is in scope (that is, until we exit the part of the program where the variable is accessible). Hence, our code is shorter and one source of mistakes (i.e., incorrect memory management) has been removed. Rcpp provides many other wrappers around R objects, such as:

The pattern is always the same, upon passing from R to a C++ function, wrap all the SEXP inputs via one of the wrapper classes defined in Rcpp, so that you will be able manipulate them in C++ without thinking about R memory management. As we will see in the next section, the conversion from SEXP to an Rcpp wrapper can be handled automatically by Rcpp.

Rcpp via sourceCpp()

Consider the following C++ function:

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export(name = "expSmoothRcpp")]]
NumericVector expSmoothRcpp_I(const NumericVector y, const double a)
{
  int ni = y.size();
  NumericVector ys(ni);
  
  ys[0] = y[0];
  for(int i = 1; i < ni; i++){
    ys[i] = a * ys[i-1] + (1 - a) * y[i];
  }
  
  return ys;
}

This is very similar to expSmoothRcpp_manual (see above). In fact, if we leave aside the comment:

// [[Rcpp::export(name = "expSmoothRcpp")]]

for a moment, the main difference is that now the function signature has changed from:

RcppExport SEXP expSmoothRcpp_manual(SEXP ySEXP, SEXP aSEXP)

to:

NumericVector expSmoothRcpp_I(const NumericVector y, const double a)

Hence, the function is not using RcppExport and its inputs and output are not of class SEXP, which means that this function cannot be directly called from R using .Call(). To make this function accessible from R, we can do:

sourceCpp("expSmoothRcpp_I.cpp")

which compiles expSmoothRcpp_I, loads the corresponding dynamic library in R (using dyn.load) and creates an R wrapper called expSmoothRcpp. Hence, to test if this worked we do:

max( abs(expSmoothRcpp(y, 0.9) - .Call("expSmooth2", y, 0.9)) )
## [1] 0

Perfect agreement.

Now, to understand how this works, note that:

expSmoothRcpp
## function (y, a) 
## .Call(<pointer: 0x7f80bfb018f0>, y, a)

shows that expSmoothRcpp is a simple wrapper around a .Call to a C++ function called <pointer: 0x7ff6985a0520> (or similar). The R function and the C++ function called in the internal call to .Call have been generated by the call to sourceCpp. We can get more details about what sourceCpp is doing, by using the verbose argument:

sourceCpp("expSmoothRcpp_I.cpp", rebuild = TRUE, verbose = TRUE)
## 
## Generated extern "C" functions 
## --------------------------------------------------------
## 
## 
## #include <Rcpp.h>
## #ifdef RCPP_USE_GLOBAL_ROSTREAM
## Rcpp::Rostream<true>&  Rcpp::Rcout = Rcpp::Rcpp_cout_get();
## Rcpp::Rostream<false>& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get();
## #endif
## 
## // expSmoothRcpp_I
## NumericVector expSmoothRcpp_I(const NumericVector y, const double a);
## RcppExport SEXP sourceCpp_1_expSmoothRcpp_I(SEXP ySEXP, SEXP aSEXP) {
## BEGIN_RCPP
##     Rcpp::RObject rcpp_result_gen;
##     Rcpp::RNGScope rcpp_rngScope_gen;
##     Rcpp::traits::input_parameter< const NumericVector >::type y(ySEXP);
##     Rcpp::traits::input_parameter< const double >::type a(aSEXP);
##     rcpp_result_gen = Rcpp::wrap(expSmoothRcpp_I(y, a));
##     return rcpp_result_gen;
## END_RCPP
## }
## 
## Generated R functions 
## -------------------------------------------------------
## 
## `.sourceCpp_1_DLLInfo` <- dyn.load('/tmp/RtmpNhjugc/sourceCpp-x86_64-pc-linux-gnu-1.0.10/sourcecpp_c42f4bb5317e/sourceCpp_3.so')
## 
## expSmoothRcpp <- Rcpp:::sourceCppFunction(function(y, a) {}, FALSE, `.sourceCpp_1_DLLInfo`, 'sourceCpp_1_expSmoothRcpp_I')
## 
## rm(`.sourceCpp_1_DLLInfo`)
## 
## Building shared library
## --------------------------------------------------------
## 
## DIR: /tmp/RtmpNhjugc/sourceCpp-x86_64-pc-linux-gnu-1.0.10/sourcecpp_c42f4bb5317e
## 
## /opt/R/4.2.2/lib/R/bin/R CMD SHLIB --preclean -o 'sourceCpp_3.so' 'expSmoothRcpp_I.cpp'

The first part of the text output states that some extern "C" functions have been generated. In particular, we see the definition of the function:

RcppExport SEXP sourceCpp_1_expSmoothRcpp_I(SEXP ySEXP, SEXP aSEXP)

The fact that the definition of sourceCpp_1_expSmoothRcpp is preceded by RcppExport and that the inputs and output are of class SEXP means that the function can be called from R using .Call (remember that this is a requirement of .Call). The fact that sourceCpp_1_expSmoothRcpp is compatible with .Call and that it contains a call to expSmoothRcpp_I (our C++ code, above) shows that sourceCpp_1_expSmoothRcpp is a wrapper around expSmoothRcpp_I, making it accessible from R. We now examine the function body:

So, to make our expSmoothRcpp_I function accessible from R, sourceCpp has written a C++ wrapper which is compatible with .Call and does all the necessary conversions from SEXP to NumericVector, double and potentially many other Rcpp or C++ classes for us. The next lines in the text output show that sourceCpp has:

The last thing that needs to be explained is the purpose of the comment preceding our C++ function definition, that is:

// [[Rcpp::export(name = "expSmoothRcpp")]]
NumericVector expSmoothRcpp_I(const NumericVector y, const double a)

Here Rcpp::exportRcpp is an example of the Rcpp export attribute, used to declare that the expSmoothRcpp_I C++ function should be callable from R via the automatically generated R and C++ code we have described above. The optional (name = "expSmoothRcpp") argument is used to specify the name of the R function that should be used to access the C++ code. In fact, simply doing:

// [[Rcpp::export]]
NumericVector expSmoothRcpp_I(const NumericVector y, const double a)

would have created an R wrapper called expSmoothRcpp_I. For more on attributes, see the Rcpp vignette on attributes, which also clarifies what are the requirements that must met by a C++ function to be exportable to R via Rcpp::exportRcpp.

Inline Rcpp via cppFunction() and evalCpp()

While it is generally preferable to have C++ code stored in its own .cpp source files, Rcpp provides mechanisms for defining and executing C++ functions within an R script. In particular, consider the following standard Rcpp example:

sourceCpp(code = '
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export(name = "fiboInline1")]]
int fibonacci(const int n) {
 if (n < 2)
  return n;
 else
  return (fibonacci(n-1)) + fibonacci(n-2);
}')

which creates a C++ function for computing the \(n\)-th element of the Fibonacci sequence:

fiboInline1(6)
## [1] 8

This is equivalent to having defined the function in, say, fibo.cpp and having sourced it using sourceCpp("fibo.cpp"). An even more compact approach involves using the cppFunction function:

cppFunction('
int fiboInline2(const int x) {
 if (x < 2)
  return x;
 else
  return (fiboInline2(x-1)) + fiboInline2(x-2);
}
')

fiboInline2(6)
## [1] 8

Using Rcpp in an inline fashion can be convenient we dealing with very small chunks of C++ code. Simple C++ expression can be evaluated via evalCpp, for instance:

evalCpp('1.0 + 1.0')
## [1] 2

As for sourceCpp, we can use the verbose argument to look at the code automatically generated by Rcpp:

evalCpp('1.0 + 1.0', verbose = TRUE)
## 
## Generated code for function definition: 
## --------------------------------------------------------
## 
## #include <Rcpp.h>
## 
## using namespace Rcpp;
## 
## // [[Rcpp::export]]
## SEXP get_value(){ return wrap( 1.0 + 1.0 ) ; }
## 
## No rebuild required (use rebuild = TRUE to force a rebuild)
## [1] 2

Notice the use of Rcpp::wrap when the get_value function returns. wrap allows to transfer objects from C++ to R, for example here it converts the double 2.0 to a numeric vector in R. An important point is that, when converting from R to C++ using Rcpp::as<T>, we need to specify the C++ type via the T argument, while the conversion from C++ to R using wrap is handled automatically.

References