In this chapter we’ll learn how to rewrite R code in C++ to make it faster using the Rcpp package. The Rcpp package makes it simple to connect C++ to R! With C++ you can fix:
Loops that can’t be easily vectorised because subsequent iterations depend on previous ones.
Recursive functions, or problems which involve calling functions millions of times. The overhead of calling a function in C++ is much lower than in R.
Problems that require advanced data structures and algorithms that R doesn’t provide. Through the standard template library (STL), C++ has efficient implementations of many important data structures, from ordered maps to double-ended queue
Like how?
Getting started with C++
library(Rcpp)
Install a C++ compiler:
Rtools, on Windows
Xcode, on Mac
Sudo apt-get install r-base-dev or similar, on Linux.
First example
Rcpp compiling the C++ code:
cppFunction('int add(int x, int y, int z) { int sum = x + y + z; return sum;}')# add works like a regular R functionadd
#> function (x, y, z)
#> .Call(<pointer: 0x00007ff9d7e315f0>, x, y, z)
add(1, 2, 3)
#> [1] 6
Some things to note:
The syntax to create a function is different.
Types of inputs and outputs must be explicitly declared
Use = for assignment, not <-.
Every statement is terminated by a ;
C++ has it’s own name for the types we are used to:
scalar types are int, double, bool and String
vector types (for Rcpp) are IntegerVector, NumericVector, LogicalVector and CharacterVector
Other R types are available in C++: List, Function, DataFrame, and more.
Explicitly use a return statement to return a value from a function.
Source stand-alone C++ files into R using sourceCpp()
C++ files have extension .cpp
#include <Rcpp.h>
using namespace Rcpp;
And for each function that you want available within R, you need to prefix it with:
// [[Rcpp::export]]
Inside a cpp file you can include R code using special comments
/*** R
rcode here
*/
Example
This block in Rmarkdown uses {Rcpp} as a short hand for engine = “Rcpp”.
#include <Rcpp.h>usingnamespace Rcpp;// [[Rcpp::export]]double meanC(NumericVector x){int n = x.size();double total =0;for(int i =0; i < n;++i){ total += x[i];}return total / n;}/*** Rx <- runif(1e5)bench::mark( mean(x), meanC(x))*/
NOTE: For some reason although the r code above runs, knit doesn’t include the output. Why?
#> List of 4
#> $ : num NA
#> $ : int NA
#> $ : logi NA
#> $ : chr NA
Standard Template Library
STL provides powerful data structures and algorithms for C++.
Iterators
Iterators are used extensively in the STL to abstract away details of underlying data structures.
If you an iterator it, you can:
Get the value by ‘dereferencing’ with *it
Advance to the next value with ++it
Compare iterators (locations) with ==
Algorithms
The real power of iterators comes from using them with STL algorithms.
A good reference is [https://en.cppreference.com/w/cpp/algorithm]
Book provides examples using accumulate and upper_buond
Another Example:
#include <algorithm>#include <Rcpp.h>usingnamespace Rcpp;// Explicit iterator version// [[Rcpp::export]]NumericVector square_C_it(NumericVector x){ NumericVector out(x.size());// Each container has its own iterator type NumericVector::iterator in_it; NumericVector::iterator out_it;for(in_it = x.begin(), out_it = out.begin(); in_it != x.end();++in_it,++out_it){*out_it = pow(*in_it,2);}return out;}// Use algorithm 'transform'// [[Rcpp::export]]NumericVector square_C(NumericVector x){ NumericVector out(x.size());std::transform(x.begin(),x.end(), out.begin(),[](double v)->double{return v*v;});return out;}
square_C(c(1.0,2.0,3.0))
#> [1] 1 4 9
square_C_it(c(1.0,2.0,3.0))
#> [1] 1 4 9
Data Structures
STL provides a large set of data structures. Some of the most important:
std::vector - like an R vector, except knows how to grow efficiently
std::unordered_set - unique set of values. Ordered version std::set. Unordered is more efficient.
std::map - Moslty similar to R lists, provide an association between a key and a value. There is also an unordered version.
A quick example illustrating the map:
#include <Rcpp.h>usingnamespace Rcpp;// [[Rcpp::export]]std::map<double,int> tableC(NumericVector x){// Note the types are <key, value>std::map<double,int> counts;int n = x.size();for(int i =0; i < n; i++){ counts[x[i]]++;}return counts;}
res =tableC(c(1,1,2,1,4,5))res
#> 1 2 4 5
#> 3 1 1 1
Note that the map is converted to a named vector in this case on return
To learn more about the STL data structures see containers at cppreference
Case Studies
Case Study
Real life uses of C++ to replace slow R code.
Case study 1: Gibbs sampler
The Gibbs sampler is a method for estimating parameters expectations. It is a MCMC algorithm that has been adapted to sample from multidimensional target distributions. Gibbs sampling generates a Markov chain of samples, each of which is correlated with nearby samples.
“Darren Wilkinson stresses the rather pragmatic aspects of how fast and/or easy it is to write the code, rather than just the mere runtime.
R code:
gibbs_r <-function(N, thin) { mat <-matrix(nrow = N, ncol =2) x <- y <-0for (i in1:N) {for (j in1:thin) { x <-rgamma(1, 3, y * y +4) y <-rnorm(1, 1/ (x +1), 1/sqrt(2* (x +1))) } mat[i, ] <-c(x, y) } mat}
Actions to convert R to C++:
Add type declarations to all variables
Use ( instead of [ to index into the matrix
Subscript the results of rgamma and rnorm to convert from a vector into a scalar.
#include <Rcpp.h>usingnamespace Rcpp;// [[Rcpp::export]]NumericMatrix gibbs_cpp(int N,int thin){ NumericMatrix mat(N,2);double x =0, y =0;for(int i =0; i < N; i++){for(int j =0; j < thin; j++){ x = rgamma(1,3,1/(y * y +4))[0]; y = rnorm(1,1/(x +1),1/ sqrt(2*(x +1)))[0];} mat(i,0)= x; mat(i,1)= y;}return(mat);}
