06 - Rolling Functions
In cpp4r: Header-Only 'C++' and 'R' Interface

Notes

These examples were adapted from @vaughan24.
These functions ignore NA values for now. Adjustments for handling NA values are covered in a separate vignette. *R already provides efficient versions of the functions covered here. This is just to illustrate how to use C++ code.

Cumulative sum

Base R provides the cumsum() function to compute the cumulative sum of a vector:

cumsum(1:5)

> cumsum(1:5)
[1]  1  3  6 10 15

One possible C++ function to implement this is:

[[cpp4r::register]]
doubles cumsum_cpp(doubles x) {
  int n = x.size();
  writable::doubles out(n);

  out[0] = x[0];
  for (int i = 1; i < n; ++i) {
    out[i] = out[i - 1] + x[i];
  }
  return out;
}

Its R equivalent would be:

cumsum_r <- function(x) {
  n <- length(x)
  out <- numeric(n)
  out[1] <- x[1]
  for (i in 2:n) {
    out[i] <- out[i - 1] + x[i]
  }
  out
}

Add and document the functions, update the package as in the previous vignettes, and then compare the functions speed with:

# install.packages("microbenchmark")
library(microbenchmark)
library(mypkg)

set.seed(123) # for reproducibility
x <- rpois(1e6, lambda = 2) # 1,000,000 elements
x <- as.numeric(x) # integer to numeric, cumsum_cpp does not accept integers

cumsum(x[1:3])
cumsum_cpp(x[1:3])
cumsum_r(x[1:3])

microbenchmark(
  cumsum(x),
  cumsum_cpp(x),
  cumsum_r(x)
)

Unit: milliseconds
          expr   min    lq   mean median     uq    max neval cld
     cumsum(x)  2.07  2.40   3.40   2.61   4.59   7.04   100  a 
 cumsum_cpp(x)  2.35  2.62   3.84   2.85   4.93  11.00   100  a 
   cumsum_r(x) 91.20 98.50 114.00 109.00 124.00 233.00   100  b

Cumulative product

Base R provides the cumprod() function to compute the cumulative product of a vector:

cumprod(1:5)

> cumprod(1:5)
[1]   1   2   6  24 120

One possible C++ function to implement this is:

[[cpp4r::register]] doubles cumprod_cpp(doubles x) {
  int n = x.size();
  writable::doubles out(n);

  out[0] = x[0];
  for (int i = 1; i < n; ++i) {
    out[i] = out[i - 1] * x[i];
  }
  return out;
}

Its R equivalent would be:

cumprod_r <- function(x) {
  n <- length(x)
  out <- numeric(n)
  out[1] <- x[1]
  for (i in 2:n) {
    out[i] <- out[i - 1] * x[i]
  }
  out
}

To test the functions, you can run the following benchmark code in the R console:

microbenchmark(
  cumprod(x),
  cumprod_cpp(x),
  cumprod_r(x)
)

Unit: milliseconds
           expr   min    lq  mean median    uq    max neval cld
     cumprod(x)  2.50  2.66  3.19   2.79  3.04   6.08   100  a 
 cumprod_cpp(x)  2.24  2.45  2.88   2.55  2.71   6.99   100  a 
   cumprod_r(x) 79.70 83.30 93.70  93.50 99.60 125.00   100  b

Range of values

A simple example of the range() function in R is:

range(x)

> range(x)
[1]  0 13

One possible C++ function to implement this is:

[[cpp4r::register]]
doubles range_cpp(doubles x) {
  int n = x.size();
  double x1 = x[0], x2 = x[0];

  for (int i = 1; i < n; ++i) {
    x1 = std::min(x1, x[i]);
    x2 = std::max(x2, x[i]);
  }

  writable::doubles out(2);
  out[0] = x1;
  out[1] = x2;

  return out;
}

Write its R equivalent.

To verify the functions, you can run the following tests and benchmark code in the R console:

set.seed(123) # for reproducibility
x <- rnorm(10)

range(x)
range_cpp(x)

bigx <- list(
  rpois(2e6, lambda = 2),
  rpois(4e6, lambda = 2),
  rpois(8e6, lambda = 2)
)

bigx <- lapply(bigx, as.double)

lapply(
  bigx,
  function(x) {
    m <- summary(microbenchmark(range(x), range_cpp(x)))

    data.frame(
      expr = m[, 1], nobs = length(x), median = signif(m[, 5], 3),
      unit = attr(m, "unit")
    )
  }
)

[[1]]
          expr    nobs median         unit
1     range(x) 2000000  17.70 milliseconds
2 range_cpp(x) 2000000   5.42 milliseconds

[[2]]
          expr    nobs median         unit
1     range(x) 4000000   36.1 milliseconds
2 range_cpp(x) 4000000   11.0 milliseconds

[[3]]
          expr    nobs median         unit
1     range(x) 8000000   76.9 milliseconds
2 range_cpp(x) 8000000   22.0 milliseconds