---
title: "Functional Programming with cppally"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Functional Programming with cppally}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r, include = FALSE}
required <- c("bench", "brio", "callr", "cli", "cpp11",
              "decor", "desc", "ggplot2", "glue", "purrr", "readr",
              "stringr", "tidyr", "usethis", "utils", "vctrs", "withr")
if (!all(vapply(required, requireNamespace, logical(1), quietly = TRUE))) {
  knitr::knit_exit()
}
```

This vignette is mainly a gallery of examples, serving to provide intuition behind the usage of cppally functionals, as well as showcasing their utility.

### A note on writing code in C++ versus R

Working with scalars demands a different mental model to the one R typically encourages.

In R we usually vectorise everything upfront by using vectors and vectorised operations everywhere possible. For R this is both cleaner and more efficient. In C++ the reverse is true. It is arguably cleaner to write functions at the scalar level, and then **vectorise after the fact**.

As will soon be demonstrated in this vignette, `pmap()` exists to bridge that final step of taking **scalar-based** functions and vectorising[^1] them.

[^1]: Vectorising in the R sense, not the SIMD (single-instruction-multiple-data) sense. In R, vectorising literally means taking a fully or partially scalar-based function and expanding the domain of the scalar inputs (or outputs) into the vector domain.

Let's load cppally in R and include cppally in our cpp or header file to get started.

```{r setup}
library(cppally)
```

``` cpp
// Include this if you are copying the code example-by-example
#include <cppally.hpp>
using namespace cppally;
```

```{r, include=FALSE}

set.seed(42)

ggplot2::theme_set(ggplot2::theme_grey())

not_cran <- function(){
 nzchar(Sys.getenv("CI")) || identical(Sys.getenv("NOT_CRAN"), "true")
}

n_threads <- if (not_cran()) 4L else 2L # CRAN allows at most 2 threads

# bench can report an extra gc row if a GC lands after the last iteration
mark <- function(..., env = parent.frame()) {
  out <- bench::mark(..., env = env)
  out$gc <- Map(function(gc, time) gc[seq_along(time), ], out$gc, out$time)
  out
}

# Helpers to compile all examples in debug mode
cpp_source <- function(..., code, debug = TRUE, env = parent.frame()){
  preamble <- c("#include <cppally.hpp>", "using namespace cppally;")
  code <- paste(c(preamble, code), collapse = "\n")
  tryCatch(
    cppally::cpp_source(
      debug = debug, env = env, code = code, ...), 
    error = function(e) {
      
      if (not_cran()){
        stop(e)
      }
      knitr::knit_exit()
    }
  )
}

# Helpers to source and display C++/R code
chunk_impl <- function(x, language){
  paste0("```", language, "\n", x, "\n```\n")
}
as_code_chunk <- function(x, language){
  cat(chunk_impl(x, language))
}
as_cpp_chunk <- function(x){
  as_code_chunk(x, "cpp")
}
```

```{r, include=FALSE}
# Compile necessary examples in one-go
# as it's faster when building the vignette

examples <- c(
cpp_sum = '
[[cppally::register]]
r_dbl cpp_sum(r_vector<r_dbl> x){
    return x.reduce([](auto a, auto b){ return a + b; });
}
',
cpp_sum2 = '
[[cppally::register]]
r_dbl cpp_sum2(r_vector<r_dbl> x){
    return x.reduce(std::plus<>{});
}
',
cpp_sum3 = '
[[cppally::register]]
r_dbl cpp_sum3(r_vector<r_dbl> x){
    return sum(x);
}
', 
cpp_gcd = '
[[cppally::register]]
r_int cpp_gcd(r_vector<r_int> x){
    return x.reduce([](auto acc, auto curr){
    
        auto res = cppally::gcd(acc, curr); // cppally has its own NA-aware gcd
        
        if ( (res == 1).is_true() ){
            return done(res);
        } else {
            return keep(res);
        }
        
    });
}
  ',
cpp_cumsum = '
[[cppally::register]]
r_vector<r_dbl> cpp_cumsum(r_vector<r_dbl> x){
    return x.cumulative_reduce(std::plus<>{});
}
  ', 
cpp_any_all_na = '
[[cppally::register]]
bool cpp_any_na(r_vector<r_dbl> x){
    return x.reduce([](auto, auto curr){ return is_na(curr) ? done(true) : keep(false); }, /*init = */ false);
}

[[cppally::register]]
bool cpp_all_na(r_vector<r_dbl> x){
    return x.reduce([](auto, auto curr){ return is_na(curr) ? keep(true) : done(false); }, /*init = */ true);
}
  ', 
cpp_pmax = '
[[cppally::register]]
r_vector<r_dbl> cpp_pmax(r_vector<r_dbl> x, r_vector<r_dbl> y){
    return pmap(
    
    /* fn = */ [](auto a, auto b){
    
      /* expr = */ return max(a, b);
        
    }, 
    
    /* vectors = */ x, y
  );
}
  ', 
cpp_if_else = '
template <RVector T>
[[cppally::register]]
T cpp_if_else(r_vector<r_lgl> condition, T if_true, T if_false, T if_na){
    return pmap(
      [](r_lgl condition_, auto yes, auto no, auto missing) {
        if (condition_.is_true()){
            return yes;
        } else if (condition_.is_false()){
            return no;
        } else {
            return missing;
        }
    }, 
    
    condition, if_true, if_false, if_na
  );
}
',
cpp_seq_along = '
template <RVector T>
[[cppally::register]]
r_vector<r_int> cpp_seq_along(T x){
    return pmap_with_index([](r_size_t i, auto){ // 2nd arg included so function can compile
        return as<r_int>(i) + 1; // R is 1-indexed
    }, x);
}
',
cpp_lag = '
[[cppally::register]]
r_vector<r_int> cpp_lag(r_vector<r_int> x, int k){
    return pmap_with_shift([&](auto a){
        return lag(a, k);
    }, x);
}
', 
cpp_diff = '
[[cppally::register]]
r_vector<r_dbl> cpp_diff(r_vector<r_dbl> x){
    return pmap_with_shift([&](auto a){
        return curr(a) - lag(a);
    }, x);
}
', 
cpp_in_place_abs = '
[[cppally::register]]
r_vector<r_dbl> cpp_in_place_abs(r_vector<r_dbl>& x){
    x.apply([](auto a){ return abs(a); });
    return x;
}
', 
cpp_in_place_lag = '
[[cppally::register]]
r_vector<r_dbl> cpp_in_place_lag(r_vector<r_dbl>& x, int k){
    x.shift(k);
    return x;
}
'
)

benchmark_examples <- c(
  
cpp_set_threads = '

[[cppally::register]]
void cpp_set_threads(int n){
  set_threads(n);
}
',
 cpp_parallel_pmax = '
 
// Allows us to set threads and automatically restore them once thread_guard is destroyed (even if R aborts)
struct thread_guard {

  int curr_threads;
  
  // Set threads and store previous threads
  thread_guard(int n) {
    curr_threads = get_threads();
    set_threads(n);
  }
  
  // Restore threads on destruction
  ~thread_guard() {
    set_threads(curr_threads);
  }
  
};
template <typename F>
decltype(auto) with_threads(int n, F&& f) {
  thread_guard guard(n);
  return std::forward<F>(f)();
}
 
template <typename T>
requires (any<T, r_dbl, r_str>) // doubles or strings
[[cppally::register]]
r_vector<T> cpp_parallel_pmax(r_vector<T> x, r_vector<T> y, int n_threads){
    return with_threads(n_threads, [&]{
    
      return pmap_parallel(
      [](const auto& a, const auto& b){
          return max(a, b);
      }, 
      x, y
    );
  });
}
  ',
cpp_sqrt = '

[[cppally::register]]
r_vector<r_dbl> cpp_sqrt(r_vector<r_dbl> x){
    return pmap_parallel_simd(
      [](auto v){
        return sqrt(v);
      }, 
    x
  );
}
', 
cpp_ops = '
template <RNumber T>
[[cppally::register]]
r_vector<T> cpp_add(r_vector<T> x, r_vector<T> y, r_vector<T> z){
  return x + y + z;
}
template <RNumber T>
[[cppally::register]]
r_vector<T> cpp_multiply(r_vector<T> x, r_vector<T> y, r_vector<T> z){
  return x * y * z;
}
',
group_names = '
// Helper to add group names, which are the unique values of group_data 
// so that each result is associated with its group name
r_vector<r_dbl> add_group_names(r_vector<r_dbl>& result, r_vector<r_str> group_data, groups g){
    r_vector<r_str> nms = group_names(group_data, g);
    result.set_names(nms);
    return result;
}
',
grouped_sum = '
[[cppally::register]]
r_vector<r_dbl> grouped_sum(r_vector<r_dbl> x, r_vector<r_str> g){
  
    groups grps = make_groups(g);
  
    r_vector<r_dbl> sums = apply_by_group(
      x, /*groups = */ grps,
      [](const r_vector<r_dbl>& chunk){
          return sum(chunk, /*na_rm = */ true);
      });

      return add_group_names(sums, g, grps);
}
',
grouped_median = '
[[cppally::register]]
r_vector<r_dbl> grouped_median(r_vector<r_dbl> x, r_vector<r_str> g){
  
    groups grps = make_groups(g);
  
    r_vector<r_dbl> medians = apply_by_group(
      x, grps,
      [](const r_vector<r_dbl>& chunk){
          return quantile(chunk, 0.5, /*na_rm = */ true);
      });

      return add_group_names(medians, g, grps);
}
',
cpp_grouped_apply = '
[[cppally::register]]
r_vector<r_dbl> cpp_grouped_apply(r_vector<r_dbl> x, r_vector<r_str> g, r_function fn){
    groups grps = make_groups(g);
    auto result = apply_by_group<r_dbl>(x, grps, fn);
    return add_group_names(result, g, grps);
}
'
  
)


cpp_source(code = paste(examples, collapse = "\n"), debug = TRUE)
cpp_source(code = paste(benchmark_examples, collapse = "\n"), debug = FALSE)
```

### pmap

`pmap` is a C++ variadic function that allows one to apply a function across corresponding elements of multiple vectors.

**Example:** vectorised binary max

```{r, echo=FALSE, comment="", results='asis'}
as_cpp_chunk(examples[["cpp_pmax"]])
```

```{r}
x <- c(10, 20, 30)
y <- c(10, 50, 0)
cpp_pmax(x, y)

# pmap also recycles vectors

cpp_pmax(x, 15)
```

**Example:** vectorised if else

```{r, echo=FALSE, comment="", results='asis'}
as_cpp_chunk(examples[["cpp_if_else"]])
```

```{r}
cpp_if_else(c(TRUE, FALSE, NA), "yes", "no", "missing")
```

`pmap_with_index()` is a variant that allows one to capture the index as we iterate along the vector

**Example:** Integer sequence along vector

```{r, echo=FALSE, comment="", results='asis'}
as_cpp_chunk(examples[["cpp_seq_along"]])
```

```{r}
cpp_seq_along(letters)
```

### reduce

`r_vector::reduce()` is a left-fold reduction functional that successively applies a binary function along the elements of the vector (from left-to-right). It allows for returning early and explicitly continuing by calling `done()` and `keep()`. The input function is typically a lambda, but can also be a callable (struct with `operator()`).

Example of summing a vector with `reduce`

```{r, echo=FALSE, comment="", results='asis'}
as_cpp_chunk(examples[["cpp_sum"]])
```

```{r}
cpp_sum(1:10)
```

We could have also passed the callable `std::plus<>{}`, which makes it even more readable.

```{r, echo=FALSE, comment="", results='asis'}
as_cpp_chunk(examples[["cpp_sum2"]])
```

```{r}
cpp_sum2(1:10)
```

Use `cumulative_reduce` to return a vector of all the intermediate results of the reduction.

```{r, echo=FALSE, comment="", results='asis'}
as_cpp_chunk(examples[["cpp_cumsum"]])
```

```{r}
cpp_cumsum(1:10)
```

While the above examples are useful for showing how to write a sum by hand, cppally provides `sum()` for free.

```{r, echo=FALSE, comment="", results='asis'}
as_cpp_chunk(examples[["cpp_sum3"]])
```

```{r}
cpp_sum3(c(10, 20, 30))
```

#### Returning early

To perform a reduction until a condition is met, use the helpers `done()` and `keep()`

**Example:** Check vector has any `NA` or all `NA`

```{r, echo=FALSE, comment="", results='asis'}
as_cpp_chunk(examples[["cpp_any_all_na"]])
```

```{r}
x <- c(1, 2, NA, 4, 5)
cpp_any_na(x)
cpp_all_na(x)
```

Notice that the two folds are identical - only `done`/`keep` and `init` are swapped. This is the general pattern for any/all-style predicates.

**Example:** greatest-common-divisor across integer vector. The trick here is to break when the result is 1 as `gcd(1, x) = 1` for any x.

```{r, echo=FALSE, comment="", results='asis'}
as_cpp_chunk(examples[["cpp_gcd"]])
```

```{r}
cpp_gcd(c(5L, 25L, 125L))
cpp_gcd(c(5L, 25L, 1L, 125L))
```

### Other pmap helpers

There are 2 core pmap functionals, with 9 variants in total.

**Core pmap functionals**

- `pmap` - Applies a function across elements.

- `pmap_with_index` - Like `pmap` but the first argument of the lambda must be an index.

**Other pmap variants**

- `pmap_parallel` - Like `pmap` but executed using multiple threads

- `pmap_simd` - Like `pmap` but executed under OpenMP SIMD instructions

- `pmap_parallel_simd` - Like `pmap` but multi-threaded and executed under OpenMP SIMD instructions

- `pmap_parallel_with_index` - Like `pmap_with_index` but executed using multiple threads

- `pmap_simd_with_index` - Like `pmap_with_index` but executed under OpenMP SIMD instructions

- `pmap_parallel_simd_with_index` - Like `pmap_with_index` but executed using multiple threads and under OpenMP SIMD instructions

- `pmap_with_shift` - A convenient wrapper around `pmap` to work with lagged values. Positional helpers `lag(k)`, `lead(k)`, `curr()` must be used to access lagged values. `lag_exists()` and `lead_exists()` can be used to test for out-of-bounds indexing. Out-of-bounds indexing returns `NA` by default, though this default can be changed in `lag()` and `lead()`.

#### Multi-threading safety

All `pmap` variants try very hard to avoid unsafe multi-threaded calls to R C API entry-points. Multi-threaded and/or SIMD execution only applies for `RVectorisable` types, meaning that for types like `r_str` which are not inherently thread-safe, execution is **always single-threaded** and **non-SIMD**, even if you request multiple threads or SIMD.

To illustrate this, let's revisit our `pmax` from earlier but this time using `pmap_parallel` and with a template that accepts `r_dbl` or `r_str` vectors.

```{r, echo=FALSE, comment="", results='asis'}
as_cpp_chunk(benchmark_examples[["cpp_parallel_pmax"]])
```

We can detect whether or not multiple threads were used by observing benchmark time.

```{r}
library(bench)
library(ggplot2)

x <- as.double(rpois(10^5, lambda = 10))
y <- as.double(rpois(10^5, lambda = 15))

(
  mark(
    single_threaded_pmax = cpp_parallel_pmax(x, y, n_threads = 1),
    multi_threaded_pmax = cpp_parallel_pmax(x, y, n_threads = n_threads)
  ) |> 
    autoplot(type = "violin") 
) +
  labs(title = "Numeric data")
```

We can see that the same benchmarks on character vector data yield identical results between the one requesting 1 thread and the one requesting `r n_threads` threads. Even though we request `r n_threads` threads, it still gets executed as single-threaded.

```{r}

x <- as.character(x)
y <- as.character(y)

(
  mark(
    single_threaded_pmax = cpp_parallel_pmax(x, y, n_threads = 1),
    multi_threaded_pmax = cpp_parallel_pmax(x, y, n_threads = n_threads)
  ) |> 
    autoplot(type = "violin") 
) +
  labs(title = "Character data")
```

### Lagged operations

Use `pmap_with_shift` to perform lagged operations. It checks for out-of-bounds access and therefore is safer than hand-writing it via `pmap_with_index`

```{r, echo=FALSE, comment="", results='asis'}
as_cpp_chunk(examples[["cpp_lag"]])
```

```{r}
# Lags
cpp_lag(1:10, k = 1)
cpp_lag(1:10, k = 2)
cpp_lag(1:10, k = 3)

# Leads
cpp_lag(1:10, k = -1)
cpp_lag(1:10, k = -2)
cpp_lag(1:10, k = -3)
```

`pmap_with_shift` has five helpers: `lag()`, `lead()`, `curr()`, `lag_exists()`, and `lead_exists()`. These are designed to assist in performing efficient lagged operations in a vectorised context, while maintaining readability.

**Example:** Lagged differencing

```{r, echo=FALSE, comment="", results='asis'}
as_cpp_chunk(examples[["cpp_diff"]])
```

```{r}
cpp_diff(1:10)
cpp_diff(seq(10, 100, by = 5))
```

### In-place functionals

To perform in-place transformations, use `r_vector::apply` as `pmap` always allocates a fresh vector and therefore cannot do in-place modification. `apply` comes in the same flavours as `pmap` - `apply_simd`, `apply_parallel`, `apply_parallel_simd`, and the `_with_index` variants.

```{r, echo=FALSE, comment="", results='asis'}
as_cpp_chunk(examples[["cpp_in_place_abs"]])
```

```{r}
x <- c(-20, -10)
cpp_in_place_abs(x)

x # Modified in-place
```

`r_vector::shift` is a helper which can shift an entire vector in-place. It takes shift `k` and an optional `fill_value` which defaults to `NA`.

```{r, echo=FALSE, comment="", results='asis'}
as_cpp_chunk(examples[["cpp_in_place_lag"]])
```

```{r}
x <- c(1, 2, 3, 4, 5)
cpp_in_place_lag(x, k = 1)
x # lagged in-place

# keep lagging until we run out of elements to lag
cpp_in_place_lag(x, k = 1)
cpp_in_place_lag(x, k = 1)
cpp_in_place_lag(x, k = 1)
cpp_in_place_lag(x, k = 1)
```

### Vectorised math

`pmap` also makes it easy to write vectorised math functions.

**Example:** Vectorised square-root

cppally provides a scalar version of `sqrt()`, which can be easily vectorised with `pmap`.

```{r, echo=FALSE, comment="", results='asis'}
as_cpp_chunk(benchmark_examples[["cpp_sqrt"]])
```

Here we are using `pmap_parallel_simd`, a variant of `pmap` that applies the supplied transformation under multiple threads and under an OpenMP SIMD directive. SIMD (single-instruction-multiple-data) is when the machine performs the same operation on multiple data points instead of one data point at time.

Let's benchmark this against `base::sqrt()`, but first let's use `r n_threads` threads for the rest of the examples in this vignette.

```{r, echo=FALSE, comment="", results='asis'}
as_cpp_chunk(benchmark_examples[["cpp_set_threads"]])
```

```{r}
cpp_set_threads(n_threads) # Set threads for the rest of the examples
```

```{r}
x <- rnorm(10^5, mean = 50)

(
  mark(
    base_sqrt = sqrt(x),
    cpp_sqrt = cpp_sqrt(x)
  ) |> 
    autoplot(type = "violin")
) + 
  labs(title = "R vs C++ vectorised square-root")
```

For more math functions, see `cppally/math/math.h`, a header containing a rich set of cppally math functions.

cppally also provides quite a few vectorised operators out-of-the-box. The vectorised operators currently defined:

binary: `+,-,*,/,+=,-=,*=,/=,==,<=,<,>=,>,|,&`

unary: `!,-`

When it comes to arithmetic operations, both cppally and R are heavily optimised to use in-place modification where possible. This greatly improves performance when multiple arithmetic operations are chained one after the other.

Since R is heavily optimised in this case, any performance gains cppally makes over R are likely to come from using multiple threads.

```{r, echo=FALSE, comment="", results='asis'}
as_cpp_chunk(benchmark_examples[["cpp_ops"]])
```

```{r}
a <- sample(c(-1L, NA_integer_, 0L, 1L, -100L, 100L), 10^5, replace = TRUE)
b <- sample(c(NA_real_, 0, -1/3, 1/3, -10^9, 10^9), 10^5, replace = TRUE)

results <-
 mark(
   `a + a + a (base w/integers)` = a + a + a,
   `a * a * a (base w/integers)` = a * a * a,
   `b + b + b (base w/doubles)` = b + b + b,
   `b * b * b (base w/doubles)` = b * b * b,
   
   `a + a + a (cppally w/integers)` = cpp_add(a, a, a),
   `a * a * a (cppally w/integers)` = cpp_multiply(a, a, a),
   `b + b + b (cppally w/doubles)` = cpp_add(b, b, b),
   `b * b * b (cppally w/doubles)` = cpp_multiply(b, b, b),
   check = FALSE
 )

results |> 
  autoplot(type = "violin") +
  geom_violin(aes(fill = ifelse(grepl("base", expression), "base", "cppally"))) +
  labs(title = "cppally vs R arithmetic operations") + 
  scale_y_discrete(limits = rev) +
  theme(legend.title = element_blank())
  
```

From the results we can see that integer operations are faster with cppally. This is partly due to multiple threads being used, but also because cppally integer arithmetic is carefully written to be SIMD-friendly. It propagates `NA` values and also returns `NA` on integer-overflow, so it is just as safe as R. The difference is that the arithmetic is written to be as branchless as possible, allowing SIMD instructions to be executed. You can see the internal implementation in `scalar/arithmetic_ops.h`.

### Grouped operations

In modern analyses it has become more common to see grouped operations. Packages like 'dplyr' (via `group_by()`) and 'data.table' (via `by` keyword) natively support applying R functions by group.

cppally provides a similar group-by-and-apply functionality 
via `apply_by_group()`: a C++ functional that accepts 
a vector `x` of data, groups `g`, and a vector-based function `fn` 
which is applied to `x` on a group-by-group basis using `g`.

`fn` can be a C++ function or an `r_function` (R function). We will first
use `apply_by_group()` with an `r_function`, then with a C++ function, and
compare both to base R.

#### Applying an R function by group

To produce the groups, use `make_groups()`, a function that converts a 
vector into group IDs and other group metadata, wrapped in a C++ class `groups`.

```{r, echo=FALSE, comment="", results='asis'}
as_cpp_chunk(benchmark_examples[["cpp_grouped_apply"]])
```

The group names are added using a helper `add_group_names()`

```{r, echo=FALSE, comment="", results='asis'}
as_cpp_chunk(benchmark_examples[["group_names"]])
```

For comparison, a base R split-and-apply equivalent

```{r}
r_grouped_apply <- function(x, g, fn){
  vapply(split(x, factor(g, levels = unique(g))), fn, 0)
}
```

#### Applying a C++ function by group

Instead of passing R's `sum()`, we can pass a lambda calling a C++
version, `cppally::sum()`.

```{r, echo=FALSE, comment="", results='asis'}
as_cpp_chunk(benchmark_examples[["grouped_sum"]])
```

#### Benchmarks

`base::sum()` is a primitive, which means it calls C code directly, hence it
should provide a good baseline for what can be achieved with a base R
split-and-apply method.

```{r}
cpp_set_threads(1) # Run these benchmarks using a single thread
n <- 10^5
x <- rnorm(n)

# 3 groups
g <- sample(c("a", "b", "c"), n, replace = TRUE)
(
  mark(
  base_grouped_sum = r_grouped_apply(x, g, sum),
  cppally_grouped_r_sum = cpp_grouped_apply(x, g, sum),
  cppally_grouped_cpp_sum = grouped_sum(x, g),
  min_iterations = 10
) |> 
    autoplot(type = "violin")
) + 
  labs(title = "R vs cppally sum by group - 3 groups") +
  bench::scale_x_bench_time(base = 2)

# 1e04 groups
g <- as.character(sample.int(10^4, n, replace = TRUE))
(
  mark(
  base_grouped_sum = r_grouped_apply(x, g, sum),
  cppally_grouped_r_sum = cpp_grouped_apply(x, g, sum),
  cppally_grouped_cpp_sum = grouped_sum(x, g),
  min_iterations = 10
) |> 
    autoplot(type = "violin")
) + 
  labs(title = "R vs cppally sum by group - 10k groups") +
  bench::scale_x_bench_time(base = 2)

# 1e05 groups
g <- as.character(sample.int(10^5, n, replace = TRUE))
(
  mark(
  base_grouped_sum = r_grouped_apply(x, g, sum),
  cppally_grouped_r_sum = cpp_grouped_apply(x, g, sum),
  cppally_grouped_cpp_sum = grouped_sum(x, g),
  min_iterations = 10
) |> 
    autoplot(type = "violin")
) + 
  labs(title = "R vs cppally sum by group - 100k groups") +
  bench::scale_x_bench_time(base = 2)
```

Passing R's `sum()` performs similarly to `cppally::sum()` for small numbers 
of groups. Most of the time is spent grouping and moving data, which is shared by both, and calling a primitive like `sum()` from C++ adds little overhead per group.

#### By-group median

Non-primitive R functions like `median()` are much more expensive to call,
and that cost is paid once per group. cppally's `quantile(prob = 0.5)` gives the same result as `median()`, so we can use it for the comparison.

```{r, echo=FALSE, comment="", results='asis'}
as_cpp_chunk(benchmark_examples[["grouped_median"]])
```

```{r}
# 1e04 groups
g <- as.character(sample.int(10^4, n, replace = TRUE))
(
  mark(
  base_grouped_median = r_grouped_apply(x, g, median),
  cppally_grouped_r_median = cpp_grouped_apply(x, g, median),
  cppally_grouped_cpp_median = grouped_median(x, g),
  min_iterations = 3
) |> 
    autoplot(type = "violin")
) + 
  labs(title = "R vs cppally median by group - 10k groups") +
  bench::scale_x_bench_time(base = 2)
```

As we move from simpler functions like `sum()` to more complex ones like 
`median()`, we see that the calculation over many groups starts to become 
quite slow with R functions. There is a large difference between `cppally_grouped_r_median` and `cppally_grouped_cpp_median`, indicating that 
most of the time difference is spent in calling `stats::median()` itself.

Passing efficient C++ functions is hence more suitable for function calls on many groups.