---
title: "06 - Rolling Functions"
output:
  litedown::html_format:
    options:
      toc: true
      number_sections: true
vignette: >
  %\VignetteIndexEntry{06 - Rolling Functions}
  %\VignetteEngine{litedown::vignette}
  %\VignetteEncoding{UTF-8}
editor:
  markdown:
    wrap: sentence
bibliography: ["references.bib"]
---

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

```r
cumsum(1:5)
```

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

One possible C++ function to implement this is:

```cpp
[[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:

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

```r
# 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)
)
```

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

```r
cumprod(1:5)
```

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

One possible C++ function to implement this is:

```cpp
[[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:

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

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

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

```r
range(x)
```

```r
> range(x)
[1]  0 13
```

One possible C++ function to implement this is:

```cpp
[[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:

```r
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")
    )
  }
)
```

```r
[[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
```

## References