25_Rewriting_R_code_in_C++.Rmd

# Rewriting R code in C++

**Learning objectives:**

-   Learn to improve performance by rewriting bottlenecks in C++

-   Introduction to the [{Rcpp} package](https://www.rcpp.org/)

## Introduction

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

<center>Like how?</center>

<center> </center>

<center>![](https://media.giphy.com/media/vLyZk5CJo12Wk/giphy.gif)</center>

 

## Getting started with C++

```{r warning=FALSE}
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:

```{r}
cppFunction('int add(int x, int y, int z) {
  int sum = x + y + z;
  return sum;
}')
# add works like a regular R function
add

add(1, 2, 3)
```

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.

 

## Example with scalar input and output {-}

```{r}
signR <- function(x) {
  if (x > 0) {
    1
  } else if (x == 0) {
    0
  } else {
    -1
  }
}

a <- -0.5
b <- 0.5
c <- 0
signR(c)
```

Translation:

```{r}
cppFunction('int signC(int x) {
  if (x > 0) {
    return 1;
  } else if (x == 0) {
    return 0;
  } else {
    return -1;
  }
}')
```

* Note that the `if` syntax is identical! Not everything is different!

## Vector Input, Scalar output:{-}

```{r}
sumR <- function(x) {
  total <- 0
  for (i in seq_along(x)) {
    total <- total + x[i]
  }
  total
}

x<- runif(100)
sumR(x)
```

Translation:

```{r}
cppFunction('double sumC(NumericVector x) {
  int n = x.size();
  double total = 0;
  for(int i = 0; i < n; ++i) {
    total += x[i];
  }
  return total;
}')
```

Some observations:

-   vector indices *start at 0*
-   The for statement has a different syntax: for(init; check; increment)
-   Methods are called with `.`
-   `total += x[i]` is equivalent to `total = total + x[i]`.
-   other in-place operators are `-=`, `*=`, `and /=`


To check for the fastest way we can use:

```{r eval=FALSE}
?bench::mark
```

```{r}
x <- runif(1e3)
bench::mark(
  sum(x),
  sumC(x),
  sumR(x)
)
```

## Vector input and output {-}

```{r}
pdistR <- function(x, ys) {
  sqrt((x - ys) ^ 2)
}
```

```{r}
cppFunction('NumericVector pdistC(double x, NumericVector ys) {
  int n = ys.size();
  NumericVector out(n);

  for(int i = 0; i < n; ++i) {
    out[i] = sqrt(pow(ys[i] - x, 2.0));
  }
  return out;
}')
```

Note:   uses `pow()`, not `^`, for exponentiation

```{r}
y <- runif(1e6)
bench::mark(
  pdistR(0.5, y),
  pdistC(0.5, y)
)[1:6]
```

## Source your C++ code {-}

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". 

```{Rcpp}
#include <Rcpp.h>
using namespace 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;
}

/*** R
x <- 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?

```{r}
x <- runif(1e5)
bench::mark(
  mean(x),
  meanC(x)
)
```



## Data frames, functions, and attributes

### Lists and Dataframes {-}

Contrived example to illustrate how to access a dataframe from c++:

```{Rcpp}
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
double mpe(List mod) {
  if (!mod.inherits("lm")) stop("Input must be a linear model");

  NumericVector resid = as<NumericVector>(mod["residuals"]);
  NumericVector fitted = as<NumericVector>(mod["fitted.values"]);

  int n = resid.size();
  double err = 0;
  for(int i = 0; i < n; ++i) {
    err += resid[i] / (fitted[i] + resid[i]);
  }
  return err / n;
}
```

```{r}
mod <- lm(mpg ~ wt, data = mtcars)
mpe(mod)
```

- Note that you must *cast* the values to the required type. C++ needs to know the types in advance.

### Functions {-}

```{Rcpp}
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
RObject callWithOne(Function f) {
  return f(1);
}
```


```{r}
callWithOne(function(x) x + 1)
```


* Other values can be accessed from c++ including

   * attributes (use: `.attr()`. Also `.names()` is alias for name attribute.
   * `Environment`, `DottedPair`, `Language`, `Symbol` , etc. 

## Missing values

### Missing values behave differently for C++ scalers{-}

* Scalar NA's in Cpp : `NA_LOGICAL`, `NA_INTEGER`, `NA_REAL`, `NA_STRING`.

* Integers (`int`) stores R NA's as the smallest integer. Better to use length 1 `IntegerVector`
* Doubles use IEEE 754 NaN , which behaves a bit differently for logical expressions (but ok for math expressions). 

```{r}
evalCpp("NA_REAL || FALSE")
```

* Strings are a class from Rcpp, so they handle missing values fine.

* `bool` can only hold two values, so be careful. Consider using vectors of length 1 or coercing to `int`


### Vectors

* Vectors are all type introduced by RCpp and know how to handle missing values if you use the specific type for that vector.

```{Rcpp}
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
List missing_sampler() {
  return List::create(
    NumericVector::create(NA_REAL),
    IntegerVector::create(NA_INTEGER),
    LogicalVector::create(NA_LOGICAL),
    CharacterVector::create(NA_STRING)
  );
}
```

```{r}
str(missing_sampler())
```

## 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:

```{Rcpp}

#include <algorithm>
#include <Rcpp.h>

using namespace 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;
}
```

```{r}
square_C(c(1.0,2.0,3.0))
```
```{r}
square_C_it(c(1.0,2.0,3.0))
```

## 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`:

```{Rcpp}
#include <Rcpp.h>
using namespace 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;
}
```


```{r}
res = tableC(c(1,1,2,1,4,5))
res
```

* 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](https://en.cppreference.com/w/cpp/container) at `cppreference`

## Case Studies

![Case Study](images/case_study.jpg)

Real life uses of C++ to replace slow R code.


## Case study 1: Gibbs sampler {-}

The [Gibbs sampler](https://en.wikipedia.org/wiki/Gibbs_sampling) 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. 

[Example blogged by Dirk Eddelbuettel](https://dirk.eddelbuettel.com/blog/2011/07/14/), the R and C++ code is very similar but runs about 20 times faster.

> "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.


<center>![](https://media.giphy.com/media/13GIgrGdslD9oQ/giphy.gif)</center>


R code:

```{r}
gibbs_r <- function(N, thin) {
  mat <- matrix(nrow = N, ncol = 2)
  x <- y <- 0

  for (i in 1:N) {
    for (j in 1: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.

```{Rcpp}
#include <Rcpp.h>
using namespace 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);
}
```

Checking who's best:

```{r}
bench::mark(
  gibbs_r(100, 10),
  gibbs_cpp(100, 10),
  check = FALSE
)
```

## Case study 2: predict a model response from three inputs {-}

[Rcpp is smoking fast for agent based models in data frames](https://gweissman.github.io/post/rcpp-is-smoking-fast-for-agent-based-models-in-data-frames/) by Gary Weissman, MD, MSHP.

Starts with this code:

```{r}
vacc1a <- function(age, female, ily) {
  p <- 0.25 + 0.3 * 1 / (1 - exp(0.04 * age)) + 0.1 * ily
  p <- p * if (female) 1.25 else 0.75
  p <- max(0, p)
  p <- min(1, p)
  p
}
```

R code with a for loop:

```{r}
vacc1 <- function(age, female, ily) {
  n <- length(age)
  out <- numeric(n)
  for (i in seq_len(n)) {
    out[i] <- vacc1a(age[i], female[i], ily[i])
  }
  out
}
```

Vectorized R code:

```{r}
vacc2 <- function(age, female, ily) {
  p <- 0.25 + 0.3 * 1 / (1 - exp(0.04 * age)) + 0.1 * ily
  p <- p * ifelse(female, 1.25, 0.75)
  p <- pmax(0, p)
  p <- pmin(1, p)
  p
}
```

C++:

```{Rcpp}
#include <Rcpp.h>
using namespace Rcpp;

double vacc3a(double age, bool female, bool ily){
  double p = 0.25 + 0.3 * 1 / (1 - exp(0.04 * age)) + 0.1 * ily;
  p = p * (female ? 1.25 : 0.75);
  p = std::max(p, 0.0);
  p = std::min(p, 1.0);
  return p;
}

// [[Rcpp::export]]
NumericVector vacc3(NumericVector age, LogicalVector female, 
                    LogicalVector ily) {
  int n = age.size();
  NumericVector out(n);

  for(int i = 0; i < n; ++i) {
    out[i] = vacc3a(age[i], female[i], ily[i]);
  }

  return out;
}
```

Sample data:

```{r}
n <- 1000
age <- rnorm(n, mean = 50, sd = 10)
female <- sample(c(T, F), n, rep = TRUE)
ily <- sample(c(T, F), n, prob = c(0.8, 0.2), rep = TRUE)

stopifnot(
  all.equal(vacc1(age, female, ily), vacc2(age, female, ily)),
  all.equal(vacc1(age, female, ily), vacc3(age, female, ily))
)
```

<center>**Who's faster?**</center>
<center>![](https://media.giphy.com/media/l41JGlWa1xOjJSsV2/giphy.gif)</center>

```{r}
bench::mark(
  vacc1 = vacc1(age, female, ily),
  vacc2 = vacc2(age, female, ily),
  vacc3 = vacc3(age, female, ily)
)
```

## Resources

-   [Rcpp: Seamless R and C++ Integration](https:\\Rcpp.org)
-   [cpp-tutorial](https://www.learncpp.com) is often recommended. Lots of ads though!
-   [cpp-reference](https://en.cppreference.com/w/cpp)
-   [C++20 for Programmers](https://www.pearson.com/en-us/subject-catalog/p/c20-for-programmers-an-objects-natural-approach/P200000000211/9780137570461) is a newer book that covers modern c++ for people who know programming in another language.
 
## Op Success!

![Congrats!](images/we-did-it-celebration-meme.jpg)


## Meeting Videos

### Cohort 1

`r knitr::include_url("https://www.youtube.com/embed/2JDeacWl1DM")`

`r knitr::include_url("https://www.youtube.com/embed/sLWCelHpcqc")`

### Cohort 2

`r knitr::include_url("https://www.youtube.com/embed/rQwOosOJpaY")`

### Cohort 3

`r knitr::include_url("https://www.youtube.com/embed/ZWdIeR1jK9Q")`

### Cohort 4

`r knitr::include_url("https://www.youtube.com/embed/_K8DKF3Fzes")`

### Cohort 5

`r knitr::include_url("https://www.youtube.com/embed/nske4iqsgh0")`

### Cohort 6

`r knitr::include_url("https://www.youtube.com/embed/hyVK08jXiYw")`

<details>

<summary>Meeting chat log</summary>
```
00:10:13	Arthur Shaw:	Did things freeze for anyone else?
00:55:40	Federica Gazzelloni:	https://en.cppreference.com/w/cpp/container
00:57:44	Federica Gazzelloni:	https://dirk.eddelbuettel.com/blog/2011/07/14/
01:07:33	Trevin:	I don’t have experience
01:07:54	Oluwafemi Oyedele:	Same here!!!
01:11:57	Arthur Shaw:	Does anyone know any packages that use C++? The one that comes to mind for me is haven, which uses a C++ library
01:12:30	Trevin:	When I was looking, one that stood out to me was rstan
01:13:02	Arthur Shaw:	Reacted to "When I was looking, ..." with 👍
```
</details>

### Cohort 7

`r knitr::include_url("https://www.youtube.com/embed/Luu7JsixQgY")`

<details>

<summary>Meeting chat log</summary>
```
00:43:02	Gus Lipkin:	I think I found the definition for `mean`

An R call goes to *a which then calls the C function *b

*a: https://github.com/wch/r-source/blob/trunk/src/library/base/R/mean.R
*b: https://github.com/wch/r-source/blob/trunk/src/library/stats/src/cov.c#L207

It looks like the second pass only happens if `R_FINITE(mean_from_first_pass)` which tries to call `isfinite` from C++ and if it’s not there, it’ll make sure it is a number and is not positive or negative infinity.
00:49:55	Gus Lipkin:	I feel bad for dropping in on the last chapter and getting Collin’s thanks 😅 I wish I’d joined sooner.
```
</details>