9 min read

Introduction to Rcpp

In the last post I described the new memory management for vectors, which significantly improves the performance of growable vector. It is a step in the right direction in making R interpreter much faster. However, there are times when we need more speed, and R isn’t fast enough.

On the other hand, R is a great environment for joining many different tools. It might be hard or even impossible just to resign from R, especially because in most situations there are only a few critical parts of a system which must be optimized. Other stuff, like visualizations, can still be done in R and switching between different environments is not very comfortable:( And this is the point where the Rcpp comes to rescue. Rcpp is an R package which allows easily create and use C++ code within R.

Why C++?

But why should we use C++ instead of other languages? The first answer is that C++ is fast. Much faster than R. Consider this simple function from the last post. There is a simple method which copies the vector element by element (this is just dummy example, but resizing a data structure is a widely used operation):

library(microbenchmark)

fnc_assign <- function(y) {
  x <- NULL
  for (i in y) {
    x[[i]] <- i
  }
  x
}

y <- 1:1e4
microbenchmark(fnc_assign(y))
## Unit: milliseconds
##           expr      min       lq     mean   median       uq      max neval
##  fnc_assign(y) 1.935875 2.005248 2.266958 2.049815 2.149955 4.527443   100

Then we can write pretty the same function using C++ (there are some differences in the syntax, but the idea is the same).

#include <Rcpp.h>
#include <vector>

// [[Rcpp::export]]
std::vector<double> fnc_cpp(const std::vector<double>& y) {
  std::vector<double> x;
  
  for(size_t i = 0; i < y.size(); i++) {
    x.push_back(y.at(i));
  }
  return x;
}

Note that the code above is completely written in C++. There are only two non-C++ elements.

  • #include <Rcpp.h> - includes Rcpp into this file.
  • // [[Rcpp::export]] - this is the best thing in Rcpp. The function below that line is exported to the R session. Rcpp creates all the glue code under the hood.

Ok. So let see the difference in performance between C++ and R versions.

microbenchmark(fnc_assign(y), fnc_cpp(y))
## Unit: microseconds
##           expr      min        lq       mean   median       uq      max
##  fnc_assign(y) 1925.394 1954.8430 2153.30908 1980.697 2079.039 3757.623
##     fnc_cpp(y)   60.443   64.2615   99.51412   86.908   96.017 1003.241
##  neval cld
##    100   b
##    100  a
all(fnc_assign(y) == fnc_cpp(y))
## [1] TRUE

The speedup is significant! About 20 times. I think that there is no need for other explanations, why the C++ is a better choice when there’s need for speed. But what cause such big difference?

C++ is a compiled language, so before the execution, the compiler has a lot of time to look at the code and do some optimizations. It is also much closer to bare metal, and there is no intermediary between generated code and the hardware (except operating system). The R code is interpreted by the C code (R interpreter is written in C), so there is always some overhead. It is also a bit harder to do some optimizations because R is not a very strict language. For example, the user does not need to specify the type of an object because it might change during execution. In C++ the type of a variable must be declared in advance, so it is easier for the compiler to understand what is going on in that code, and how to arrange everything to get the best performance.

C++ is harder than R.

It is true. When you are writing C++ code, you need to take into account much more things like the declaring and specifying the type of variables or some memory management issues. However, in the modern C++, the differences are not as significant as it might seem.

Most of the time we still need to declare the type of a variable, but in some cases we can leave that to the compiler because it can deduce the right form from the context. See the example.

#include <Rcpp.h>
#include <vector>

// [[Rcpp::plugins(cpp11)]]
// [[Rcpp::export]]
int auto_type(const std::vector<double>& y) {

  // We assign y to x, so it is obvious that x must have the same class as y.
  auto x = y; 

  // Size method returns exact type, so the variable for storing the result can inherit that class.
  auto size = x.size() / 2;
  
  return size;
}

We also needed to include // [[Rcpp::plugins(cpp11)]] statement in the code above. It tells the Rcpp that modern version of C++ is used, so it must inform the C++ compiler about this fact.

Of course, this functions is exported to R and works pretty well;)

auto_type(rnorm(100))
## [1] 50

It is also worth to note that the C++11 has a new version of for loop. It operates in a very similar fashion as the R’s one.

#include <Rcpp.h>
#include <vector>

// [[Rcpp::plugins(cpp11)]]
// [[Rcpp::export]]
std::vector<double> fnc_cpp2(const std::vector<double>& y) {
  std::vector<double> x;
  
  for(const auto& elem : y) {
    x.push_back(elem);
  }
  
  return x;
}

There is no real change in performance between those versions of for loop.

microbenchmark(fnc_cpp(y), fnc_cpp2(y))
## Unit: microseconds
##         expr    min     lq     mean  median      uq      max neval cld
##   fnc_cpp(y) 59.905 68.523 107.1801 92.1965 93.8315 2207.083   100   a
##  fnc_cpp2(y) 59.943 89.475 118.0960 91.5580 94.0170 2158.232   100   a

R object in C++.

Rcpp is an incredible tool. Not only it magically generates all the code needed to include C++ functions inside R but also allows to use R objects directly inside C++. See the example:

#include <Rcpp.h>
#include <vector>

// [[Rcpp::plugins(cpp11)]]
// [[Rcpp::export]]
Rcpp::NumericVector fnc_cpp_rstruct(Rcpp::NumericVector y) {
  Rcpp::NumericVector x;
  
  for(const auto& elem : y) {
    x.push_back(elem);
  }
  
  return x;
}

Now the x and y are the standard R objects and behave like that, with all pros and cons.

There is a problem with that code above. Rcpp (version. 0.12.10) still uses the old memory allocation scheme, so each time the new element is inserted the whole vector is reallocated (see my last post about this issue http://www.blog.zstat.pl/2017/05/06/updated-vectors-in-r-3.4.0/). It can be easily shown using simple benchmarks, but it will possibly be changed in the next versions.

# It is a function from the last post which allocates memory in the wrong way.
fnc_combines <- function(y) {
  x <- NULL
  for (i in y) {
    x <- c(x, i)
  }
  x
}

microbenchmark(fnc_cpp(y), fnc_cpp_rstruct(y), fnc_combines(y))
## Unit: microseconds
##                expr        min         lq        mean     median
##          fnc_cpp(y)     61.353     77.339    141.5711    103.511
##  fnc_cpp_rstruct(y) 113530.535 115608.046 148654.4390 116595.240
##     fnc_combines(y) 117104.505 118327.301 129903.8500 118973.359
##          uq        max neval cld
##     109.429   1421.698   100 a  
##  221778.314 233912.565   100   c
##  120151.078 228821.605   100  b

Other R objects.

Below I show an example of how to create a named list on the C++ side. It is very useful because you can easily export export nearly everything as a list:) Particularly useful for debugging.

#include <Rcpp.h>
#include <vector>

// [[Rcpp::plugins(cpp11)]]
// [[Rcpp::export]]
Rcpp::List create_test_list() {
  
  
  Rcpp::NumericVector x = {1.0, 2.0};
  Rcpp::NumericVector y = {1.0, 2.0, 3.0};
  Rcpp::NumericVector z = {1.0, 2.0, 3.0, 4.0, 5.0};
  
  return Rcpp::List::create(Rcpp::Named("x") = x,
                          Rcpp::Named("y") = y,
                          Rcpp::Named("z") = z);
  
}
create_test_list()
## $x
## [1] 1 2
## 
## $y
## [1] 1 2 3
## 
## $z
## [1] 1 2 3 4 5

A word about memory management.

You might have heard somewhere that C++ is very compilated when it comes to memory management. There are some things like pointers, which can cause memory leaks and other scary stuff. But as long as you use standard data structures like std::vector or R’s ones like Rcpp::NumericVector, there’s no need to worry about anything. C++ does an excellent job of handling this material. So I would recommend staying away from pointers and other more compilated stuff, at least at the beginning.

I’m planning to write a post covering some basic memory management issues in C++, so stay focus;)

A simple example.

Below there is a simple example of using std::map to create a bit faster table function for a one-dimensional vector. You can find more about std::map here: http://www.cplusplus.com/reference/map/map/.

#include <Rcpp.h>
#include <vector>
#include <map>

// [[Rcpp::plugins(cpp11)]]
// [[Rcpp::export]]
Rcpp::IntegerVector fast_table(Rcpp::NumericVector x) {
  // Create map with double value as a key and size_t to  store number of occurence.
  std::map<double, size_t> tbl; // 
  
  // How many unique values does the vector have
  size_t unique = 0;
  
  for(const auto& it : x) {
    // The result of tbl.find is an iterator, but this is not important -
    // the type is deducted using auto.
    auto mp = tbl.find(it);

    // When there is no such key in the map, a function
    // returns value equal to end() iterator.
    // In such situation that new key must be added.
    if(mp == tbl.end()) {
      tbl[it] = 1;
      unique++;
      
    } else {

      // When such key exists in the map, the number of occurrences is updated.
      // Note:
      // mp->first - key.
      // mp->second - value.
      mp->second++;
    }
  }
  
  // rewriting the map content into NumericVector
  Rcpp::IntegerVector res(unique);
  Rcpp::NumericVector names(unique);
  auto it = tbl.begin();
  for(size_t i = 0; i < unique; i++) {
    res.at(i) = it->second;
    names.at(i) = it->first;
    it++;
  }
  
  res.names() = names;
  return res;
}

And here are the results - a bit faster:

set.seed(123)

x <- sample(1:1000, size = 1e5, replace = TRUE)

microbenchmark(table(x), fast_table(x))
## Unit: milliseconds
##           expr       min        lq      mean    median        uq       max
##       table(x) 37.609261 38.095833 41.494076 39.453395 39.998501 151.51534
##  fast_table(x)  9.401862  9.632002  9.874784  9.784468  9.910705  12.64195
##  neval cld
##    100   b
##    100  a
all(table(x) == fast_table(x))
## [1] TRUE

Summary.

C++ is not a mysterious language used only by experienced developers with proper programming training. It can be easily embedded into R code to speed up the critical paths using Rcpp package.

In the modern version of C++, you can use auto type to declare variable type. The compiler should deduce the proper value by itself.

You can use all C++ data structures in your code. There is a short list of the standard possibilities -http://www.cplusplus.com/reference/stl/. You can learn how to use them by looking at the examples, e.g. http://www.cplusplus.com/reference/map/map/begin/.

I’m working on new stuff covering a bit more advanced (but VERY important) aspects of Rcpp, so don’t forget to visit my blog from time to time. You can also follow me on twitter https://twitter.com/zzawadz or leave a comment below.

comments powered by Disqus