bench: 2026-07 public benchmark refresh harness with fail-closed publication gates

.gitignore

@@ -111,3 +111,9 @@ survey-did-paper-arxiv-v2/
 # Reviewer-eval harness — local run artifacts, detached worktrees, comparison bundles
 tools/reviewer-eval/runs/
 
+
+# Benchmark refresh (2026-07) - venvs and raw per-rep artifacts are local-only;
+# consolidated results JSON and the version-story report ARE committed.
+benchmarks/refresh_2026_07/venvs/
+benchmarks/refresh_2026_07/results/raw/
+benchmarks/refresh_2026_07/results/refresh_results_smoke.json

CHANGELOG.md

@@ -7,6 +7,43 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [3.7.0] - 2026-07-08
 
+### Added
+- **Public benchmark refresh harness (`benchmarks/refresh_2026_07/`).** Fair,
+  fail-closed re-run machinery for the published docs benchmarks against the
+  RELEASED diff-diff wheel (isolated `uv` venv; provenance hard-fails on
+  dev-tree shadowing, wrong backend, or unexpected `DIFF_DIFF_*` knobs) vs R
+  fixest/did/synthdid: untimed in-process warm-up on both sides (R's
+  byte-compiler JIT previously sat inside the timing window), fresh
+  subprocess per replication, strictly sequential, median of counted reps,
+  thread-count env vars stripped and R run under `--vanilla` so package
+  defaults are enforced rather than assumed. Publication gates cover ATT/SE
+  tolerances per rendered arm (CI overlap enforced independently), per-period
+  / per-(g,t) / event-study / group effect surfaces, SDID id-aligned
+  unit/time weight identity (1e-8), MPDTA known-answer, replication
+  determinism and non-finite outputs, and environment fingerprints;
+  `gen_benchmark_tables.py` refuses to render any hard-gated or
+  mixed-environment artifact into `docs/benchmarks.rst` (marker-bounded
+  regions). SDID placebo-SE parity is gated at a documented 35% Monte
+  Carlo bound (R's placebo permutation is unseeded; see the SyntheticDiD
+  registry note). Benchmark scripts (both languages) gained `--warmup`,
+  strict flag parsing, and estimator-field self-identification so a
+  mislabeled cell can never publish; the legacy `--type twfe` no-op on the
+  BasicDiD scripts now fails loudly (absorbed-FE TWFE has its own dedicated
+  script pair, and the refresh benchmarks it separately from BasicDiD).
+  Regenerated page numbers land together with the committed results
+  artifact once the timed refresh run completes.
+
+### Fixed
+- **`docs/benchmarks.rst` SDID weight-parity claim now machine-verified (and
+  its earlier wording corrected).** The page's "identical unit and time
+  weights" claim was previously unenforced; the new id-aligned benchmark
+  gate verifies it at < 1e-8 (observed ≤ 5e-13 vs R). Note for anyone
+  comparing weights manually: `SyntheticDiDResults.get_unit_weights_df()`
+  returns weights sorted by descending weight while R's `synthdid` uses
+  panel control-row order — a naive positional comparison fabricates
+  apparent ~1e-2 divergence out of pure ordering (documented in the
+  SyntheticDiD registry note).
+
 ### Removed
 - **HAD pretest helpers: deprecated `survey=` / `weights=` kwargs removed (3.7.x).**
   Completes the HAD survey-design API consolidation started with

TODO.md

@@ -45,6 +45,10 @@ generic sparse-FE, QR+SVD rank-detection redundancy, `check_finite` bypass — m
 
 ### Testing / docs
 
+| Issue | Location | Origin | Effort | Priority |
+|-------|----------|--------|--------|----------|
+| Benchmark refresh phase 2 (same PR #672): run the gated timed refresh on an idle machine, commit `benchmarks/refresh_2026_07/results/refresh_results.json`, regenerate the marker-bounded regions of `docs/benchmarks.rst` via `gen_benchmark_tables.py`, and reconcile the remaining pre-refresh prose (protocol bullets "3 replications / mean ± std", combined BasicDiD/TWFE wording, SDID note under the perf table, Key Observations, "Reproducing Benchmarks" section, `llms.txt` speedup cross-references). Row removed by the phase-2 push itself. | `benchmarks/refresh_2026_07/`, `docs/benchmarks.rst` | #672 | Mid | Medium |
+
 | Issue | Location | Origin | Effort | Priority |
 |-------|----------|--------|--------|----------|
 

benchmarks/R/benchmark_did.R

@@ -40,7 +40,8 @@ parse_args <- function(args) {
     cband = FALSE,
     pl = FALSE,
     cores = 1L,
-    faster_mode = NULL  # NULL -> use did's own default for this version
+    faster_mode = NULL,  # NULL -> use did's own default for this version
+    warmup = FALSE
   )
 
   i <- 1
@@ -72,6 +73,8 @@ parse_args <- function(args) {
       result$cores <- as.integer(val)
     } else if (flag == "--faster-mode") {
       result$faster_mode <- parse_bool(val, flag)
+    } else if (flag == "--warmup") {
+      result$warmup <- parse_bool(val, flag)
     } else {
       # Unknown flags used to be silently skipped, which turned typos into
       # silent default runs. Fail loudly instead.
@@ -110,6 +113,19 @@ data[, first_treat := as.numeric(first_treat)]
 data[first_treat == 0, first_treat := Inf]
 message(sprintf("Never-treated units (first_treat=Inf): %d", sum(is.infinite(data$first_treat))))
 
+# Warm-up: run the FULL timed pipeline (att_gt + the aggte calls) once
+# untimed so byte-compiler JIT and first-call setup stay out of the window.
+run_warmup <- function(att_gt_args, config) {
+  message("Warm-up fit (untimed)...")
+  out_w <- do.call(att_gt, att_gt_args)
+  invisible(aggte(out_w, type = "simple", bstrap = config$bstrap,
+                  biters = config$biters, cband = FALSE))
+  invisible(aggte(out_w, type = "dynamic", bstrap = config$bstrap,
+                  biters = config$biters, cband = config$cband))
+  invisible(aggte(out_w, type = "group", bstrap = config$bstrap,
+                  biters = config$biters, cband = FALSE))
+}
+
 # Run benchmark
 message("Running Callaway-Sant'Anna estimation...")
 start_time <- Sys.time()
@@ -143,6 +159,11 @@ if (!is.null(config$faster_mode)) {
     effective_faster_mode <- "unsupported-ignored"
   }
 }
+if (config$warmup) {
+  run_warmup(att_gt_args, config)
+  # Reset the clock: the warm-up must stay out of the timing window.
+  start_time <- Sys.time()
+}
 out <- do.call(att_gt, att_gt_args)
 
 estimation_time <- as.numeric(difftime(Sys.time(), start_time, units = "secs"))
@@ -216,6 +237,8 @@ results <- list(
     pl = config$pl,
     cores = config$cores,
     faster_mode = effective_faster_mode,
+    dt_threads = data.table::getDTthreads(),
+    warmup = config$warmup,
     blas = tryCatch(sessionInfo()$BLAS, error = function(e) NULL)
   )
 )

benchmarks/R/benchmark_fixest.R

@@ -11,12 +11,20 @@ library(data.table)
 # Parse command line arguments
 args <- commandArgs(trailingOnly = TRUE)
 
+parse_bool <- function(x, flag) {
+  v <- tolower(x)
+  if (v %in% c("true", "t", "1", "yes")) return(TRUE)
+  if (v %in% c("false", "f", "0", "no")) return(FALSE)
+  stop(sprintf("Invalid boolean for %s: '%s' (use true/false)", flag, x))
+}
+
 parse_args <- function(args) {
   result <- list(
     data = NULL,
     output = NULL,
     cluster = "unit",
-    type = "twfe"  # "basic" for simple 2x2, "twfe" for two-way FE
+    type = "basic",  # only "basic" runs here; benchmark_twfe.R does TWFE
+    warmup = FALSE
   )
 
   i <- 1
@@ -33,24 +41,45 @@ parse_args <- function(args) {
     } else if (args[i] == "--type") {
       result$type <- args[i + 1]
       i <- i + 2
+    } else if (args[i] == "--warmup") {
+      result$warmup <- parse_bool(args[i + 1], "--warmup")
+      i <- i + 2
     } else {
-      i <- i + 1
+      # Unknown flags used to be silently skipped, which turned typos into
+      # silent protocol changes. Fail loudly instead.
+      stop(sprintf("Unknown flag: %s", args[i]))
     }
   }
 
   if (is.null(result$data) || is.null(result$output)) {
-    stop("Usage: Rscript benchmark_fixest.R --data <path> --output <path> [--cluster unit|time] [--type basic|twfe]")
+    stop("Usage: Rscript benchmark_fixest.R --data <path> --output <path> [--cluster unit|time] [--type basic] (absorbed-FE TWFE lives in benchmark_twfe.R)")
   }
 
   return(result)
 }
 
 config <- parse_args(args)
 
+if (config$type == "twfe") {
+  stop(paste(
+    "--type twfe is not implemented by benchmark_fixest.R (it always runs",
+    "the interaction OLS). Use benchmark_twfe.R for absorbed-FE TWFE."
+  ))
+} else if (config$type != "basic") {
+  stop(sprintf("Invalid --type '%s' (only 'basic' is supported here)", config$type))
+}
+
 # Load data
 message(sprintf("Loading data from: %s", config$data))
 data <- fread(config$data)
 
+# Warm-up: run the full estimation once untimed so byte-compiler JIT and
+# first-call setup costs stay out of the timing window.
+if (config$warmup) {
+  message("Warm-up fit (untimed)...")
+  invisible(feols(outcome ~ treated * post, data = data, cluster = config$cluster))
+}
+
 # Run benchmark
 message(sprintf("Running %s DiD estimation...", toupper(config$type)))
 start_time <- Sys.time()
@@ -138,6 +167,9 @@ results <- list(
   metadata = list(
     r_version = R.version.string,
     fixest_version = as.character(packageVersion("fixest")),
+    nthreads = fixest::getFixest_nthreads(),
+    dt_threads = data.table::getDTthreads(),
+    warmup = config$warmup,
     n_units = length(unique(data$unit)),
     n_periods = length(unique(data$time)),
     n_obs = nrow(data)

benchmarks/R/benchmark_multiperiod.R

@@ -12,14 +12,22 @@ library(data.table)
 # Parse command line arguments
 args <- commandArgs(trailingOnly = TRUE)
 
+parse_bool <- function(x, flag) {
+  v <- tolower(x)
+  if (v %in% c("true", "t", "1", "yes")) return(TRUE)
+  if (v %in% c("false", "f", "0", "no")) return(FALSE)
+  stop(sprintf("Invalid boolean for %s: '%s' (use true/false)", flag, x))
+}
+
 parse_args <- function(args) {
   result <- list(
     data = NULL,
     output = NULL,
     cluster = "unit",
     n_pre = NULL,
     n_post = NULL,
-    reference_period = NULL
+    reference_period = NULL,
+    warmup = FALSE
   )
 
   i <- 1
@@ -42,8 +50,13 @@ parse_args <- function(args) {
     } else if (args[i] == "--reference-period") {
       result$reference_period <- as.integer(args[i + 1])
       i <- i + 2
+    } else if (args[i] == "--warmup") {
+      result$warmup <- parse_bool(args[i + 1], "--warmup")
+      i <- i + 2
     } else {
-      i <- i + 1
+      # Unknown flags used to be silently skipped, which turned typos into
+      # silent protocol changes. Fail loudly instead.
+      stop(sprintf("Unknown flag: %s", args[i]))
     }
   }
 
@@ -75,6 +88,15 @@ message(sprintf("n_pre: %d, n_post: %d", config$n_pre, config$n_post))
 # Create factor for time with reference level
 data[, time_f := relevel(factor(time), ref = as.character(ref_period))]
 
+# Warm-up: run the full estimation once untimed so byte-compiler JIT and
+# first-call setup costs stay out of the timing window.
+cluster_formula_warm <- as.formula(paste0("~", config$cluster))
+if (config$warmup) {
+  message("Warm-up fit (untimed)...")
+  invisible(feols(outcome ~ treated * time_f | unit, data = data,
+                  cluster = cluster_formula_warm))
+}
+
 # Run benchmark
 message("Running MultiPeriodDiD estimation (fixest::feols)...")
 start_time <- Sys.time()
@@ -186,6 +208,9 @@ results <- list(
   metadata = list(
     r_version = R.version.string,
     fixest_version = as.character(packageVersion("fixest")),
+    nthreads = fixest::getFixest_nthreads(),
+    dt_threads = data.table::getDTthreads(),
+    warmup = config$warmup,
     n_units = length(unique(data$unit)),
     n_periods = length(unique(data$time)),
     n_obs = nrow(data),

benchmarks/R/benchmark_synthdid.R

@@ -11,10 +11,19 @@ library(data.table)
 # Parse command line arguments
 args <- commandArgs(trailingOnly = TRUE)
 
+parse_bool <- function(x, flag) {
+  v <- tolower(x)
+  if (v %in% c("true", "t", "1", "yes")) return(TRUE)
+  if (v %in% c("false", "f", "0", "no")) return(FALSE)
+  stop(sprintf("Invalid boolean for %s: '%s' (use true/false)", flag, x))
+}
+
 parse_args <- function(args) {
   result <- list(
     data = NULL,
-    output = NULL
+    output = NULL,
+    warmup = FALSE,
+    skip_jackknife = FALSE
   )
 
   i <- 1
@@ -25,8 +34,16 @@ parse_args <- function(args) {
     } else if (args[i] == "--output") {
       result$output <- args[i + 1]
       i <- i + 2
+    } else if (args[i] == "--warmup") {
+      result$warmup <- parse_bool(args[i + 1], "--warmup")
+      i <- i + 2
+    } else if (args[i] == "--skip-jackknife") {
+      result$skip_jackknife <- parse_bool(args[i + 1], "--skip-jackknife")
+      i <- i + 2
     } else {
-      i <- i + 1
+      # Unknown flags used to be silently skipped, which turned typos into
+      # silent protocol changes. Fail loudly instead.
+      stop(sprintf("Unknown flag: %s", args[i]))
     }
   }
 
@@ -60,6 +77,15 @@ setup <- panel.matrices(
   treatment = "treatment"
 )
 
+# Warm-up: run the FULL timed pipeline (estimation + placebo vcov) once
+# untimed so byte-compiler JIT and first-call setup stay out of the window.
+if (config$warmup) {
+  message("Warm-up fit (untimed)...")
+  tau_warm <- synthdid_estimate(setup$Y, setup$N0, setup$T0)
+  invisible(vcov(tau_warm, method = "placebo"))
+  rm(tau_warm)
+}
+
 # Run benchmark
 message("Running Synthetic DiD estimation...")
 start_time <- Sys.time()
@@ -80,12 +106,20 @@ se_matrix <- vcov(tau_hat, method = "placebo")
 se <- as.numeric(sqrt(se_matrix[1, 1]))  # Extract scalar SE
 se_time <- as.numeric(difftime(Sys.time(), se_start, units = "secs"))
 
-# Compute SE via jackknife (Algorithm 3)
-message("Computing jackknife standard errors...")
-se_jk_start <- Sys.time()
-se_jk_matrix <- vcov(tau_hat, method = "jackknife")
-se_jackknife <- as.numeric(sqrt(se_jk_matrix[1, 1]))
-se_jk_time <- as.numeric(difftime(Sys.time(), se_jk_start, units = "secs"))
+# Compute SE via jackknife (Algorithm 3). Not part of total_seconds (which
+# is estimation + placebo SE, matching the Python arm); skippable because it
+# wastes hours of wall-clock at large scales in timing-focused runs.
+if (config$skip_jackknife) {
+  message("Skipping jackknife standard errors (--skip-jackknife)")
+  se_jackknife <- NA_real_
+  se_jk_time <- NA_real_
+} else {
+  message("Computing jackknife standard errors...")
+  se_jk_start <- Sys.time()
+  se_jk_matrix <- vcov(tau_hat, method = "jackknife")
+  se_jackknife <- as.numeric(sqrt(se_jk_matrix[1, 1]))
+  se_jk_time <- as.numeric(difftime(Sys.time(), se_jk_start, units = "secs"))
+}
 
 total_time <- estimation_time + se_time  # placebo only, matches `se` field
 
@@ -107,9 +141,12 @@ results <- list(
   se = se,
   se_jackknife = se_jackknife,
 
-  # Weights (full precision)
+  # Weights (full precision; omega is in panel control-row order, lambda in
+  # pre-period column order - emit the ids so comparison aligns by key)
   unit_weights = as.numeric(unit_weights),
+  unit_weight_ids = rownames(setup$Y)[seq_len(N0)],
   time_weights = as.numeric(time_weights),
+  time_weight_ids = colnames(setup$Y)[seq_len(T0)],
 
   # Regularization parameters
   noise_level = noise_level,
@@ -128,6 +165,9 @@ results <- list(
   metadata = list(
     r_version = R.version.string,
     synthdid_version = as.character(packageVersion("synthdid")),
+    dt_threads = data.table::getDTthreads(),
+    warmup = config$warmup,
+    skip_jackknife = config$skip_jackknife,
     n_control = N0,
     n_treated = N1,
     n_pre_periods = T0,