---
title: "Benchmark with a reproducible receipt"
name: benchmarks
description: Run a scoped micro-comparison and report enough context to reproduce it.
output:
  rmarkdown::html_vignette:
    toc: yes
    toc_depth: 2
    css: albers.css
    includes:
      in_header: albers-header.html
params:
  family: red
  base_size: 13
  content_width: 72
  preset: homage
resource_files:
- albers.css
- albers.js
- albers-header.html
- fonts
vignette: |
  %\VignetteIndexEntry{Benchmark with a reproducible receipt}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE, comment = "#>", message = FALSE, warning = FALSE,
  fig.align = "center", fig.retina = 2, out.width = "100%",
  fig.width = 7, fig.asp = 0.56
)
if (requireNamespace("albersdown", quietly = TRUE) &&
    requireNamespace("ggplot2", quietly = TRUE)) {
  ggplot2::theme_set(albersdown::theme_albers(family = params$family, preset = params$preset))
}
library(neurocluster)
```

```{r albers-classes, echo=FALSE, results='asis'}
cat(sprintf(
  paste0('<script>document.addEventListener("DOMContentLoaded",function(){',
         'document.body.classList.add("palette-%s","preset-%s");',
         '});</script><style>:root{--content:%sch}</style>'),
  params$family, params$preset, params$content_width
))
```

A benchmark is evidence only for its measured workload. This article does not
read a checked-in `results.csv` or present old machine timings as current. It
runs a deliberately small comparison during rendering and shows the receipt
needed to reproduce it. Use the pattern—not the tiny timings—to design a study
for your data.

## What exactly is measured here?

The estimand is elapsed wall-clock seconds for a complete `cluster4d()` call on
one generated 400-voxel, 18-time-point volume, after one warm-up call. We keep
three repetitions for each method and require every fit to be valid and to
produce the requested K. These are smoke-scale timings, not whole-brain
performance claims.

```{r benchmark-data, include=FALSE}
syn <- generate_synthetic_volume(
  scenario = "gaussian_blobs", dims = c(10, 10, 4),
  n_clusters = 4, n_time = 18, noise_sd = 0.06, seed = 713
)
fit_one <- function(method) {
  if (method == "snic") {
    cluster4d(syn$vec, syn$mask, 4, method = "snic", spatial_weight = 0.4)
  } else {
    cluster4d(
      syn$vec, syn$mask, 4, method = "slic", spatial_weight = 0.4,
      max_iterations = 6, preserve_k = TRUE, parallel = FALSE
    )
  }
}
invisible(fit_one("snic"))
invisible(fit_one("slic"))
```

```{r benchmark-run}
methods <- c("snic", "slic")
benchmark_order <- c("snic", "slic", "slic", "snic", "snic", "slic")
repetitions <- as.integer(ave(
  seq_along(benchmark_order), benchmark_order, FUN = seq_along
))
data.frame(run = seq_along(benchmark_order), method = benchmark_order,
           repetition = repetitions)

timings <- do.call(rbind, lapply(seq_along(benchmark_order), function(run) {
  method <- benchmark_order[run]
  elapsed <- system.time(fit <- fit_one(method))[["elapsed"]]
  validation <- validate_cluster4d(fit, syn$vec, syn$mask)
  data.frame(
    run = run, method = method, repetition = repetitions[run],
    elapsed_seconds = unname(elapsed),
    actual_k = fit$actual_k, valid = validation$valid
  )
}))
timings
```

```{r benchmark-check, include=FALSE}
stopifnot(nrow(timings) == 6L)
stopifnot(all(timings$valid), all(timings$actual_k == 4L))
stopifnot(all(is.finite(timings$elapsed_seconds)))
stopifnot(all(timings$elapsed_seconds >= 0))
stopifnot(identical(timings$method, benchmark_order))
stopifnot(all(table(timings$method) == 3L))
benchmark_summary <- aggregate(elapsed_seconds ~ method, timings, median)
benchmark_summary <- benchmark_summary[match(methods, benchmark_summary$method), ]
```

```{r benchmark-figure, echo=FALSE, fig.cap='Three elapsed-time repetitions per method on the vignette’s tiny generated workload; diamonds are medians.', fig.alt='Jittered timing points and median diamonds for SNIC and serial SLIC on a 400-voxel example.'}
palette_cb <- c(snic = "#0072B2", slic = "#D55E00")
x <- match(timings$method, methods)
offset <- (timings$repetition - 2) * 0.08
plot(x + offset, timings$elapsed_seconds, pch = 21, cex = 1.5,
     bg = palette_cb[timings$method], col = "white",
     xaxt = "n", xlab = "method", ylab = "elapsed seconds",
     xlim = c(0.5, 2.5), ylim = c(0, max(timings$elapsed_seconds) * 1.15 + 1e-6))
axis(1, at = seq_along(methods), labels = toupper(methods))
points(seq_along(methods), benchmark_summary$elapsed_seconds,
       pch = 23, cex = 1.8, bg = "#009E73", col = "white")
box()
```

**What to notice.** Repeated timings vary even on a tiny fixed input. The plot
supports only this render’s workload and environment; it does not rank methods
for a different volume, K, parameter set, or machine.

## What belongs in the receipt?

```{r benchmark-receipt}
list(
  package_version = as.character(packageVersion("neurocluster")),
  r_version = paste(R.version$major, R.version$minor, sep = "."),
  platform = R.version$platform,
  dimensions = syn$dims,
  included_voxels = length(syn$truth),
  time_points = ncol(syn$patterns),
  requested_k = 4L,
  repetitions = 3L,
  run_order = benchmark_order
)
```

For a real benchmark, also record the commit, CPU, RAM, operating system,
thread backend and count, full method calls, warm-up policy, allocation or peak
RSS measurement method, and every repetition. Predefine quality estimands when
comparing approximation settings; elapsed time alone cannot establish a useful
clustering. See [Performance and memory](performance-memory.html) and
[Validate and compare clusterings](validate-compare.html).