# Benchmarks

Each program runs as plain R on GNU R, as the same plain R on ry, and typed in
[strict mode](/type-checking/strict-mode.md) on ry. All three compute the same result.

## CPU

| program | GNU R | ry | ry, typed and strict | speed-up |
| --- | ---: | ---: | ---: | ---: |
| `array_dot` | 0.278 s | 0.043 s | 0.047 s | 6.5x |
| `array_filter` | 0.246 s | 0.043 s | 0.038 s | 5.7x |
| `array_map` | 0.385 s | 0.054 s | 0.067 s | 7.1x |
| `array_sum` | 0.175 s | 0.033 s | 0.033 s | 5.3x |
| `closures` | 1.100 s | 0.570 s | 0.607 s | 1.9x |
| `drawdown` | 0.335 s | 0.074 s | 0.079 s | 4.5x |
| `fib` | 1.312 s | 0.010 s | 0.010 s | 131.2x |
| `log_rank` | 0.128 s | 0.137 s | 0.136 s | 0.9x |
| `loop_sum` | 0.312 s | 0.069 s | 0.068 s | 4.5x |
| `vasicek` | 0.055 s | 0.040 s | 0.040 s | 1.4x |
| `vector_arith` | 0.290 s | 0.242 s | 0.249 s | 1.2x |

## Memory

| program | GNU R | ry | ry, typed and strict | ry / GNU R |
| --- | ---: | ---: | ---: | ---: |
| `array_dot` | 81 MiB | 103 MiB | 124 MiB | 1.3x |
| `array_filter` | 81 MiB | 102 MiB | 123 MiB | 1.3x |
| `array_map` | 94 MiB | 123 MiB | 141 MiB | 1.3x |
| `array_sum` | 81 MiB | 102 MiB | 123 MiB | 1.3x |
| `closures` | 99 MiB | 145 MiB | 148 MiB | 1.5x |
| `drawdown` | 73 MiB | 133 MiB | 141 MiB | 1.8x |
| `fib` | 54 MiB | 43 MiB | 43 MiB | 0.8x |
| `log_rank` | 107 MiB | 56 MiB | 58 MiB | 0.5x |
| `loop_sum` | 50 MiB | 82 MiB | 106 MiB | 1.6x |
| `vasicek` | 106 MiB | 64 MiB | 70 MiB | 0.6x |
| `vector_arith` | 109 MiB | 157 MiB | 164 MiB | 1.4x |

Peak memory is the peak resident set size of the whole process, R's startup included.

Recursion and loops gain the most, while vectorized builtins already run in C on GNU R. Types do not
make code faster yet, and ry usually peaks higher in memory than GNU R.

## The programs

Each program exists in `tools/bench/why/` as plain R and typed.

### array_dot

Element loops over a million doubles: a dot product of `x[[i]] * y[[i]]`. The typed program stores the same values in Arrays, which hold no NA and no attributes.

Plain R:

```r
array_dot <- function(x, y, times) {
  total <- 0
  n <- length(x)
  for (k in seq_len(times)) {
    for (i in seq_len(n)) total <- total + x[[i]] * y[[i]]
  }
  total
}
x <- seq(0, 1, length.out = 1000000L)
y <- rev(x)
measure("array_dot", function() array_dot(x, y, 5L))
```

Typed:

```r
array_dot <- function(x: Array<double>, y: Array<double>, times: integer) -> double {
  total <- 0
  n <- List.len(x)
  for (k in seq_len(times)) {
    for (i in seq_len(n)) total <- total + x[[i]] * y[[i]]
  }
  total
}
values: double[] <- seq(0, 1, length.out = 1000000L)
reversed: double[] <- rev(values)
x <- as.Array(values)
y <- as.Array(reversed)
measure("array_dot", function() array_dot(x, y, 5L))
```

### array_filter

Element loops over a million doubles: a filter counting the items above one half. The typed program stores the same values in Arrays, which hold no NA and no attributes.

Plain R:

```r
array_filter <- function(x, y, times) {
  total <- 0
  n <- length(x)
  for (k in seq_len(times)) {
    for (i in seq_len(n)) if (x[[i]] > 0.5) total <- total + 1
  }
  total
}
x <- seq(0, 1, length.out = 1000000L)
y <- rev(x)
measure("array_filter", function() array_filter(x, y, 5L))
```

Typed:

```r
array_filter <- function(x: Array<double>, y: Array<double>, times: integer) -> double {
  total <- 0
  n <- List.len(x)
  for (k in seq_len(times)) {
    for (i in seq_len(n)) if (x[[i]] > 0.5) total <- total + 1
  }
  total
}
values: double[] <- seq(0, 1, length.out = 1000000L)
reversed: double[] <- rev(values)
x <- as.Array(values)
y <- as.Array(reversed)
measure("array_filter", function() array_filter(x, y, 5L))
```

### array_map

Element loops over a million doubles: a map writing `out[[i]] <- x[[i]] * 2 + 1` into a preallocated vector. The typed program stores the same values in Arrays, which hold no NA and no attributes.

Plain R:

```r
array_map <- function(x, y, times) {
  total <- 0
  n <- length(x)
  out <- numeric(n)
  for (k in seq_len(times)) {
    for (i in seq_len(n)) out[[i]] <- x[[i]] * 2 + 1
    total <- total + sum(out)
  }
  total
}
x <- seq(0, 1, length.out = 1000000L)
y <- rev(x)
measure("array_map", function() array_map(x, y, 5L))
```

Typed:

```r
array_map <- function(x: Array<double>, y: Array<double>, times: integer) -> double {
  total <- 0
  n <- List.len(x)
  out <- as.Array(numeric(n))
  for (k in seq_len(times)) {
    for (i in seq_len(n)) out[[i]] <- x[[i]] * 2 + 1
    total <- total + List.sum(out)
  }
  total
}
values: double[] <- seq(0, 1, length.out = 1000000L)
reversed: double[] <- rev(values)
x <- as.Array(values)
y <- as.Array(reversed)
measure("array_map", function() array_map(x, y, 5L))
```

### array_sum

Element loops over a million doubles: a sum of `x[[i]]`. The typed program stores the same values in Arrays, which hold no NA and no attributes.

Plain R:

```r
array_sum <- function(x, y, times) {
  total <- 0
  n <- length(x)
  for (k in seq_len(times)) {
    for (i in seq_len(n)) total <- total + x[[i]]
  }
  total
}
x <- seq(0, 1, length.out = 1000000L)
y <- rev(x)
measure("array_sum", function() array_sum(x, y, 5L))
```

Typed:

```r
array_sum <- function(x: Array<double>, y: Array<double>, times: integer) -> double {
  total <- 0
  n <- List.len(x)
  for (k in seq_len(times)) {
    for (i in seq_len(n)) total <- total + x[[i]]
  }
  total
}
values: double[] <- seq(0, 1, length.out = 1000000L)
reversed: double[] <- rev(values)
x <- as.Array(values)
y <- as.Array(reversed)
measure("array_sum", function() array_sum(x, y, 5L))
```

### closures

Closures and lapply: small closures called in a loop, then lapply and vapply over a sequence.

Plain R:

```r
closures <- function(n) {
  add <- function(a, b) a + b
  scale <- 3
  times_scale <- function(a) a * scale
  s <- 0
  for (i in 1:n) s <- add(s, times_scale(i))
  squares <- lapply(seq_len(n), function(x) x * (x + 0.5))
  halves <- vapply(seq_len(n), function(x) x / 2, numeric(1))
  squared <- unlist(squares)
  s + sum(squared) + sum(halves)
}
measure("closures", function() closures(300000L))
```

Typed:

```r
closures <- function(n: integer) -> double {
  add <- function(a: double, b: double) -> double a + b
  scale <- 3
  times_scale <- function(a: double) -> double a * scale
  s <- 0
  for (i in 1:n) s <- add(s, times_scale(i))
  squares <- lapply(seq_len(n), function(x: integer) -> double x * (x + 0.5))
  halves <- vapply(seq_len(n), function(x: integer) -> double x / 2, numeric(1))
  squared: double[] <- unlist(squares)
  s + sum(squared) + sum(halves)
}
measure("closures", function() closures(300000L))
```

### drawdown

A hot path from the finance code base: maximum drawdown and EWMA volatility, scalar state carried through one pass over a return series.

Plain R:

```r
risk_report <- function(returns, lambda = 0.94) {
  wealth <- 1
  peak <- 1
  worst <- 0
  variance <- returns[1] * returns[1]
  for (i in seq_along(returns)) {
    r <- returns[i]
    wealth <- wealth * (1 + r)
    if (wealth > peak) peak <- wealth
    drawdown <- (peak - wealth) / peak
    if (drawdown > worst) worst <- drawdown
    variance <- lambda * variance + (1 - lambda) * r * r
  }
  c(max_drawdown = worst, volatility = sqrt(variance * 252))
}
returns <- sin(seq_len(2000000L) * 0.37) * 0.01 + 0.0002
measure("drawdown", function() sum(risk_report(returns)))
```

Typed:

```r
risk_report <- function(returns: double[], lambda: double = 0.94) -> double[] {
  wealth <- 1
  peak <- 1
  worst <- 0
  variance <- returns[1] * returns[1]
  for (i in seq_along(returns)) {
    r <- returns[i]
    wealth <- wealth * (1 + r)
    if (wealth > peak) peak <- wealth
    drawdown <- (peak - wealth) / peak
    if (drawdown > worst) worst <- drawdown
    variance <- lambda * variance + (1 - lambda) * r * r
  }
  c(max_drawdown = worst, volatility = sqrt(variance * 252))
}
returns <- sin(seq_len(2000000L) * 0.37) * 0.01 + 0.0002
measure("drawdown", function() sum(risk_report(returns)))
```

### fib

Recursion: about 2.7 million calls of a small function.

Plain R:

```r
fib <- function(n) if (n < 2) n else fib(n - 1) + fib(n - 2)
measure("fib", function() fib(30))
```

Typed:

```r
fib <- function(n: double) -> double if (n < 2) n else fib(n - 1) + fib(n - 2)
measure("fib", function() fib(30))
```

### log_rank

A hot path from the biostatistics code base: the log-rank statistic, a loop over event times with whole-vector comparisons in its body.

Plain R:

```r
log_rank <- function(data) {
  times <- sort(unique(data$time[data$event == 1L]))
  observed <- 0
  expected <- 0
  variance <- 0
  for (t in times) {
    risk <- data$time >= t
    n <- sum(risk)
    n1 <- sum(risk & data$arm == "treated")
    d <- sum(data$time == t & data$event == 1L)
    d1 <- sum(data$time == t & data$event == 1L & data$arm == "treated")
    observed <- observed + d1
    expected <- expected + d * n1 / n
    if (n > 1) variance <- variance + d * (n1 / n) * (1 - n1 / n) * (n - d) / (n - 1)
  }
  (observed - expected)^2 / variance
}
subjects <- 4000L
data <- list(
  time = (seq_len(subjects) * 7919L) %% 1000L + 1L,
  event = as.integer(seq_len(subjects) %% 3L != 0L),
  arm = ifelse(seq_len(subjects) %% 2L == 0L, "treated", "control")
)
measure("log_rank", function() log_rank(data))
```

Typed:

```r
log_rank <- function(data: list{time: integer[], event: integer[], arm: character[]}) -> double {
  times <- sort(unique(data$time[data$event == 1L]))
  observed <- 0
  expected <- 0
  variance <- 0
  for (t in times) {
    risk <- data$time >= t
    n <- sum(risk)
    n1 <- sum(risk & data$arm == "treated")
    d <- sum(data$time == t & data$event == 1L)
    d1 <- sum(data$time == t & data$event == 1L & data$arm == "treated")
    observed <- observed + d1
    expected <- expected + d * n1 / n
    if (n > 1) variance <- variance + d * (n1 / n) * (1 - n1 / n) * (n - d) / (n - 1)
  }
  (observed - expected)^2 / variance
}
subjects <- 4000L
data <- list(
  time = (seq_len(subjects) * 7919L) %% 1000L + 1L,
  event = as.integer(seq_len(subjects) %% 3L != 0L),
  arm = ifelse(seq_len(subjects) %% 2L == 0L, "treated", "control")
)
measure("log_rank", function() log_rank(data))
```

### loop_sum

A numeric loop: a scalar double accumulator in a counted loop.

Plain R:

```r
loop_sum <- function(n) {
  s <- 0
  for (i in 1:n) s <- s + i * 0.5
  s
}
measure("loop_sum", function() loop_sum(10000000L))
```

Typed:

```r
loop_sum <- function(n: integer) -> double {
  s <- 0
  for (i in 1:n) s <- s + i * 0.5
  s
}
measure("loop_sum", function() loop_sum(10000000L))
```

### vasicek

A simulation kernel from the finance code base: Vasicek interest-rate paths, one matrix row per step.

Plain R:

```r
vasicek_paths <- function(paths = 200L, steps = 2000L, r0 = 0.03, kappa = 0.5, theta = 0.04, sigma = 0.01, dt = 1 / 12, seed = 9L) {
  set.seed(seed)
  rates <- matrix(r0, nrow = steps + 1L, ncol = paths)
  for (t in seq_len(steps)) {
    shock <- stats::rnorm(paths)
    previous <- rates[t, ]
    rates[t + 1L, ] <- previous + kappa * (theta - previous) * dt + sigma * sqrt(dt) * shock
  }
  rates
}
measure("vasicek", function() sum(vasicek_paths()))
```

Typed:

```r
vasicek_paths <- function(paths: integer = 200L, steps: integer = 2000L, r0: double = 0.03, kappa: double = 0.5, theta: double = 0.04, sigma: double = 0.01, dt: double = 1 / 12, seed: integer = 9L) -> matrix {
  set.seed(seed)
  rates <- matrix(r0, nrow = steps + 1L, ncol = paths)
  for (t in seq_len(steps)) {
    shock <- stats::rnorm(paths)
    previous: double[] <- rates[t, ]
    rates[t + 1L, ] <- previous + kappa * (theta - previous) * dt + sigma * sqrt(dt) * shock
  }
  rates
}
measure("vasicek", function() sum(vasicek_paths()))
```

### vector_arith

Vector code: whole-vector arithmetic on a million doubles, repeated in a loop.

Plain R:

```r
vector_arith <- function(x, times) {
  total <- 0
  for (k in seq_len(times)) {
    y <- x * 2 + x / 3 - 1
    total <- total + sum(y)
  }
  total
}
x <- seq(0, 1, length.out = 1000000L)
measure("vector_arith", function() vector_arith(x, 20L))
```

Typed:

```r
vector_arith <- function(x: double[], times: integer) -> double {
  total <- 0
  for (k in seq_len(times)) {
    y <- x * 2 + x / 3 - 1
    total <- total + sum(y)
  }
  total
}
x: double[] <- seq(0, 1, length.out = 1000000L)
measure("vector_arith", function() vector_arith(x, 20L))
```

### harness

Shared by every benchmark: runs `run()` once to warm up (GNU R's JIT compiles, ry tiers up), then three timed times, and prints the checksum and the best time in seconds.

Plain R:

```r
measure <- function(name, run) {
  value <- run()
  best <- Inf
  for (k in 1:3) {
    started <- proc.time()[[3]]
    again <- run()
    finished <- proc.time()[[3]]
    best <- min(best, finished - started)
    if (!identical(again, value)) stop("the benchmark is not deterministic")
  }
  cat(sprintf("%s %s %.4f\n", name, format(value, digits = 15), best))
}
```

Typed:

```r
measure <- function(name: character, run: fn() -> double) -> NULL {
  value <- run()
  best: double <- Inf
  for (k in 1:3) {
    started: double <- proc.time()[[3]]
    again <- run()
    finished: double <- proc.time()[[3]]
    best <- min(best, finished - started)
    if (!identical(again, value)) stop("the benchmark is not deterministic")
  }
  cat(sprintf("%s %s %.4f\n", name, format(value, digits = 15), best))
}
```

## How it was measured

Each program times its own work with the harness above, so startup and the static check are left
out. The harness runs once to warm up, then three times, and keeps the fastest run.

`just bench-why --write` runs every program twice per engine, alternating engines, and stores the
fastest time and the largest memory peak in `tools/bench/why/results.json`, which every table on
this site reads.

Measured on Intel(R) Xeon(R) Processor @ 2.80GHz, 4 cores, Linux 6.18.44-fc-v185, against GNU R 4.3.3.