Intro to Rcpp

PHS 7045: Advanced Programming

George G. Vega Yon, Ph.D.

Intro

Before we start

You need to have Rcpp installed in your system:
```
install.packages("Rcpp")
```
You need to have a compiler
- Windows: You can download Rtools from here.
- MacOS: It is a bit complicated… Here are some options:
  - CRAN’s manual to get the clang, clang++, and gfortran compilers here.
  - A great guide by the coatless professor here

And that’s it!

R is great, but…

The problem:
- As we saw, R is very fast… once vectorized
- What to do if your model cannot be vectorized?
The solution: Use C/C++/Fotran! It works with R!
The problem to the solution: What R user knows any of those!?

R is great, but… (cont’d)

R has had an API (application programming interface) for integrating C/C++ code with R for a long time.
Unfortunately, it is not very straightforward

Enter Rcpp

One of the most important R packages on CRAN.
As of January 22, 2023, about 50% of CRAN packages depend on it (directly or not).
From the package description:

The ‘Rcpp’ package provides R functions as well as C++ classes which offer a seamless integration of R and C++

Why bother?

To draw ten numbers from a normal distribution with sd = 100.0 using R C API:

SEXP stats = PROTECT(R_FindNamespace(mkString("stats")));
SEXP rnorm = PROTECT(findVarInFrame(stats, install("rnorm")));
SEXP call = PROTECT(
  LCONS( rnorm, CONS(ScalarInteger(10), CONS(ScalarReal(100.0),
  R_NilValue))));
SET_TAG(CDDR(call),install("sd"));
SEXP res = PROTECT(eval(call, R_GlobalEnv));
UNPROTECT(4);
return res;

Using Rcpp:

Environment stats("package:stats");
Function rnorm = stats["rnorm"];
return rnorm(10, Named("sd", 100.0));

Example 1: Looping over a vector

#include<Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
NumericVector add1(NumericVector x) {
  NumericVector ans(x.size());
  for (int i = 0; i < x.size(); ++i)
    ans[i] = x[i] + 1;
  return ans;
}

add1(1:10)

 [1]  2  3  4  5  6  7  8  9 10 11

Example 1: Looping over a vector (vers 2)

Make it sweeter by adding some “sugar” (the Rcpp kind)

#include<Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
NumericVector add1Cpp(NumericVector x) {
  return x + 1;
}

add1Cpp(1:10)

 [1]  2  3  4  5  6  7  8  9 10 11

How much fast?

Compared to this:

add1R <- function(x) {
  for (i in 1:length(x))
    x[i] <- x[i] + 1
  x
}

microbenchmark::microbenchmark(add1R(1:1000), add1Cpp(1:1000))

Unit: microseconds
            expr    min     lq     mean  median      uq     max neval cld
   add1R(1:1000) 29.725 29.889 39.77246 30.3195 31.2625 911.266   100  a 
 add1Cpp(1:1000)  1.271  1.804  5.52475  1.9270  2.3575 320.333   100   b

C++ in R

Main differences between R and C++

One is compiled, and the other interpreted
Indexing objects: In C++ the indices range from 0 to (n - 1), whereas in R is from 1 to n.
All expressions end with a ; (optional in R).
In C++ object need to be declared, in R not (dynamic).

C++/Rcpp fundamentals: Types

Besides C-like data types (double, int, char, and bool), we can use the following types of objects with Rcpp:

Matrices: NumericMatrix, IntegerMatrix, LogicalMatrix, CharacterMatrix
Vectors: NumericVector, IntegerVector, LogicalVector, CharacterVector
And more!: DataFrame, List, Function, Environment

Parts of “an Rcpp program”

#include<Rcpp.h>
using namespace Rcpp
// [[Rcpp::export]]
NumericVector add1(NumericVector x) {
  NumericVector ans(x.size());
  for (int i = 0; i < x.size(); ++i)
    ans[i] = x[i] + 1;
  return ans;
}

Line by line, we see the following:

The #include<Rcpp.h> is similar to library(...) in R, it brings in all that we need to write C++ code for Rcpp.
using namespace Rcpp is somewhat similar to detach(...). This simplifies syntax. If we don’t include this, all calls to Rcpp members need to be explicit, e.g., instead of typing NumericVector, we would need to type Rcpp::NumericVector

The // starts a comment in C++, in this case, the // [[Rcpp::export]] comment is a flag Rcpp uses to “export” this C++ function to R.
It is the first part of the function definition. We are creating a function that returns a NumericVector, is called add1, has a single input element named x that is also a NumericVector.

Parts of “an Rcpp program” (cont’d)

#include<Rcpp.h>
using namespace Rcpp
// [[Rcpp::export]]
NumericVector add1(NumericVector x) {
  NumericVector ans(x.size());
  for (int i = 0; i < x.size(); ++i)
    ans[i] = x[i] + 1;
  return ans;
}

Here, we are declaring an object called ans, which is a NumericVector with an initial size equal to the size of x. Notice that .size() is called a “member function” of the x object, which is of class NumericVector.

We are declaring a for-loop (three parts):
1. int i = 0 We declare the variable i, an integer, and initialize it at 0.
2. i < x.size() This loop will end when i’s value is at or above the length of x.
3. ++i At each iteration, i will increment in one unit.
ans[i] = x[i] + 1 set the i-th element of ans equal to the i-th element of x plus 1.
return ans exists the function returning the vector ans.

C++/Rcpp fundamentals (cont’d 2)

Now, where to execute/run this?

You can use the sourceCpp function from the Rcpp package to run .cpp scripts (this is what I do most of the time).

There’s also cppFunction, which allows compiling a single function.

Write an R package that works with Rcpp.

For now, let’s use the first option.

Example running .cpp file

Imagine that we have the following file named norm.cpp

#| label: simple-file
#| cahe: true
#| echo: true
#include <Rcpp.h>

using namespace Rcpp;

// [[Rcpp::export]]
double normRcpp(NumericVector x) {
  
  return sqrt(sum(pow(x, 2.0)));
  
}

We can compile and obtain this function using this line Rcpp::sourceCpp("norm.cpp"). Once compiled, a function called normRcpp will be available in the current R session.

Now, get ready for some Rcpp action!

Your turn

Problem 1: Adding vectors

Using what you have just learned about Rcpp, write a function to add two vectors of the same length. Use the following template

#include <Rcpp.h>

using namespace Rcpp;

// [[Rcpp::export]]
NumericVector add_vectors([declare vector 1], [declare vector 2]) {
  
  ... magick ...
  
  return [something];
}

Now, we have to check for lengths. Use the stop function to make sure lengths match. Add the following lines in your code

if ([some condition])
  stop("an arbitrary error message :)");

Problem 2: Fibonacci series

Each element of the sequence is determined by the following:

\[ F(n) = \left\{\begin{array}{ll} n, & \mbox{ if }n \leq 1\\ F(n - 1) + F(n - 2), & \mbox{otherwise} \end{array}\right. \]

Using recursions, we can implement this algorithm in R as follows:

fibR <- function(n) {
  if (n <= 1)
    return(n)
  fibR(n - 1) + fibR(n - 2)
}

# Is it working?
c(
  fibR(0), fibR(1), fibR(2),
  fibR(3), fibR(4), fibR(5),
  fibR(6)
)

[1] 0 1 1 2 3 5 8

Now, let’s translate this code into Rcpp and see how much speed boost we get!

Problem 2: Fibonacci series (solution)

microbenchmark::microbenchmark(fibR(20), fibCpp(20))

Unit: microseconds
       expr      min        lq       mean    median       uq      max neval cld
   fibR(20) 2684.844 2771.6205 2991.03323 2828.3850 2992.139 4894.129   100  a 
 fibCpp(20)   18.204   19.0445   28.50812   20.3155   22.796  711.842   100   b

Fin

─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.4.1 (2024-06-14)
 os       macOS Sonoma 14.6.1
 system   aarch64, darwin23.4.0
 ui       unknown
 language (EN)
 collate  en_US.UTF-8
 ctype    en_US.UTF-8
 tz       America/Denver
 date     2024-08-19
 pandoc   3.2.1 @ /opt/homebrew/bin/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────
 package        * version  date (UTC) lib source
 cachem           1.1.0    2024-05-16 [1] CRAN (R 4.4.1)
 cli              3.6.3    2024-06-21 [1] CRAN (R 4.4.1)
 codetools        0.2-20   2024-03-31 [2] CRAN (R 4.4.1)
 devtools         2.4.5    2022-10-11 [1] CRAN (R 4.4.1)
 digest           0.6.36   2024-06-23 [1] CRAN (R 4.4.1)
 ellipsis         0.3.2    2021-04-29 [1] CRAN (R 4.4.1)
 evaluate         0.24.0   2024-06-10 [1] CRAN (R 4.4.1)
 fastmap          1.2.0    2024-05-15 [1] CRAN (R 4.4.1)
 fs               1.6.4    2024-04-25 [1] CRAN (R 4.4.1)
 glue             1.7.0    2024-01-09 [1] CRAN (R 4.4.1)
 htmltools        0.5.8.1  2024-04-04 [1] CRAN (R 4.4.1)
 htmlwidgets      1.6.4    2023-12-06 [1] CRAN (R 4.4.1)
 httpuv           1.6.15   2024-03-26 [1] CRAN (R 4.4.1)
 jsonlite         1.8.8    2023-12-04 [1] CRAN (R 4.4.1)
 knitr            1.47     2024-05-29 [1] CRAN (R 4.4.1)
 later            1.3.2    2023-12-06 [1] CRAN (R 4.4.1)
 lattice          0.22-6   2024-03-20 [2] CRAN (R 4.4.1)
 lifecycle        1.0.4    2023-11-07 [1] CRAN (R 4.4.1)
 magrittr         2.0.3    2022-03-30 [1] CRAN (R 4.4.1)
 MASS             7.3-60.2 2024-04-26 [2] CRAN (R 4.4.1)
 Matrix           1.7-0    2024-04-26 [2] CRAN (R 4.4.1)
 memoise          2.0.1    2021-11-26 [1] CRAN (R 4.4.1)
 microbenchmark   1.4.10   2023-04-28 [1] RSPM (R 4.4.1)
 mime             0.12     2021-09-28 [1] CRAN (R 4.4.1)
 miniUI           0.1.1.1  2018-05-18 [1] CRAN (R 4.4.1)
 multcomp         1.4-25   2023-06-20 [1] CRAN (R 4.4.1)
 mvtnorm          1.2-5    2024-05-21 [1] CRAN (R 4.4.1)
 pkgbuild         1.4.4    2024-03-17 [1] CRAN (R 4.4.1)
 pkgload          1.4.0    2024-06-28 [1] CRAN (R 4.4.1)
 profvis          0.3.8    2023-05-02 [1] CRAN (R 4.4.1)
 promises         1.3.0    2024-04-05 [1] CRAN (R 4.4.1)
 purrr            1.0.2    2023-08-10 [1] CRAN (R 4.4.1)
 R6               2.5.1    2021-08-19 [1] CRAN (R 4.4.1)
 Rcpp             1.0.12   2024-01-09 [1] CRAN (R 4.4.1)
 remotes          2.5.0    2024-03-17 [1] CRAN (R 4.4.1)
 rlang            1.1.4    2024-06-04 [1] CRAN (R 4.4.1)
 rmarkdown        2.27     2024-05-17 [1] CRAN (R 4.4.1)
 sandwich         3.1-0    2023-12-11 [1] CRAN (R 4.4.1)
 sessioninfo      1.2.2    2021-12-06 [1] CRAN (R 4.4.1)
 shiny            1.8.1.1  2024-04-02 [1] CRAN (R 4.4.1)
 stringi          1.8.4    2024-05-06 [1] CRAN (R 4.4.1)
 stringr          1.5.1    2023-11-14 [1] CRAN (R 4.4.1)
 survival         3.6-4    2024-04-24 [2] CRAN (R 4.4.1)
 TH.data          1.1-2    2023-04-17 [1] CRAN (R 4.4.1)
 urlchecker       1.0.1    2021-11-30 [1] CRAN (R 4.4.1)
 usethis          2.2.3    2024-02-19 [1] CRAN (R 4.4.1)
 vctrs            0.6.5    2023-12-01 [1] CRAN (R 4.4.1)
 xfun             0.45     2024-06-16 [1] CRAN (R 4.4.1)
 xtable           1.8-4    2019-04-21 [1] CRAN (R 4.4.1)
 yaml             2.3.8    2023-12-11 [1] CRAN (R 4.4.1)
 zoo              1.8-12   2023-04-13 [1] CRAN (R 4.4.1)

 [1] /opt/homebrew/lib/R/4.4/site-library
 [2] /opt/homebrew/Cellar/r/4.4.1/lib/R/library

──────────────────────────────────────────────────────────────────────────────