diff --git a/.gitignore b/.gitignore index 239debff..e4c627d7 100644 --- a/.gitignore +++ b/.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 diff --git a/CHANGELOG.md b/CHANGELOG.md index d00f6239..9d6f1ed3 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/TODO.md b/TODO.md index 38cdcc82..a6872c79 100644 --- a/TODO.md +++ b/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 | |-------|----------|--------|--------|----------| diff --git a/benchmarks/R/benchmark_did.R b/benchmarks/R/benchmark_did.R index 210d98b2..4579c8aa 100644 --- a/benchmarks/R/benchmark_did.R +++ b/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) ) ) diff --git a/benchmarks/R/benchmark_fixest.R b/benchmarks/R/benchmark_fixest.R index dd055a44..0d80a94f 100644 --- a/benchmarks/R/benchmark_fixest.R +++ b/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,13 +41,18 @@ 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 --output [--cluster unit|time] [--type basic|twfe]") + stop("Usage: Rscript benchmark_fixest.R --data --output [--cluster unit|time] [--type basic] (absorbed-FE TWFE lives in benchmark_twfe.R)") } return(result) @@ -47,10 +60,26 @@ parse_args <- function(args) { 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) diff --git a/benchmarks/R/benchmark_multiperiod.R b/benchmarks/R/benchmark_multiperiod.R index 71324883..9dbb62a8 100644 --- a/benchmarks/R/benchmark_multiperiod.R +++ b/benchmarks/R/benchmark_multiperiod.R @@ -12,6 +12,13 @@ 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, @@ -19,7 +26,8 @@ parse_args <- function(args) { 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), diff --git a/benchmarks/R/benchmark_synthdid.R b/benchmarks/R/benchmark_synthdid.R index 496be3a7..bd18bc3d 100644 --- a/benchmarks/R/benchmark_synthdid.R +++ b/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, diff --git a/benchmarks/R/benchmark_twfe.R b/benchmarks/R/benchmark_twfe.R index e1c9771a..7a0e3c8f 100644 --- a/benchmarks/R/benchmark_twfe.R +++ b/benchmarks/R/benchmark_twfe.R @@ -14,10 +14,18 @@ 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 ) i <- 1 @@ -28,8 +36,13 @@ 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 { - 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])) } } @@ -46,6 +59,13 @@ config <- parse_args(args) 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 | unit + post, data = data, cluster = ~unit)) +} + # Run benchmark message("Running TWFE estimation with absorbed FE...") start_time <- Sys.time() @@ -119,6 +139,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$post)), n_obs = nrow(data) diff --git a/benchmarks/python/benchmark_basic.py b/benchmarks/python/benchmark_basic.py index e173d6f0..12a010c2 100644 --- a/benchmarks/python/benchmark_basic.py +++ b/benchmarks/python/benchmark_basic.py @@ -29,11 +29,18 @@ def _get_backend_from_args(): import numpy as np import pandas as pd -# Add parent to path for imports -sys.path.insert(0, str(Path(__file__).parent.parent.parent)) +# Add repo root to path for benchmarks.python.utils. When +# DIFF_DIFF_BENCH_USE_INSTALLED=1 (isolated-venv refresh runs), APPEND so the +# venv's installed diff-diff wheel wins over the dev tree; otherwise PREPEND +# (historical behavior: benchmark the working tree). +_REPO_ROOT = str(Path(__file__).parent.parent.parent) +if os.environ.get("DIFF_DIFF_BENCH_USE_INSTALLED") == "1": + sys.path.append(_REPO_ROOT) +else: + sys.path.insert(0, _REPO_ROOT) from diff_diff import DifferenceInDifferences, HAS_RUST_BACKEND -from benchmarks.python.utils import Timer +from benchmarks.python.utils import Timer, collect_provenance def parse_args(): @@ -44,13 +51,18 @@ def parse_args(): "--cluster", default="unit", help="Column to cluster standard errors on" ) parser.add_argument( - "--type", default="twfe", choices=["basic", "twfe"], - help="Estimator type (basic or twfe, default: twfe)" + "--type", default="basic", choices=["basic", "twfe"], + help="Accepted for backward compatibility; only 'basic' runs here " + "(use benchmark_twfe.py for absorbed-FE TWFE)" ) parser.add_argument( "--backend", default="auto", choices=["auto", "python", "rust"], help="Backend to use: auto (default), python (pure Python), rust (Rust backend)" ) + parser.add_argument( + "--warmup", action="store_true", + help="Run one untimed fit before the timed fit (JIT/cache warm-up)" + ) return parser.parse_args() @@ -70,9 +82,22 @@ def main(): print(f"Loading data from: {args.data}") data = pd.read_csv(args.data) + if args.type == "twfe": + raise SystemExit( + "--type twfe is not implemented by benchmark_basic.py (it always " + "runs the interaction OLS). Use benchmark_twfe.py for the " + "absorbed-FE TwoWayFixedEffects estimator." + ) + # Run benchmark print("Running DiD estimation...") + if args.warmup: + print("Warm-up fit (untimed)...") + DifferenceInDifferences(robust=True, cluster=args.cluster).fit( + data, formula="outcome ~ treated * post" + ) + # Use DifferenceInDifferences with formula to match R's fixest::feols did = DifferenceInDifferences(robust=True, cluster=args.cluster) @@ -116,7 +141,10 @@ def main(): "n_units": len(data["unit"].unique()), "n_periods": len(data["time"].unique()), "n_obs": len(data), + "warmup": args.warmup, }, + # Wheel/backend provenance (refresh runs hard-fail on mismatch) + "provenance": collect_provenance(), } # Write output diff --git a/benchmarks/python/benchmark_callaway.py b/benchmarks/python/benchmark_callaway.py index ba999086..5fb959d8 100644 --- a/benchmarks/python/benchmark_callaway.py +++ b/benchmarks/python/benchmark_callaway.py @@ -29,11 +29,18 @@ def _get_backend_from_args(): import numpy as np import pandas as pd -# Add parent to path for imports -sys.path.insert(0, str(Path(__file__).parent.parent.parent)) +# Add repo root to path for benchmarks.python.utils. When +# DIFF_DIFF_BENCH_USE_INSTALLED=1 (isolated-venv refresh runs), APPEND so the +# venv's installed diff-diff wheel wins over the dev tree; otherwise PREPEND +# (historical behavior: benchmark the working tree). +_REPO_ROOT = str(Path(__file__).parent.parent.parent) +if os.environ.get("DIFF_DIFF_BENCH_USE_INSTALLED") == "1": + sys.path.append(_REPO_ROOT) +else: + sys.path.insert(0, _REPO_ROOT) from diff_diff import CallawaySantAnna, HAS_RUST_BACKEND -from benchmarks.python.utils import BenchmarkResult, Timer +from benchmarks.python.utils import BenchmarkResult, Timer, collect_provenance def parse_args(): @@ -58,6 +65,10 @@ def parse_args(): "--backend", default="auto", choices=["auto", "python", "rust"], help="Backend to use: auto (default), python (pure Python), rust (Rust backend)" ) + parser.add_argument( + "--warmup", action="store_true", + help="Run one untimed fit before the timed fit (JIT/cache warm-up)" + ) return parser.parse_args() @@ -81,6 +92,22 @@ def main(): print("Running Callaway-Sant'Anna estimation...") # Use analytical SE (n_bootstrap=0) - matches R's did package after # influence function aggregation fix (accounts for covariance) + if args.warmup: + print("Warm-up fit (untimed)...") + CallawaySantAnna( + estimation_method=args.method, + control_group=args.control_group, + n_bootstrap=0, + seed=42, + ).fit( + df, + outcome="outcome", + time="time", + unit="unit", + first_treat="first_treat", + aggregate="all", + ) + cs = CallawaySantAnna( estimation_method=args.method, control_group=args.control_group, @@ -163,7 +190,10 @@ def main(): "n_units": n_units, "n_periods": n_periods, "n_obs": n_obs, + "warmup": args.warmup, }, + # Wheel/backend provenance (refresh runs hard-fail on mismatch) + "provenance": collect_provenance(), } # Write output diff --git a/benchmarks/python/benchmark_multiperiod.py b/benchmarks/python/benchmark_multiperiod.py index 5507aa47..2122f242 100644 --- a/benchmarks/python/benchmark_multiperiod.py +++ b/benchmarks/python/benchmark_multiperiod.py @@ -29,11 +29,18 @@ def _get_backend_from_args(): # NOW import diff_diff and other dependencies (will see the env var) import pandas as pd -# Add parent to path for imports -sys.path.insert(0, str(Path(__file__).parent.parent.parent)) +# Add repo root to path for benchmarks.python.utils. When +# DIFF_DIFF_BENCH_USE_INSTALLED=1 (isolated-venv refresh runs), APPEND so the +# venv's installed diff-diff wheel wins over the dev tree; otherwise PREPEND +# (historical behavior: benchmark the working tree). +_REPO_ROOT = str(Path(__file__).parent.parent.parent) +if os.environ.get("DIFF_DIFF_BENCH_USE_INSTALLED") == "1": + sys.path.append(_REPO_ROOT) +else: + sys.path.insert(0, _REPO_ROOT) from diff_diff import MultiPeriodDiD, HAS_RUST_BACKEND -from benchmarks.python.utils import Timer +from benchmarks.python.utils import Timer, collect_provenance def parse_args(): @@ -57,6 +64,10 @@ def parse_args(): "--backend", default="auto", choices=["auto", "python", "rust"], help="Backend to use: auto (default), python (pure Python), rust (Rust backend)" ) + parser.add_argument( + "--warmup", action="store_true", + help="Run one untimed fit before the timed fit (JIT/cache warm-up)" + ) return parser.parse_args() @@ -89,6 +100,18 @@ def main(): # Run benchmark print("Running MultiPeriodDiD estimation...") + if args.warmup: + print("Warm-up fit (untimed)...") + MultiPeriodDiD(robust=True, cluster=args.cluster).fit( + data, + outcome="outcome", + treatment="treated", + time="time", + post_periods=post_periods, + reference_period=ref_period, + absorb=["unit"], + ) + did = MultiPeriodDiD(robust=True, cluster=args.cluster) with Timer() as timer: @@ -142,7 +165,10 @@ def main(): "n_obs": len(data), "n_pre": n_pre, "n_post": len(post_periods), + "warmup": args.warmup, }, + # Wheel/backend provenance (refresh runs hard-fail on mismatch) + "provenance": collect_provenance(), } # Write output diff --git a/benchmarks/python/benchmark_synthdid.py b/benchmarks/python/benchmark_synthdid.py index 4392bf45..53f6e01f 100644 --- a/benchmarks/python/benchmark_synthdid.py +++ b/benchmarks/python/benchmark_synthdid.py @@ -12,6 +12,7 @@ import sys from pathlib import Path + # IMPORTANT: Parse --backend and set environment variable BEFORE importing diff_diff # This ensures the backend configuration is respected by all modules def _get_backend_from_args(): @@ -21,6 +22,7 @@ def _get_backend_from_args(): args, _ = parser.parse_known_args() return args.backend + _requested_backend = _get_backend_from_args() if _requested_backend in ("python", "rust"): os.environ["DIFF_DIFF_BACKEND"] = _requested_backend @@ -29,11 +31,18 @@ def _get_backend_from_args(): import numpy as np import pandas as pd -# Add parent to path for imports -sys.path.insert(0, str(Path(__file__).parent.parent.parent)) +# Add repo root to path for benchmarks.python.utils. When +# DIFF_DIFF_BENCH_USE_INSTALLED=1 (isolated-venv refresh runs), APPEND so the +# venv's installed diff-diff wheel wins over the dev tree; otherwise PREPEND +# (historical behavior: benchmark the working tree). +_REPO_ROOT = str(Path(__file__).parent.parent.parent) +if os.environ.get("DIFF_DIFF_BENCH_USE_INSTALLED") == "1": + sys.path.append(_REPO_ROOT) +else: + sys.path.insert(0, _REPO_ROOT) from diff_diff import SyntheticDiD, HAS_RUST_BACKEND -from benchmarks.python.utils import Timer +from benchmarks.python.utils import Timer, collect_provenance def parse_args(): @@ -44,13 +53,22 @@ def parse_args(): "--n-bootstrap", type=int, default=200, help="Number of bootstrap iterations" ) parser.add_argument( - "--variance-method", type=str, default="placebo", + "--variance-method", + type=str, + default="placebo", choices=["bootstrap", "jackknife", "placebo"], - help="Variance estimation method (default: placebo to match R)" + help="Variance estimation method (default: placebo to match R)", + ) + parser.add_argument( + "--backend", + default="auto", + choices=["auto", "python", "rust"], + help="Backend to use: auto (default), python (pure Python), rust (Rust backend)", ) parser.add_argument( - "--backend", default="auto", choices=["auto", "python", "rust"], - help="Backend to use: auto (default), python (pure Python), rust (Rust backend)" + "--warmup", + action="store_true", + help="Run one untimed fit before the timed fit (JIT/cache warm-up)", ) return parser.parse_args() @@ -77,11 +95,23 @@ def main(): # Run benchmark print("Running Synthetic DiD estimation...") - sdid = SyntheticDiD( - n_bootstrap=args.n_bootstrap, - variance_method=args.variance_method, - seed=42 - ) + if args.warmup: + # Full pipeline warm-up: estimation + variance are both inside fit() + print("Warm-up fit (untimed)...") + SyntheticDiD( + n_bootstrap=args.n_bootstrap, + variance_method=args.variance_method, + seed=42, + ).fit( + data, + outcome="outcome", + treatment="treated", + unit="unit", + time="time", + post_periods=post_periods, + ) + + sdid = SyntheticDiD(n_bootstrap=args.n_bootstrap, variance_method=args.variance_method, seed=42) with Timer() as timer: results = sdid.fit( @@ -95,9 +125,11 @@ def main(): total_time = timer.elapsed - # Get weights - unit_weights_df = results.get_unit_weights_df() - time_weights_df = results.get_time_weights_df() + # Get weights. NOTE: get_unit_weights_df() returns weights sorted by + # DESCENDING WEIGHT, not panel order - sort by unit/period id and emit + # the ids so cross-implementation comparison can align by key. + unit_weights_df = results.get_unit_weights_df().sort_values("unit") + time_weights_df = results.get_time_weights_df().sort_values("period") # Build output output = { @@ -106,9 +138,11 @@ def main(): # Point estimate and SE "att": float(results.att), "se": float(results.se), - # Weights (full precision) + # Weights (full precision, ordered by ascending unit/period id) "unit_weights": unit_weights_df["weight"].tolist(), + "unit_weight_ids": unit_weights_df["unit"].tolist(), "time_weights": time_weights_df["weight"].tolist(), + "time_weight_ids": time_weights_df["period"].tolist(), # Regularization parameters "noise_level": float(results.noise_level) if results.noise_level is not None else None, "zeta_omega": float(results.zeta_omega) if results.zeta_omega is not None else None, @@ -126,7 +160,10 @@ def main(): "n_post_periods": len(post_periods), "n_bootstrap": args.n_bootstrap, "variance_method": args.variance_method, + "warmup": args.warmup, }, + # Wheel/backend provenance (refresh runs hard-fail on mismatch) + "provenance": collect_provenance(), } # Write output diff --git a/benchmarks/python/benchmark_twfe.py b/benchmarks/python/benchmark_twfe.py index 22cadacf..ce72266b 100644 --- a/benchmarks/python/benchmark_twfe.py +++ b/benchmarks/python/benchmark_twfe.py @@ -30,11 +30,18 @@ def _get_backend_from_args(): # NOW import diff_diff and other dependencies (will see the env var) import pandas as pd -# Add parent to path for imports -sys.path.insert(0, str(Path(__file__).parent.parent.parent)) +# Add repo root to path for benchmarks.python.utils. When +# DIFF_DIFF_BENCH_USE_INSTALLED=1 (isolated-venv refresh runs), APPEND so the +# venv's installed diff-diff wheel wins over the dev tree; otherwise PREPEND +# (historical behavior: benchmark the working tree). +_REPO_ROOT = str(Path(__file__).parent.parent.parent) +if os.environ.get("DIFF_DIFF_BENCH_USE_INSTALLED") == "1": + sys.path.append(_REPO_ROOT) +else: + sys.path.insert(0, _REPO_ROOT) from diff_diff import TwoWayFixedEffects, HAS_RUST_BACKEND -from benchmarks.python.utils import Timer +from benchmarks.python.utils import Timer, collect_provenance def parse_args(): @@ -45,6 +52,10 @@ def parse_args(): "--backend", default="auto", choices=["auto", "python", "rust"], help="Backend to use: auto (default), python (pure Python), rust (Rust backend)" ) + parser.add_argument( + "--warmup", action="store_true", + help="Run one untimed fit before the timed fit (JIT/cache warm-up)" + ) return parser.parse_args() @@ -66,6 +77,12 @@ def main(): # Run benchmark using TwoWayFixedEffects (within-transformation approach) print("Running TWFE estimation...") + if args.warmup: + print("Warm-up fit (untimed)...") + TwoWayFixedEffects(robust=True).fit( + data, outcome="outcome", treatment="treated", time="post", unit="unit" + ) + twfe = TwoWayFixedEffects(robust=True) # auto-clusters at unit level with Timer() as timer: @@ -111,7 +128,10 @@ def main(): "n_units": len(data["unit"].unique()), "n_periods": len(data["post"].unique()), "n_obs": len(data), + "warmup": args.warmup, }, + # Wheel/backend provenance (refresh runs hard-fail on mismatch) + "provenance": collect_provenance(), } # Write output diff --git a/benchmarks/python/utils.py b/benchmarks/python/utils.py index 3ccca803..6d484423 100644 --- a/benchmarks/python/utils.py +++ b/benchmarks/python/utils.py @@ -3,6 +3,8 @@ """ import json +import os +import sys import time from dataclasses import dataclass, field from pathlib import Path @@ -12,6 +14,57 @@ import pandas as pd +# DIFF_DIFF_* knobs whose VALUES are safe to persist in committed benchmark +# artifacts. Any other DIFF_DIFF_* variable is recorded as present-but-redacted +# so a sensitive value can never leak into results JSON. +_BENIGN_PROVENANCE_ENV_KEYS = { + "DIFF_DIFF_BENCH_USE_INSTALLED", + "DIFF_DIFF_BACKEND", + "DIFF_DIFF_DEMEAN_CHUNK_COLS", + "DIFF_DIFF_SOLVE_OLS_FASTPATH", +} + + +def collect_provenance() -> Dict[str, Any]: + """ + Collect provenance metadata for the diff_diff package actually imported. + + Used by the isolated-venv benchmark refresh to hard-fail runs that + silently resolve diff_diff from the dev tree instead of the pinned + wheel. All imports are lazy and defensive so this works on older + released versions too. + """ + import diff_diff + + try: + from diff_diff import HAS_RUST_BACKEND + + has_rust = bool(HAS_RUST_BACKEND) + except ImportError: + has_rust = None + try: + from diff_diff._backend import rust_backend_info + + rbi = rust_backend_info() + except Exception: + rbi = None + return { + "diff_diff_version": getattr(diff_diff, "__version__", None), + "diff_diff_path": getattr(diff_diff, "__file__", None), + "has_rust_backend": has_rust, + "rust_backend_info": rbi, + "python_version": sys.version.split()[0], + "python_executable": sys.executable, + "numpy_version": np.__version__, + "pandas_version": pd.__version__, + "diff_diff_env": { + k: (v if k in _BENIGN_PROVENANCE_ENV_KEYS else "") + for k, v in os.environ.items() + if k.startswith("DIFF_DIFF_") + }, + } + + def compute_timing_stats( timings: List[float], exclude_first: bool = True, diff --git a/benchmarks/refresh_2026_07/gen_benchmark_tables.py b/benchmarks/refresh_2026_07/gen_benchmark_tables.py new file mode 100644 index 00000000..303171ca --- /dev/null +++ b/benchmarks/refresh_2026_07/gen_benchmark_tables.py @@ -0,0 +1,469 @@ +#!/usr/bin/env python3 +""" +Regenerate the numeric tables in docs/benchmarks.rst from refresh results. + +Modeled on benchmarks/speed_review/gen_findings_tables.py: every regenerable +numeric surface on the published page lives between marker comments + + .. refresh-table-start: + .. refresh-table-end: + +and this script rewrites ONLY the content between markers, from the +committed benchmarks/refresh_2026_07/results/refresh_results.json. Prose +stays hand-written. Regions with no data in the results JSON are left +untouched (so a partial run never wipes published tables). + +Usage: + python gen_benchmark_tables.py # rewrite docs/benchmarks.rst + python gen_benchmark_tables.py --check # print diff, write nothing + python gen_benchmark_tables.py --results path.json +""" + +import argparse +import difflib +import json +import re +import sys +from pathlib import Path +from typing import Any, Dict, List, Optional + +REFRESH_DIR = Path(__file__).parent +REPO_ROOT = REFRESH_DIR.parent.parent +DEFAULT_RESULTS = REFRESH_DIR / "results" / "refresh_results.json" +DOCS_RST = REPO_ROOT / "docs" / "benchmarks.rst" + +MARKER_START = ".. refresh-table-start: {key}" +MARKER_END = ".. refresh-table-end: {key}" + +# Hard-gate classification is shared with run_refresh.py via refresh_common +# (single source of truth - the two consumers can never drift). +from refresh_common import hard_flags as _hard_flags # noqa: E402 + + +def assert_uniform_environment(payload: Dict[str, Any]) -> None: + """ + Every rendered cell must come from the same environment/protocol as the + payload's run metadata - a partial --only rerun after an R package + upgrade (or protocol change) can never silently mix with older cells. + """ + expected = payload.get("env_fingerprint") + stale = sorted( + key + for key, cell in payload.get("cells", {}).items() + if cell.get("env_fingerprint") != expected + ) + if expected is None or stale: + raise SystemExit( + "REFUSING generation - cells captured under a different " + f"environment/protocol than the current run metadata: " + f"{stale or 'payload has no env_fingerprint'}. Re-run those " + "cells (or the full refresh) so all rendered numbers share one " + "environment." + ) + + +def assert_no_hard_flags(payload: Dict[str, Any]) -> None: + """A payload with correctness-gated cells must never render as tables.""" + flagged = {} + for key, cell in payload.get("cells", {}).items(): + hard = _hard_flags(cell.get("flags", [])) + if hard: + flagged[key] = hard + if flagged: + detail = "; ".join(f"{k}: {v}" for k, v in sorted(flagged.items())) + raise SystemExit( + "REFUSING to generate tables - hard-gated cells present " + f"(resolve or re-run them first): {detail}" + ) + + +# Ordered estimator specs: (cell key prefix, display name, R package label) +ESTIMATORS = [ + ("basic", "BasicDiD", "R fixest"), + ("twfe", "TWFE (absorbed FE)", "R fixest"), + ("multiperiod", "MultiPeriodDiD", "R fixest"), + ("callaway", "CallawaySantAnna", "R did"), + ("synthdid", "SyntheticDiD", "R synthdid"), +] + +PERF_SCALES = { + "basic": ["small", "1k", "5k", "10k", "20k"], + "twfe": ["small", "1k", "5k", "10k", "20k"], + "callaway": ["small", "1k", "5k", "10k", "20k"], + "synthdid": ["small", "1k", "5k"], +} + + +def fmt_time(t: float) -> str: + if t >= 100: + return f"{t:.0f}" + if t >= 10: + return f"{t:.1f}" + if t >= 1: + return f"{t:.2f}" + if t >= 0.01: + return f"{t:.3f}" + return f"{t:.4f}" + + +def fmt_ratio(r: float) -> str: + if r >= 10: + return f"**{r:.0f}x**" + if r >= 2: + return f"**{r:.1f}x**" + if r >= 1: + return f"{r:.1f}x" + return f"{r:.1f}x" + + +class _MissingData(Exception): + """Raised when a cell/arm/field needed for a region is absent.""" + + +def _arm_median(cell: Dict[str, Any], arm: str) -> float: + entry = cell.get(arm) + if not entry: + raise _MissingData(arm) + return float(entry["timing"]["stats"]["median"]) + + +def _arm_val(cell: Dict[str, Any], arm: str, *keys: str) -> float: + entry = cell.get(arm) + if not entry: + raise _MissingData(arm) + res = entry["result"] + for k in keys: + if k in res and res[k] is not None: + return float(res[k]) + raise _MissingData(f"{arm}:{keys}") + + +def _list_table(widths: List[int], header: List[str], rows: List[List[str]]) -> List[str]: + lines = [ + ".. list-table::", + " :header-rows: 1", + f" :widths: {' '.join(str(w) for w in widths)}", + "", + ] + for row in [header] + rows: + lines.append(f" * - {row[0]}") + for col in row[1:]: + lines.append(f" - {col}") + return lines + + +def _cell_max_arm_diffs(cell: Dict[str, Any]) -> tuple: + """Max ATT abs diff / SE rel diff across BOTH rendered Python arms vs R + (the rust-only comparison would understate a larger-but-passing pure + difference). Falls back to the comparison record if arm data is absent. + """ + try: + r_att = _arm_val(cell, "r", "overall_att", "att") + r_se = _arm_val(cell, "r", "overall_se", "se") + atts = [] + ses = [] + for arm in ("python_pure", "python_rust"): + atts.append(abs(_arm_val(cell, arm, "overall_att", "att") - r_att)) + if r_se: + ses.append(abs(_arm_val(cell, arm, "overall_se", "se") - r_se) / abs(r_se)) + return max(atts), max(ses) if ses else 0.0 + except _MissingData: + return cell["comparison"]["att_diff"], cell["comparison"]["se_rel_diff"] + + +def gen_summary(cells: Dict[str, Any]) -> Optional[List[str]]: + rows = [] + for key, display, _ in ESTIMATORS: + scales = [c for k, c in cells.items() if k.startswith(f"{key}/")] + if not scales: + return None + diffs = [_cell_max_arm_diffs(c) for c in scales] + att = max(d[0] for d in diffs) + se = max(d[1] for d in diffs) + overlap = all(c["comparison"]["ci_overlap"] for c in scales) + passed = all(c["comparison"]["passed"] for c in scales) + rows.append( + [ + display, + f"< {att:.0e}" if att > 0 else "0", + f"{se:.1%}", + "Yes" if overlap else "No", + "**PASS**" if passed else "**FAIL**", + ] + ) + table = _list_table( + [25, 20, 20, 15, 20], ["Estimator", "ATT Diff", "SE Rel Diff", "CI Overlap", "Status"], rows + ) + return table + [ + "", + "SyntheticDiD SE differences reflect Monte Carlo dispersion of the", + "placebo variance (R's placebo permutation is unseeded, so the two", + "implementations agree in distribution, not draw-by-draw): its SE is", + "gated at a 35% relative bound with R's rep-to-rep SE values recorded", + "in the committed results artifact, while the deterministic", + "Frank-Wolfe ATT is gated at 1e-8. All other estimators use analytical", + "SEs gated at the tolerances above. See the SyntheticDiD methodology", + "registry note (benchmark SE gate is Monte Carlo-bounded).", + ] + + +def gen_accuracy(cells: Dict[str, Any], key: str, r_label: str) -> Optional[List[str]]: + cell = cells.get(f"{key}/small") + if not cell: + return None + try: + return _gen_accuracy_inner(cell, r_label) + except _MissingData: + return None + + +def _gen_accuracy_inner(cell: Dict[str, Any], r_label: str) -> List[str]: + pure_att = _arm_val(cell, "python_pure", "overall_att", "att") + rust_att = _arm_val(cell, "python_rust", "overall_att", "att") + r_att = _arm_val(cell, "r", "overall_att", "att") + pure_se = _arm_val(cell, "python_pure", "overall_se", "se") + rust_se = _arm_val(cell, "python_rust", "overall_se", "se") + r_se = _arm_val(cell, "r", "overall_se", "se") + t_pure = _arm_median(cell, "python_pure") + t_rust = _arm_median(cell, "python_rust") + t_r = _arm_median(cell, "r") + + att_diff = max(abs(pure_att - r_att), abs(rust_att - r_att)) + se_rel = max( + abs(pure_se - r_se) / r_se if r_se else 0.0, + abs(rust_se - r_se) / r_se if r_se else 0.0, + ) + t_best = min(t_pure, t_rust) + best_label = "pure" if t_pure <= t_rust else "rust" + if t_r >= t_best: + time_diff = f"**{t_r / t_best:.1f}x faster** ({best_label})" + else: + time_diff = f"{t_best / t_r:.1f}x slower ({best_label})" + + rows = [ + [ + "ATT", + f"{pure_att:.3f}", + f"{rust_att:.3f}", + f"{r_att:.3f}", + f"< {att_diff:.0e}" if att_diff > 0 else "0", + ], + ["SE", f"{pure_se:.3f}", f"{rust_se:.3f}", f"{r_se:.3f}", f"{se_rel:.1%}"], + ["Time (s)", fmt_time(t_pure), fmt_time(t_rust), fmt_time(t_r), time_diff], + ] + return _list_table( + [16, 21, 21, 21, 21], + ["Metric", "diff-diff (Pure)", "diff-diff (Rust)", r_label, "Difference"], + rows, + ) + + +def gen_perf(cells: Dict[str, Any], key: str) -> Optional[List[str]]: + rows = [] + for scale in PERF_SCALES[key]: + cell = cells.get(f"{key}/{scale}") + if not cell: + return None + try: + t_r = _arm_median(cell, "r") + t_pure = _arm_median(cell, "python_pure") + t_rust = _arm_median(cell, "python_rust") + except _MissingData: + return None + rows.append( + [ + scale, + fmt_time(t_r), + fmt_time(t_pure), + fmt_time(t_rust), + fmt_ratio(t_r / t_pure), + fmt_ratio(t_r / t_rust), + ] + ) + return _list_table( + [12, 15, 18, 18, 12, 12], + ["Scale", "R (s)", "Python Pure (s)", "Python Rust (s)", "Pure/R", "Rust/R"], + rows, + ) + + +def gen_mpdta(cells: Dict[str, Any]) -> Optional[List[str]]: + cell = cells.get("mpdta/real") + if not cell: + return None + try: + py_att = _arm_val(cell, "python_rust", "overall_att", "att") + r_att = _arm_val(cell, "r", "overall_att", "att") + py_se = _arm_val(cell, "python_rust", "overall_se", "se") + r_se = _arm_val(cell, "r", "overall_se", "se") + t_py = _arm_median(cell, "python_rust") + t_r = _arm_median(cell, "r") + except _MissingData: + return None + att_diff = abs(py_att - r_att) + n_reps = cell["python_rust"]["timing"]["n_reps"] + rows = [ + [ + "ATT", + f"{py_att:.6f}", + f"{r_att:.6f}", + "**0** (exact match)" if att_diff < 5e-7 else f"{att_diff:.1e}", + ], + [ + "SE (analytical)", + f"{py_se:.4f}", + f"{r_se:.4f}", + f"**< {max(abs(py_se - r_se) / r_se, 0.001):.1%}**", + ], + [ + f"Time (median of {n_reps})", + f"{fmt_time(t_py)}s", + f"{fmt_time(t_r)}s", + f"**{t_r / t_py:.1f}x faster**", + ], + ] + return _list_table([25, 25, 25, 25], ["Metric", "diff-diff", "R did", "Difference"], rows) + + +def gen_environment(payload: Dict[str, Any]) -> Optional[List[str]]: + meta = payload.get("run_metadata") + cells = payload.get("cells", {}) + if not meta or not cells: + return None + # Pull wheel provenance from any python arm. + prov = {} + r_meta = {} + for cell in cells.values(): + entry = cell.get("python_rust") + if entry and entry["result"].get("provenance"): + prov = entry["result"]["provenance"] + r_entry = cell.get("r") + if r_entry and r_entry["result"].get("metadata"): + r_meta = r_entry["result"]["metadata"] + if prov and r_meta: + break + pkgs = meta.get("r_packages", {}) + pkg_str = ", ".join(f"{k} {v}" for k, v in pkgs.items() if k in ("fixest", "did", "synthdid")) + hw = meta["hardware"] + lines = [ + ".. rubric:: Benchmark environment (2026-07 refresh)", + "", + f"- **Captured**: {meta['date_utc']}", + f"- **Hardware**: {hw['cpu']}, {hw['memory_gb']} GB RAM, " f"{hw['os']} ({hw['arch']})", + f"- **diff-diff**: {payload.get('headline_pin')} released wheel " + f"from PyPI (Rust backend + Apple Accelerate), Python " + f"{prov.get('python_version', '?')}, NumPy " + f"{prov.get('numpy_version', '?')}, pandas " + f"{prov.get('pandas_version', '?')}", + f"- **R**: {meta.get('r_version', '?')}; {pkg_str} " f"(installed at capture)", + f"- **Threads**: {meta['thread_policy']}", + f"- **Protocol**: {meta['protocol']}", + ] + return lines + + +def gen_perf_basic_twfe(cells: Dict[str, Any]) -> Optional[List[str]]: + """ + The perf_basic docs region renders BOTH the simple interaction OLS + (BasicDiD) and the genuine absorbed-FE TWFE pair - the legacy page + labeled a single interaction-OLS table "BasicDiD/TWFE". + """ + basic = gen_perf(cells, "basic") + twfe = gen_perf(cells, "twfe") + if basic is None or twfe is None: + return None + return ( + ["**BasicDiD (interaction OLS, clustered):**", ""] + + basic + + ["", "**TWFE (absorbed unit + post fixed effects, clustered):**", ""] + + twfe + ) + + +def build_regions(payload: Dict[str, Any]) -> Dict[str, Optional[List[str]]]: + cells = payload.get("cells", {}) + return { + "summary": gen_summary(cells), + "accuracy_basic": gen_accuracy(cells, "basic", "R fixest"), + "accuracy_multiperiod": gen_accuracy(cells, "multiperiod", "R fixest"), + "accuracy_synthdid": gen_accuracy(cells, "synthdid", "R synthdid"), + "accuracy_callaway": gen_accuracy(cells, "callaway", "R did"), + "environment": gen_environment(payload), + "perf_basic": gen_perf_basic_twfe(cells), + "perf_callaway": gen_perf(cells, "callaway"), + "perf_synthdid": gen_perf(cells, "synthdid"), + "mpdta": gen_mpdta(cells), + } + + +def splice(src: str, key: str, lines: List[str]) -> str: + start = MARKER_START.format(key=key) + end = MARKER_END.format(key=key) + pattern = re.compile( + re.escape(start) + r"\n.*?" + re.escape(end), + flags=re.DOTALL, + ) + if not pattern.search(src): + raise SystemExit(f"marker pair not found in {DOCS_RST}: {key}") + replacement = start + "\n\n" + "\n".join(lines) + "\n\n" + end + return pattern.sub(lambda _: replacement, src, count=1) + + +def main() -> int: + ap = argparse.ArgumentParser(description=__doc__) + ap.add_argument("--results", type=Path, default=DEFAULT_RESULTS) + ap.add_argument( + "--check", action="store_true", help="print diff and exit 1 if the page would change" + ) + ap.add_argument( + "--allow-partial", + action="store_true", + help="tolerate regions with no data (local partial runs); default " + "publication mode fails so fresh and stale tables can never mix", + ) + args = ap.parse_args() + + with open(args.results) as f: + payload = json.load(f) + + assert_no_hard_flags(payload) + if not args.allow_partial: + assert_uniform_environment(payload) + + original = DOCS_RST.read_text() + updated = original + skipped = [] + for key, lines in build_regions(payload).items(): + if lines is None: + skipped.append(key) + continue + updated = splice(updated, key, lines) + if skipped: + if not args.allow_partial: + raise SystemExit( + "REFUSING partial generation - these regions have no data in " + f"{args.results.name}: {', '.join(skipped)}. Re-run the " + "missing cells or pass --allow-partial for local iteration." + ) + print(f"regions left untouched (no data): {', '.join(skipped)}", file=sys.stderr) + + if args.check: + if updated == original: + print("docs/benchmarks.rst is up to date") + return 0 + diff = difflib.unified_diff( + original.splitlines(keepends=True), + updated.splitlines(keepends=True), + fromfile="docs/benchmarks.rst (current)", + tofile="docs/benchmarks.rst (regenerated)", + ) + sys.stdout.writelines(diff) + return 1 + + DOCS_RST.write_text(updated) + print(f"updated {DOCS_RST}") + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/benchmarks/refresh_2026_07/refresh_common.py b/benchmarks/refresh_2026_07/refresh_common.py new file mode 100644 index 00000000..f74d9e03 --- /dev/null +++ b/benchmarks/refresh_2026_07/refresh_common.py @@ -0,0 +1,876 @@ +#!/usr/bin/env python3 +""" +Shared machinery for the 2026-07 public benchmark refresh. + +Used by run_refresh.py (headline docs/benchmarks.rst refresh: released +diff-diff wheel vs R) and run_version_story.py (internal 3.5.3-vs-3.7.0 +optimization story). Both runners execute every replication in a fresh +subprocess, strictly sequentially, against pinned diff-diff wheels installed +in isolated uv venvs - never against the dev tree. + +Key invariants enforced here: +- Provenance hard-fail: a Python arm whose imported diff_diff does not come + from the expected venv/version, or whose backend does not match the arm, + aborts the run (guards against dev-tree sys.path shadowing). +- Sequential execution: one benchmark subprocess at a time, cool-downs + between arms, load-average guard before timed runs. +- Warm-up inside every subprocess (--warmup) so R byte-compiler JIT and + first-call setup stay out of the timing window on BOTH sides. +""" + +import hashlib +import json +import math +import os +import platform +import subprocess +import sys +import time +from datetime import datetime, timezone +from pathlib import Path +from typing import Any, Dict, List, Optional + +REFRESH_DIR = Path(__file__).parent +BENCHMARK_DIR = REFRESH_DIR.parent +REPO_ROOT = BENCHMARK_DIR.parent +VENVS_DIR = REFRESH_DIR / "venvs" +RESULTS_DIR = REFRESH_DIR / "results" +RAW_DIR = RESULTS_DIR / "raw" +DATA_DIR = BENCHMARK_DIR / "data" / "synthetic" # gitignored (*.csv) +MPDTA_SOURCE = BENCHMARK_DIR / "data" / "real" / "mpdta.csv" # committed + +sys.path.insert(0, str(REPO_ROOT)) + +from benchmarks.python.utils import ( # noqa: E402 + compute_timing_stats, + generate_basic_did_data, + generate_multiperiod_data, + generate_sdid_data, + generate_staggered_data, + save_benchmark_data, +) + +UV_BIN = os.environ.get("UV_BIN", os.path.expanduser("~/.local/bin/uv")) +PYTHON_PIN = "3.14" + +# Coefficient-of-variation gate (std/mean of counted reps), matching the +# CV_FLAG convention in benchmarks/speed_review/bench_fe_absorption.py. +CV_FLAG = 0.10 +COOLDOWN_SECONDS = 5 +LOAD_GUARD_1MIN = 2.0 + +# Published MPDTA overall ATT (docs/benchmarks.rst) - known-answer gate. +MPDTA_KNOWN_ATT = -0.039951 + +# Flag tokens that mark a cell as CORRECTNESS-failed (R-parity, known +# answers, detailed effect surfaces) - as opposed to timing-noise flags. +# Single source of truth: run_refresh.py fails the run on these and +# gen_benchmark_tables.py refuses to render payloads containing them. +# Thread-count knobs that would silently change per-arm performance while +# the artifacts claim package defaults - stripped from every benchmark +# subprocess (Python AND R) so "defaults" is enforced, not assumed. +THREAD_ENV_VARS = ( + "RAYON_NUM_THREADS", + "OMP_NUM_THREADS", + "OPENBLAS_NUM_THREADS", + "VECLIB_MAXIMUM_THREADS", + "MKL_NUM_THREADS", + "R_DATATABLE_NUM_THREADS", +) + +HARD_FLAG_TOKENS = ( + "parity_fail", + "known_answer_fail", + "att_gate_fail", + "se_gate_fail", + "ci_gate_fail", + "weights_gate_fail", + "keys_mismatch", +) + + +def hard_flags(flags: List[str]) -> List[str]: + # Subset of flags that are hard correctness gates. + return [fl for fl in flags if any(tok in fl for tok in HARD_FLAG_TOKENS)] + + +def utc_now() -> str: + return datetime.now(timezone.utc).isoformat(timespec="seconds") + + +def env_fingerprint( + meta: Dict[str, Any], pin: str, python_env: Optional[Dict[str, Any]] = None +) -> str: + """ + Short hash of everything that must be CONSTANT across all cells rendered + onto the docs page together: package versions, hardware, pin, protocol. + Excludes timestamps/load. Publication refuses to mix fingerprints, so a + partial --only rerun under a changed environment can never silently blend + with older cells. + """ + basis = { + "pin": pin, + "r_version": meta.get("r_version"), + "r_packages": meta.get("r_packages"), + "cpu": (meta.get("hardware") or {}).get("cpu"), + "os": (meta.get("hardware") or {}).get("os"), + "orchestrator_python": meta.get("orchestrator_python"), + "protocol": meta.get("protocol"), + "thread_policy": meta.get("thread_policy"), + # Python-arm provenance (venv python/numpy/pandas + BLAS linkage): + # a dependency upgrade between partial reruns must change the + # fingerprint so mixed-environment cells can never blend. + "python_env": python_env, + } + return hashlib.sha256(json.dumps(basis, sort_keys=True).encode()).hexdigest()[:16] + + +def compare_weight_vectors( + py_weights: Any, + r_weights: Any, + atol: float = 1e-8, + py_ids: Any = None, + r_ids: Any = None, + require_ids: bool = False, +) -> Dict[str, Any]: + """ + Compare SDID unit/time weight vectors, aligning by unit/period id when + both sides provide ids (ordering-robust: Python's + get_unit_weights_df() sorts by descending weight, R emits panel order). + With require_ids=True (the SDID publication gate), missing ids on + either side FAILS the comparison - the documented id-alignment contract + can never silently degrade to positional order if a benchmark script + regresses. Without it, positional comparison is a permitted fallback + for legacy artifacts. Fail-closed on length/key mismatch, duplicates, + or non-finite entries; metrics are committed for audit. + """ + py = list(py_weights or []) + r = list(r_weights or []) + metrics: Dict[str, Any] = { + "n_python": len(py), + "n_r": len(r), + "max_abs_diff": None, + "aligned_by_ids": False, + "ok": False, + } + if require_ids and (py_ids is None or r_ids is None): + return metrics + if not py or len(py) != len(r): + return metrics + + def _key(x): + return round(float(x), 9) + + if py_ids is not None and r_ids is not None: + py_ids = list(py_ids) + r_ids = list(r_ids) + if len(py_ids) != len(py) or len(r_ids) != len(r): + return metrics + try: + py_map = {_key(k): w for k, w in zip(py_ids, py)} + r_map = {_key(k): w for k, w in zip(r_ids, r)} + except (TypeError, ValueError): + return metrics + if len(py_map) != len(py) or len(r_map) != len(r): + return metrics # duplicate ids + if set(py_map) != set(r_map): + return metrics # different unit/period sets + pairs = [(py_map[k], r_map[k]) for k in sorted(py_map)] + metrics["aligned_by_ids"] = True + else: + pairs = list(zip(py, r)) + + diffs = [] + for a, b in pairs: + try: + a_f, b_f = float(a), float(b) + except (TypeError, ValueError): + return metrics + if not (math.isfinite(a_f) and math.isfinite(b_f)): + return metrics + diffs.append(abs(a_f - b_f)) + metrics["max_abs_diff"] = max(diffs) + metrics["ok"] = metrics["max_abs_diff"] < atol + return metrics + + +def validate_estimator_field(result: Dict[str, Any], expected: str) -> None: + """ + Hard-fail when a benchmark subprocess ran a different estimator than the + cell spec claims (e.g. a silently ignored --type flag). The emitted + "estimator" field is distinct per script/model, so miswiring can never + publish numbers under the wrong label. + """ + actual = result.get("estimator") + if actual != expected: + raise RuntimeError( + f"estimator mismatch: spec expects {expected!r}, " + f"benchmark subprocess reported {actual!r}" + ) + + +def headline_gate_flags( + arm_name: str, + py_att: Any, + py_se: Any, + r_att: Any, + r_se: Any, + att_atol: float, + se_rtol: float, +) -> List[str]: + """ + Strict headline gates for one Python arm vs R: ATT and SE against the + documented tolerances, plus the published CI-overlap criterion - each + enforced independently (overlapping CIs never waive an out-of-tolerance + SE, and every RENDERED arm is gated, not just the rust-vs-R comparison). + Non-finite values always gate. Returned flags carry HARD_FLAG_TOKENS + substrings. + """ + + def _finite(x) -> Optional[float]: + try: + v = float(x) + except (TypeError, ValueError): + return None + return v if math.isfinite(v) else None + + flags: List[str] = [] + py_att_f, r_att_f = _finite(py_att), _finite(r_att) + py_se_f, r_se_f = _finite(py_se), _finite(r_se) + if py_att_f is None or r_att_f is None: + flags.append(f"headline_att_gate_fail:{arm_name}:nonfinite") + else: + att_diff = abs(py_att_f - r_att_f) + att_rel = att_diff / (abs(r_att_f) + 1e-10) + if not (att_diff < att_atol or att_rel < 0.01): + flags.append(f"headline_att_gate_fail:{arm_name}:{att_diff:.2e}") + if py_se_f is None or r_se_f is None: + flags.append(f"headline_se_gate_fail:{arm_name}:nonfinite") + else: + se_rel = abs(py_se_f - r_se_f) / (abs(r_se_f) + 1e-10) + if se_rel >= se_rtol: + flags.append(f"headline_se_gate_fail:{arm_name}:{se_rel:.3f}") + if None not in (py_att_f, py_se_f, r_att_f, r_se_f): + # Published criterion: 95% CIs must overlap (same formula as + # compare_estimates); non-finite cases are already gated above. + py_lo, py_hi = py_att_f - 1.96 * py_se_f, py_att_f + 1.96 * py_se_f + r_lo, r_hi = r_att_f - 1.96 * r_se_f, r_att_f + 1.96 * r_se_f + overlap = (py_lo <= r_att_f <= py_hi) or (r_lo <= py_att_f <= r_hi) + if not overlap: + flags.append(f"ci_gate_fail:{arm_name}") + return flags + + +# --------------------------------------------------------------------------- +# Venv management +# --------------------------------------------------------------------------- + + +def venv_python(venv_name: str) -> Path: + return VENVS_DIR / venv_name / "bin" / "python" + + +def setup_venv(venv_name: str, pin: str) -> None: + """Create an isolated venv with a pinned diff-diff release wheel.""" + vdir = VENVS_DIR / venv_name + py = venv_python(venv_name) + if not py.exists(): + print(f"[setup] creating venv {vdir} (python {PYTHON_PIN})") + subprocess.run( + [UV_BIN, "venv", str(vdir), "--python", PYTHON_PIN], + check=True, + ) + print(f"[setup] installing diff-diff=={pin} into {venv_name}") + subprocess.run( + [UV_BIN, "pip", "install", "--python", str(py), f"diff-diff=={pin}"], + check=True, + ) + + +def preflight_venv(venv_name: str, pin: str) -> Dict[str, Any]: + """Import diff_diff inside the venv and assert wheel provenance.""" + code = ( + "import json, diff_diff\n" + "from diff_diff import HAS_RUST_BACKEND\n" + "try:\n" + " from diff_diff._backend import rust_backend_info\n" + " rbi = rust_backend_info()\n" + "except Exception:\n" + " rbi = None\n" + "import numpy, pandas, sys\n" + "print(json.dumps({'version': diff_diff.__version__," + " 'path': diff_diff.__file__, 'has_rust': bool(HAS_RUST_BACKEND)," + " 'rust_backend_info': rbi, 'python': sys.version.split()[0]," + " 'numpy': numpy.__version__, 'pandas': pandas.__version__}))\n" + ) + env = _child_env() + proc = subprocess.run( + [str(venv_python(venv_name)), "-c", code], + capture_output=True, + text=True, + env=env, + cwd=str(RAW_DIR), + ) + if proc.returncode != 0: + raise RuntimeError(f"preflight failed for {venv_name}: {proc.stderr.strip()}") + info = json.loads(proc.stdout.strip().splitlines()[-1]) + if info["version"] != pin: + raise RuntimeError(f"{venv_name}: expected diff-diff=={pin}, got {info['version']}") + if str(VENVS_DIR / venv_name) not in info["path"]: + raise RuntimeError(f"{venv_name}: diff_diff resolved OUTSIDE the venv: {info['path']}") + if not info["has_rust"]: + raise RuntimeError(f"{venv_name}: wheel has no Rust backend") + rbi = info.get("rust_backend_info") + if rbi is not None and not rbi.get("accelerate"): + raise RuntimeError(f"{venv_name}: wheel is not linked against Accelerate: {rbi}") + print( + f"[preflight] {venv_name}: diff-diff {info['version']} " + f"(python {info['python']}, numpy {info['numpy']}, " + f"pandas {info['pandas']}, rust+accelerate OK)" + ) + return info + + +# --------------------------------------------------------------------------- +# Subprocess execution +# --------------------------------------------------------------------------- + + +def _child_env() -> Dict[str, str]: + """ + Child env: no PYTHONPATH (installed wheel wins over the dev tree) and no + inherited DIFF_DIFF_* knobs - a stray DIFF_DIFF_DEMEAN_CHUNK_COLS or + similar in the parent shell would silently change benchmark code paths + while the artifacts claim wheel defaults. The benchmark script itself + sets DIFF_DIFF_BACKEND inside the subprocess from --backend. + """ + env = { + k: v + for k, v in os.environ.items() + if not k.startswith("DIFF_DIFF_") and k not in THREAD_ENV_VARS + } + env.pop("PYTHONPATH", None) + env["DIFF_DIFF_BENCH_USE_INSTALLED"] = "1" + return env + + +def run_python_rep( + script_name: str, + data_path: Path, + out_path: Path, + backend: str, + venv: str, + extra_args: Optional[List[str]] = None, + timeout: Optional[int] = None, + warmup: bool = True, +) -> Dict[str, Any]: + """One Python replication in a fresh subprocess using a pinned venv.""" + cmd = [ + str(venv_python(venv)), + str(BENCHMARK_DIR / "python" / script_name), + "--data", + str(data_path), + "--output", + str(out_path), + "--backend", + backend, + ] + if warmup: + cmd.append("--warmup") + if extra_args: + cmd.extend(extra_args) + proc = subprocess.run( + cmd, + capture_output=True, + text=True, + timeout=timeout, + env=_child_env(), + cwd=str(RAW_DIR), + ) + if proc.returncode != 0: + raise RuntimeError( + f"python rep failed ({script_name}, backend={backend}, venv={venv}):\n" + f"stdout: {proc.stdout[-2000:]}\nstderr: {proc.stderr[-2000:]}" + ) + with open(out_path) as f: + return json.load(f) + + +def run_r_rep( + script_name: str, + data_path: Path, + out_path: Path, + extra_args: Optional[List[str]] = None, + timeout: Optional[int] = None, + warmup: bool = True, +) -> Dict[str, Any]: + """One R replication in a fresh Rscript subprocess.""" + cmd = [ + "Rscript", + # --vanilla: no user/site .Rprofile or .Renviron - those could set + # fixest/data.table thread counts behind the enforced-defaults + # claim. Package discovery is unaffected (built-in R_LIBS_USER + # default applies without .Renviron). + "--vanilla", + str(BENCHMARK_DIR / "R" / script_name), + "--data", + str(data_path), + "--output", + str(out_path), + "--warmup", + "true" if warmup else "false", + ] + if extra_args: + cmd.extend(extra_args) + proc = subprocess.run( + cmd, + capture_output=True, + text=True, + timeout=timeout, + env=_child_env(), + cwd=str(RAW_DIR), + ) + if proc.returncode != 0: + raise RuntimeError( + f"R rep failed ({script_name}):\n" + f"stdout: {proc.stdout[-2000:]}\nstderr: {proc.stderr[-2000:]}" + ) + with open(out_path) as f: + return json.load(f) + + +def validate_python_provenance( + result: Dict[str, Any], + expected_version: str, + expected_venv: str, + backend_arm: str, +) -> None: + """Hard-fail if the subprocess imported the wrong diff_diff.""" + prov = result.get("provenance") + if not prov: + raise RuntimeError("no provenance block in Python result JSON") + if prov.get("diff_diff_version") != expected_version: + raise RuntimeError( + f"provenance: version {prov.get('diff_diff_version')} != pin " f"{expected_version}" + ) + path = prov.get("diff_diff_path") or "" + if str(VENVS_DIR / expected_venv) not in path: + raise RuntimeError( + f"provenance: diff_diff imported from OUTSIDE venv " + f"{expected_venv}: {path} (dev-tree shadowing?)" + ) + has_rust = prov.get("has_rust_backend") + if backend_arm == "python" and has_rust: + raise RuntimeError("provenance: pure-python arm ran WITH rust backend") + if backend_arm == "rust" and not has_rust: + raise RuntimeError("provenance: rust arm ran WITHOUT rust backend") + if backend_arm == "rust": + rbi = prov.get("rust_backend_info") + if rbi is not None and not rbi.get("accelerate"): + raise RuntimeError(f"provenance: rust arm not linked against Accelerate: {rbi}") + # No unexpected DIFF_DIFF_* knobs may reach the benchmark subprocess: + # the artifacts claim wheel defaults (plus the explicit backend arm). + env_rec = prov.get("diff_diff_env") or {} + allowed = {"DIFF_DIFF_BENCH_USE_INSTALLED", "DIFF_DIFF_BACKEND"} + unexpected = sorted(set(env_rec) - allowed) + if unexpected: + raise RuntimeError( + f"provenance: unexpected DIFF_DIFF_* knobs in benchmark env: {unexpected}" + ) + if env_rec.get("DIFF_DIFF_BENCH_USE_INSTALLED") != "1": + raise RuntimeError( + "provenance: DIFF_DIFF_BENCH_USE_INSTALLED not active - " + "the dev-tree path guard was off" + ) + backend_env = env_rec.get("DIFF_DIFF_BACKEND") + if backend_arm == "python" and backend_env != "python": + raise RuntimeError( + f"provenance: pure arm expected DIFF_DIFF_BACKEND=python, got {backend_env!r}" + ) + if backend_arm == "rust" and backend_env not in (None, "rust"): + raise RuntimeError( + f"provenance: rust arm expected DIFF_DIFF_BACKEND unset/rust, got {backend_env!r}" + ) + + +# --------------------------------------------------------------------------- +# Arm-level replication loop +# --------------------------------------------------------------------------- + + +def run_arm( + label: str, + rep_fn, + n_reps: int, + allow_cv_rerun: bool = True, +) -> Dict[str, Any]: + """ + Run n_reps sequential subprocess replications of one arm. + + rep_fn() -> parsed result JSON of a single replication. Returns a dict + with the last result, timing stats over all reps (first excluded), the + per-rep ATT determinism check, and CV flag status. + """ + + def _one_pass() -> Dict[str, Any]: + timings: List[float] = [] + atts: List[float] = [] + ses: List[float] = [] + n_nonfinite_att = 0 + n_nonfinite_se = 0 + result = None + for rep in range(n_reps): + result = rep_fn() + t = result["timing"]["total_seconds"] + timings.append(t) + att = result.get("overall_att", result.get("att")) + if att is not None: + att_f = float(att) + if math.isfinite(att_f): + atts.append(att_f) + else: + # Fail-closed: a NaN/Inf replication is a correctness + # failure even if the LAST rep (which feeds the headline + # gates) happens to be finite. + n_nonfinite_att += 1 + se = result.get("overall_se", result.get("se")) + if se is not None: + se_f = float(se) + if math.isfinite(se_f): + ses.append(se_f) + else: + n_nonfinite_se += 1 + print(f" [{label}] rep {rep + 1}/{n_reps}: {t:.3f}s") + stats = compute_timing_stats(timings) + att_spread = (max(atts) - min(atts)) if atts else 0.0 + cv = stats["stats"]["std"] / stats["stats"]["mean"] if stats["stats"].get("mean") else 0.0 + return { + "result": result, + "timing": stats, + "cv": cv, + "att_spread": att_spread, + "n_nonfinite_att": n_nonfinite_att, + "n_nonfinite_se": n_nonfinite_se, + # Rep-to-rep SE values: seeded arms are constant; the unseeded R + # synthdid placebo varies - recorded so the documented SE gate + # tolerance is auditable against actual Monte Carlo dispersion. + "se_values": ses, + } + + def _correctness_flags(o: Dict[str, Any]) -> List[str]: + f: List[str] = [] + if o["att_spread"] > 1e-12: + # Same data + seeded estimators: ATT must be identical across + # reps. Non-deterministic point estimates are a correctness + # failure - the token makes this a HARD publication gate. + f.append(f"rep_att_gate_fail:{o['att_spread']:.2e}") + if o["n_nonfinite_att"]: + f.append(f"rep_att_gate_fail:nonfinite:{o['n_nonfinite_att']}") + if o["n_nonfinite_se"]: + f.append(f"rep_se_gate_fail:nonfinite:{o['n_nonfinite_se']}") + return f + + out = _one_pass() + flags: List[str] = _correctness_flags(out) + if out["cv"] > CV_FLAG and allow_cv_rerun and n_reps > 2: + if flags: + # A hard correctness failure was already observed; a timing + # rerun would DISCARD that evidence if the second pass came + # back clean. Keep the failed pass and let the gates fire. + print( + f" [{label}] CV {out['cv']:.1%} > {CV_FLAG:.0%} but " + f"correctness flags present - skipping rerun: {flags}" + ) + flags.append(f"cv_flag:{out['cv']:.3f}") + else: + print(f" [{label}] CV {out['cv']:.1%} > {CV_FLAG:.0%} - rerunning arm once") + time.sleep(COOLDOWN_SECONDS) + out = _one_pass() + flags.extend(_correctness_flags(out)) + if out["cv"] > CV_FLAG: + flags.append(f"cv_flag:{out['cv']:.3f}") + print(f" [{label}] still flagged (CV {out['cv']:.1%})") + elif out["cv"] > CV_FLAG: + flags.append(f"cv_flag:{out['cv']:.3f}") + out["flags"] = flags + return out + + +def compare_effect_arrays( + py_effects: List[Dict[str, Any]], + r_effects: List[Dict[str, Any]], + join_keys: List[str], + se_rtol: float, + att_atol: float = 1e-6, +) -> Dict[str, Any]: + """ + Align two per-effect arrays (period/group-time/event-study/group) on + join_keys and compare att/se. Returns compact metrics for the committed + results JSON; the caller turns failures into hard gate flags. + + This is what enforces the docs' detailed parity claims (per-period and + per-(g,t)/event-time effects matching R) - the headline ATT/SE comparison + alone cannot certify them. + """ + + def key_of(row: Dict[str, Any]): + return tuple(round(float(row[k]), 9) for k in join_keys) + + def build(rows): + out = {} + dups = 0 + dropped = 0 + for row in rows: + att = row.get("att") + if att is None or not math.isfinite(float(att)): + dropped += 1 # fail-closed: a non-finite effect is a defect + continue + k = key_of(row) + if k in out: + dups += 1 # duplicate join key would silently overwrite + out[k] = row + return out, dups, dropped + + py_map, py_dups, py_dropped = build(py_effects or []) + r_map, r_dups, r_dropped = build(r_effects or []) + common = sorted(set(py_map) & set(r_map)) + keys_clean = ( + len(common) > 0 + and len(py_map) == len(common) + and len(r_map) == len(common) + and py_dups == r_dups == 0 + and py_dropped == r_dropped == 0 + ) + metrics: Dict[str, Any] = { + "n_python": len(py_map), + "n_r": len(r_map), + "n_common": len(common), + "n_only_python": len(py_map) - len(common), + "n_only_r": len(r_map) - len(common), + "n_dup_python": py_dups, + "n_dup_r": r_dups, + "n_dropped_python": py_dropped, + "n_dropped_r": r_dropped, + "n_se_compared": 0, + "max_att_diff": None, + "max_se_rel_diff": None, + "keys_match": keys_clean, + "att_ok": False, + "se_ok": False, + } + if not common: + return metrics + max_att = max(abs(float(py_map[k]["att"]) - float(r_map[k]["att"])) for k in common) + se_rels = [] + for k in common: + py_se = py_map[k].get("se") + r_se = r_map[k].get("se") + if py_se is None or r_se is None: + continue + py_se, r_se = float(py_se), float(r_se) + if math.isfinite(py_se) and math.isfinite(r_se) and r_se != 0.0: + se_rels.append(abs(py_se - r_se) / abs(r_se)) + max_se_rel = max(se_rels) if se_rels else 0.0 + metrics["n_se_compared"] = len(se_rels) + metrics["max_att_diff"] = max_att + metrics["max_se_rel_diff"] = max_se_rel + metrics["att_ok"] = max_att < att_atol + # Fail-closed: every common row must contribute a comparable finite SE. + metrics["se_ok"] = len(se_rels) == len(common) and max_se_rel < se_rtol + return metrics + + +def slim_result(result: Dict[str, Any]) -> Dict[str, Any]: + """Drop bulky per-observation fields before committing to results JSON.""" + slim = dict(result) + prov = slim.get("provenance") + if isinstance(prov, dict): + # Full paths are needed for runtime provenance validation (which + # runs BEFORE slimming); the committed artifact keeps only the + # trailing components so local usernames never land in the repo. + prov = dict(prov) + for key in ("diff_diff_path", "python_executable"): + val = prov.get(key) + if isinstance(val, str) and val: + prov[key] = ".../" + "/".join(Path(val).parts[-4:]) + slim["provenance"] = prov + for key in ( + "unit_weights", + "unit_weight_ids", + "time_weights", + "time_weight_ids", + "group_time_effects", + "event_study", + "group_effects", + "period_effects", + "coefficients", + ): + if key in slim: + val = slim.pop(key) + try: + slim[f"n_{key}"] = len(val) + except TypeError: + pass + return slim + + +# --------------------------------------------------------------------------- +# Machine state / metadata +# --------------------------------------------------------------------------- + + +def check_load(force: bool, enforce: bool = True) -> float: + load1 = os.getloadavg()[0] + if load1 > LOAD_GUARD_1MIN: + msg = f"1-min load average {load1:.2f} > {LOAD_GUARD_1MIN} - " f"the machine is not idle" + if enforce and not force: + raise SystemExit(f"ABORT: {msg}. Re-run with --force to override.") + print(f"WARNING: {msg}") + return load1 + + +def _cmd_out(cmd: List[str]) -> Optional[str]: + try: + return subprocess.run( + cmd, capture_output=True, text=True, check=True, env=_child_env() + ).stdout.strip() + except Exception: + return None + + +def collect_run_metadata() -> Dict[str, Any]: + r_pkgs = _cmd_out( + [ + "Rscript", + "--vanilla", + "-e", + 'for (p in c("did","fixest","synthdid","jsonlite","data.table"))' + ' cat(p, as.character(packageVersion(p)), "\\n")', + ] + ) + packages = {} + if r_pkgs: + for line in r_pkgs.splitlines(): + parts = line.split() + if len(parts) == 2: + packages[parts[0]] = parts[1] + mem_bytes = _cmd_out(["sysctl", "-n", "hw.memsize"]) + return { + "date_utc": datetime.now(timezone.utc).isoformat(timespec="seconds"), + "hardware": { + "cpu": _cmd_out(["sysctl", "-n", "machdep.cpu.brand_string"]), + "memory_gb": round(int(mem_bytes) / 2**30) if mem_bytes else None, + "os": f"macOS {platform.mac_ver()[0]}", + "arch": platform.machine(), + }, + "orchestrator_python": sys.version.split()[0], + "repo_git_sha": _cmd_out(["git", "-C", str(REPO_ROOT), "rev-parse", "HEAD"]), + "r_version": _cmd_out(["Rscript", "--vanilla", "-e", "cat(R.version.string)"]), + "r_packages": packages, + "thread_policy": ( + "No arm is thread-restricted: R runs at fixest/data.table " + "defaults; diff-diff wheels run at Accelerate/rayon defaults. " + "Thread-count env vars (RAYON/OMP/OPENBLAS/VECLIB/MKL/" + "data.table) are stripped from every benchmark subprocess and " + "R runs under --vanilla (no user/site .Rprofile/.Renviron), so " + "package defaults are enforced, not assumed. Per-arm thread " + "counts are recorded in each result's metadata." + ), + "loadavg_at_start": list(os.getloadavg()), + "protocol": ( + "Each replication is a fresh subprocess run strictly " + "sequentially (one benchmark process on the machine at a time) " + "with an untimed in-process warm-up fit before the timed fit. " + "The first replication is additionally excluded from statistics. " + "Published statistic: median of the counted replications. " + f"Arms with CV > {CV_FLAG:.0%} are rerun once and flagged if " + "still noisy." + ), + } + + +# --------------------------------------------------------------------------- +# Data generation +# --------------------------------------------------------------------------- + + +def sha256_file(path: Path) -> str: + h = hashlib.sha256() + with open(path, "rb") as f: + for chunk in iter(lambda: f.read(1 << 20), b""): + h.update(chunk) + return h.hexdigest() + + +_GENERATORS = { + "basic": lambda cfg: generate_basic_did_data(treatment_effect=5.0, seed=42, **cfg), + "staggered": lambda cfg: generate_staggered_data(treatment_effect=2.0, seed=42, **cfg), + "sdid": lambda cfg: generate_sdid_data(treatment_effect=4.0, seed=42, **cfg), + "multiperiod": lambda cfg: generate_multiperiod_data(treatment_effect=3.0, seed=42, **cfg), +} + + +def ensure_dataset(kind: str, scale: str, cfg: Dict[str, int]) -> Path: + """ + Deterministically (re)generate a synthetic dataset CSV shared by all arms. + + The dataset is ALWAYS regenerated in memory and compared against any + existing file; a stale on-disk CSV (older generator, different config) + can never be silently benchmarked under seed-42 claims. Byte-identical + regeneration keeps the file untouched. + """ + path = DATA_DIR / f"{kind}_{scale}.csv" + df = _GENERATORS[kind](cfg) + csv_bytes = df.to_csv(index=False).encode() + digest = hashlib.sha256(csv_bytes).hexdigest() + if path.exists() and sha256_file(path) == digest: + return path + if path.exists(): + print(f"[data] STALE {path.name} - overwriting with deterministic regeneration") + else: + print(f"[data] generating {path.name} ({cfg})") + path.parent.mkdir(parents=True, exist_ok=True) + path.write_bytes(csv_bytes) + return path + + +def ensure_mpdta() -> Path: + """ + Normalize the committed MPDTA CSV (R `did` schema) to the staggered + benchmark schema so the existing callaway/did benchmark pair runs it + unmodified. Asserts the known shape so a wrong/partial source file can + never silently produce a fake benchmark. + """ + import pandas as pd + + path = DATA_DIR / "mpdta_refresh.csv" + df = pd.read_csv(MPDTA_SOURCE) + df = df.rename( + columns={ + "countyreal": "unit", + "year": "time", + "lemp": "outcome", + "first.treat": "first_treat", + } + )[["unit", "time", "outcome", "first_treat"]] + if len(df) != 2500: + raise RuntimeError(f"MPDTA rows: {len(df)} != 2500") + if df["unit"].nunique() != 500: + raise RuntimeError("MPDTA counties != 500") + if set(df["first_treat"].unique()) != {0, 2004, 2006, 2007}: + raise RuntimeError("MPDTA cohorts unexpected") + save_benchmark_data(df, path) + return path + + +def data_determinism_check() -> None: + """Same seed -> byte-identical CSV content (checked at small scale).""" + a = _GENERATORS["basic"]({"n_units": 100, "n_periods": 4}).to_csv(index=False) + b = _GENERATORS["basic"]({"n_units": 100, "n_periods": 4}).to_csv(index=False) + if hashlib.sha256(a.encode()).hexdigest() != hashlib.sha256(b.encode()).hexdigest(): + raise RuntimeError("data generation is NOT deterministic at seed 42") + print("[preflight] data generation determinism OK") + + +def cooldown(seconds: int = COOLDOWN_SECONDS) -> None: + time.sleep(seconds) diff --git a/benchmarks/refresh_2026_07/run_refresh.py b/benchmarks/refresh_2026_07/run_refresh.py new file mode 100644 index 00000000..474bb7b3 --- /dev/null +++ b/benchmarks/refresh_2026_07/run_refresh.py @@ -0,0 +1,502 @@ +#!/usr/bin/env python3 +""" +Headline orchestrator for the 2026-07 public benchmark refresh. + +Re-runs the published docs/benchmarks.rst comparisons - BasicDiD/TWFE vs +fixest, MultiPeriodDiD vs fixest, CallawaySantAnna vs R did, SyntheticDiD vs +R synthdid, and the MPDTA real-data validation - with a fair protocol: + +- diff-diff arm = the RELEASED wheel (pip install diff-diff==) in an + isolated uv venv, never the dev tree (provenance hard-fail). +- R arm at package defaults (fixest/data.table threads untouched), latest + CRAN versions recorded in metadata. +- Untimed in-process warm-up on BOTH sides; fresh subprocess per replication; + strictly sequential; median of counted reps published; CV > 10% flagged. +- SDID bootstrap parity: Python placebo n_bootstrap=200 == R vcov placebo + default 200 (the old orchestrator passed 50 - unfair to R). + +Usage: + python run_refresh.py --setup # create venvs + preflight only + python run_refresh.py --smoke # small-scale pipeline validation + python run_refresh.py # FULL timed run (idle machine!) + python run_refresh.py --only synthdid # one estimator +""" + +import argparse +import json +import os +import sys +import time +from dataclasses import asdict +from pathlib import Path +from typing import Any, Dict, List + +import refresh_common as rc + +sys.path.insert(0, str(rc.REPO_ROOT)) + +from benchmarks.compare_results import compare_estimates # noqa: E402 +from benchmarks.run_benchmarks import ( # noqa: E402 + SCALE_CONFIGS, + TIMEOUT_CONFIGS, +) + +HEADLINE_PIN = "3.7.0" +HEADLINE_VENV = "dd370" + +RESULTS_PATH = rc.RESULTS_DIR / "refresh_results.json" +SMOKE_RESULTS_PATH = rc.RESULTS_DIR / "refresh_results_smoke.json" + +# Published-page scope (user-locked 2026-07-10): existing tables only. +BENCH_SPECS: Dict[str, Dict[str, Any]] = { + # The legacy harness labeled this cell "BasicDiD/TWFE" while both arms + # actually ran the simple interaction OLS (--type twfe is accepted but + # changes neither model). The refresh splits it honestly: "basic" is the + # 2x2 interaction regression, "twfe" is the genuine absorbed-FE pair. + "basic": { + "display": "BasicDiD", + "py_script": "benchmark_basic.py", + "r_script": "benchmark_fixest.R", + "dataset": "basic", + "scales": ["small", "1k", "5k", "10k", "20k"], + "py_extra": ["--type", "basic"], + "r_extra": ["--type", "basic"], + "py_estimator": "diff_diff.DifferenceInDifferences", + "r_estimator": "fixest::feols", + "se_rtol": 0.01, + "se_gate_rtol": 0.01, + }, + "twfe": { + "display": "TWFE (absorbed FE)", + "py_script": "benchmark_twfe.py", + "r_script": "benchmark_twfe.R", + "dataset": "basic", + "scales": ["small", "1k", "5k", "10k", "20k"], + "py_estimator": "diff_diff.TwoWayFixedEffects", + "r_estimator": "fixest::feols (absorbed FE)", + "se_rtol": 0.01, + "se_gate_rtol": 0.01, + }, + "multiperiod": { + "display": "MultiPeriodDiD", + "py_script": "benchmark_multiperiod.py", + "r_script": "benchmark_multiperiod.R", + "dataset": "multiperiod", + # The published page carries MultiPeriodDiD at small scale only. + "scales": ["small"], + "py_estimator": "diff_diff.MultiPeriodDiD", + "r_estimator": "fixest::feols (multiperiod)", + "se_rtol": 0.01, + "se_gate_rtol": 0.01, + }, + "callaway": { + "display": "CallawaySantAnna", + "py_script": "benchmark_callaway.py", + "r_script": "benchmark_did.R", + "dataset": "staggered", + "scales": ["small", "1k", "5k", "10k", "20k"], + "py_estimator": "diff_diff.CallawaySantAnna", + "r_estimator": "did::att_gt", + "se_rtol": 0.10, + "se_gate_rtol": 0.10, + }, + "mpdta": { + "display": "CallawaySantAnna (MPDTA real data)", + "py_script": "benchmark_callaway.py", + "r_script": "benchmark_did.R", + "dataset": "mpdta", + "scales": ["real"], + "py_estimator": "diff_diff.CallawaySantAnna", + "r_estimator": "did::att_gt", + "se_rtol": 0.02, + "se_gate_rtol": 0.02, + }, + # Slowest last: SDID R cells dominate wall-clock. + "synthdid": { + "display": "SyntheticDiD", + "py_script": "benchmark_synthdid.py", + "r_script": "benchmark_synthdid.R", + "dataset": "sdid", + "scales": ["small", "1k", "5k"], + # Placebo-bootstrap parity with R's vcov(method="placebo") default. + "py_extra": ["--n-bootstrap", "200"], + # Jackknife is excluded from total_seconds anyway; skip the wasted + # wall-clock at scale. + "r_extra": ["--skip-jackknife", "true"], + "py_estimator": "diff_diff.SyntheticDiD", + "r_estimator": "synthdid::synthdid_estimate", + "se_rtol": 0.10, + # SDID SE gate is wider than the reporting rtol: both sides estimate + # the placebo variance by Monte Carlo and R's placebo draw is + # UNSEEDED, so its SE varies run to run (rep-to-rep values recorded + # in se_values for audit). 0.35 bounds MC dispersion at the small-N0 + # scale; ATT remains gated at 1e-8 (deterministic Frank-Wolfe). + "se_gate_rtol": 0.35, + "slow_scales": {"1k", "5k"}, + }, +} + +ARMS = ("python_pure", "python_rust", "r") + +# Detailed parity surfaces the docs claim match R beyond the headline ATT/SE +# (period effects, group-time effects, event study, group aggregation). Each +# is aligned on its join keys and hard-gated; compact metrics land in the +# committed results JSON under "detail_comparison". +DETAIL_SURFACES: Dict[str, List] = { + "multiperiod": [("period_effects", ["period"])], + "callaway": [ + ("group_time_effects", ["group", "time"]), + ("event_study", ["event_time"]), + ("group_effects", ["group"]), + ], + "mpdta": [ + ("group_time_effects", ["group", "time"]), + ("event_study", ["event_time"]), + ("group_effects", ["group"]), + ], +} + + +def cell_timeouts(estimator: str, scale: str) -> Dict[str, int]: + key = "small" if scale == "real" else scale + t = dict(TIMEOUT_CONFIGS.get(key, TIMEOUT_CONFIGS["small"])) + if estimator == "synthdid": + # Warm-up doubles per-process work; R 5k budget 2x -> 7200s. + t["r"] = t["r"] * 2 + t["python"] = t["python"] * 2 + return t + + +def multiperiod_extra(scale: str) -> List[str]: + cfg = SCALE_CONFIGS[scale]["multiperiod"] + return ["--n-pre", str(cfg["n_pre"]), "--n-post", str(cfg["n_post"])] + + +def dataset_for(estimator: str, scale: str) -> Path: + spec = BENCH_SPECS[estimator] + if spec["dataset"] == "mpdta": + return rc.ensure_mpdta() + cfg = SCALE_CONFIGS[scale][spec["dataset"]] + return rc.ensure_dataset(spec["dataset"], scale, cfg) + + +def run_cell( + estimator: str, + scale: str, + n_reps: int, + smoke: bool, +) -> Dict[str, Any]: + spec = BENCH_SPECS[estimator] + data_path = dataset_for(estimator, scale) + timeouts = cell_timeouts(estimator, scale) + py_extra = list(spec.get("py_extra", [])) + r_extra = list(spec.get("r_extra", [])) + if estimator == "multiperiod": + py_extra += multiperiod_extra(scale) + r_extra += multiperiod_extra(scale) + + reps = n_reps + if scale in spec.get("slow_scales", set()) and not smoke: + reps = min(n_reps, 4) + + cell: Dict[str, Any] = { + "estimator": spec["display"], + "scale": scale, + "data_file": data_path.name, + "data_sha256": rc.sha256_file(data_path), + "n_reps_requested": reps, + "captured_at_utc": rc.utc_now(), + "loadavg_at_cell": list(os.getloadavg()), + "flags": [], + } + print(f"\n=== {spec['display']} @ {scale} " f"(reps={reps}, data={data_path.name}) ===") + + arm_results: Dict[str, Dict[str, Any]] = {} + for arm in ARMS: + raw_out = rc.RAW_DIR / f"{estimator}_{scale}_{arm}.json" + if arm == "r": + + def _r_rep(raw_out=raw_out): + return rc.run_r_rep( + spec["r_script"], + data_path, + raw_out, + extra_args=r_extra, + timeout=timeouts["r"], + ) + + rep_fn = _r_rep + else: + backend = "python" if arm == "python_pure" else "rust" + + def _py_rep(backend=backend, raw_out=raw_out): + res = rc.run_python_rep( + spec["py_script"], + data_path, + raw_out, + backend=backend, + venv=HEADLINE_VENV, + extra_args=py_extra, + timeout=timeouts["python"], + ) + rc.validate_python_provenance(res, HEADLINE_PIN, HEADLINE_VENV, backend) + return res + + rep_fn = _py_rep + + arm_out = rc.run_arm( + f"{estimator}/{scale}/{arm}", + rep_fn, + reps, + allow_cv_rerun=not smoke, + ) + rc.validate_estimator_field( + arm_out["result"], + spec["r_estimator"] if arm == "r" else spec["py_estimator"], + ) + # compare_estimates reads timing from result["timing"]["stats"]. + arm_out["result"]["timing"] = arm_out["timing"] + arm_results[arm] = arm_out + cell["flags"].extend(f"{arm}:{fl}" for fl in arm_out["flags"]) + rc.cooldown(1 if smoke else rc.COOLDOWN_SECONDS) + + # Parity: rust arm vs R (matching the historical "python" arm choice). + comparison = compare_estimates( + arm_results["python_rust"]["result"], + arm_results["r"]["result"], + spec["display"], + se_rtol=spec["se_rtol"], + scale=scale, + python_pure_results=arm_results["python_pure"]["result"], + python_rust_results=arm_results["python_rust"]["result"], + ) + print( + f" parity: {'PASS' if comparison.passed else 'FAIL'} " + f"(ATT diff {comparison.att_diff:.2e}, " + f"SE rel {comparison.se_rel_diff:.2%})" + ) + if not comparison.passed: + cell["flags"].append("parity_fail") + # CI overlap (published criterion) is gated PER RENDERED ARM inside + # headline_gate_flags() below - both python_pure and python_rust. + + # Detailed parity surfaces (docs claim per-period / per-(g,t) / + # event-time / group effects match R - gate them, not just the headline). + detail: Dict[str, Any] = {} + for surface, join_keys in DETAIL_SURFACES.get(estimator, []): + for py_arm in ("python_pure", "python_rust"): + py_rows = arm_results[py_arm]["result"].get(surface) + r_rows = arm_results["r"]["result"].get(surface) + metrics = rc.compare_effect_arrays(py_rows, r_rows, join_keys, se_rtol=spec["se_rtol"]) + detail[f"{surface}:{py_arm}"] = metrics + tag = f"{surface}:{py_arm}" + if not metrics["keys_match"]: + cell["flags"].append(f"detail_keys_mismatch:{tag}") + if not metrics["att_ok"]: + cell["flags"].append(f"detail_att_gate_fail:{tag}:{metrics['max_att_diff']}") + if not metrics["se_ok"]: + cell["flags"].append(f"detail_se_gate_fail:{tag}:{metrics['max_se_rel_diff']}") + if detail: + cell["detail_comparison"] = detail + n_bad = sum(1 for fl in cell["flags"] if fl.startswith("detail_")) + print( + f" detail surfaces: {len(detail)} compared, " + f"{'ALL OK' if n_bad == 0 else f'{n_bad} GATE FAILURES'}" + ) + + # Strict headline gates: BOTH rendered Python arms must independently + # satisfy the documented ATT/SE tolerances vs R. Unlike + # compare_estimates(), CI overlap is NOT an escape hatch here. + def _vals(arm: str): + res = arm_results[arm]["result"] + return ( + res.get("overall_att", res.get("att")), + res.get("overall_se", res.get("se")), + ) + + r_att_v, r_se_v = _vals("r") + for py_arm in ("python_pure", "python_rust"): + py_att_v, py_se_v = _vals(py_arm) + cell["flags"].extend( + rc.headline_gate_flags( + py_arm, + py_att_v, + py_se_v, + r_att_v, + r_se_v, + att_atol=1e-4, + se_rtol=spec["se_gate_rtol"], + ) + ) + + # Pure-vs-rust backend agreement on the point estimate (hard gate). + def _att(arm: str) -> float: + res = arm_results[arm]["result"] + return float(res.get("overall_att", res.get("att"))) + + pure_rust_gap = abs(_att("python_pure") - _att("python_rust")) + if pure_rust_gap > 1e-8: + cell["flags"].append(f"pure_rust_att_gate_fail:{pure_rust_gap:.2e}") + + # Known-answer gates. + if estimator == "mpdta": + for arm in ARMS: + res = arm_results[arm]["result"] + att = float(res.get("overall_att", res.get("att"))) + if abs(att - rc.MPDTA_KNOWN_ATT) > 1e-5: + cell["flags"].append(f"mpdta_known_answer_fail:{arm}:{att:.6f}") + if estimator == "synthdid": + gap = abs(_att("python_rust") - _att("r")) + if gap > 1e-8: + # Same Frank-Wolfe algorithm on both sides - this is exactly the + # check that catches the invalidated-table class of bug. + cell["flags"].append(f"sdid_att_gate_fail:{gap:.2e}") + # Weight-vector gates with auditable committed metrics (raw vectors + # are slimmed out). Both omega and lambda reproduce R at machine + # precision once aligned by unit/period id (Python's accessor sorts + # by descending weight; R emits panel order - the id-keyed + # comparison makes ordering irrelevant). Tight 1e-8 gates on both. + weight_detail: Dict[str, Any] = cell.setdefault("detail_comparison", {}) + r_res = arm_results["r"]["result"] + for py_arm in ("python_pure", "python_rust"): + py_res = arm_results[py_arm]["result"] + for surface in ("unit_weights", "time_weights"): + id_field = surface[:-1] + "_ids" + metrics = rc.compare_weight_vectors( + py_res.get(surface), + r_res.get(surface), + atol=1e-8, + py_ids=py_res.get(id_field), + r_ids=r_res.get(id_field), + # Documented contract: the publication gate is id-aligned + # and may never silently degrade to positional order. + require_ids=True, + ) + weight_detail[f"{surface}:{py_arm}"] = metrics + if not metrics["ok"]: + cell["flags"].append( + f"sdid_weights_gate_fail:{surface}:{py_arm}:" f"{metrics['max_abs_diff']}" + ) + + for arm in ARMS: + entry = dict(arm_results[arm]) + entry["result"] = rc.slim_result(entry["result"]) + cell[arm] = entry + cell["comparison"] = asdict(comparison) + return cell + + +def collect_hard_failures(payload: Dict[str, Any]) -> List[str]: + """ + Hard failures across EVERY cell in the payload - including stale cells + preserved by merge-on-write during --only reruns. The runner's exit + status must reflect the artifact that would be published, not merely + the cells executed in this invocation. + """ + failures = [] + for key in sorted(payload.get("cells", {})): + hard = rc.hard_flags(payload["cells"][key].get("flags", [])) + if hard: + failures.append(f"{key}: {hard}") + return failures + + +def main() -> int: + ap = argparse.ArgumentParser(description=__doc__) + ap.add_argument( + "--setup", action="store_true", help="Create/refresh venvs + preflight, then exit" + ) + ap.add_argument("--smoke", action="store_true", help="Small-scale pipeline validation (reps=2)") + ap.add_argument( + "--only", + action="append", + default=None, + metavar="ESTIMATOR", + help="Limit to one estimator (repeatable): " f"{', '.join(BENCH_SPECS)}", + ) + ap.add_argument( + "--reps", + type=int, + default=8, + help="Subprocess reps per fast arm (first is excluded " + "from stats; slow SDID cells cap at 4)", + ) + ap.add_argument("--force", action="store_true", help="Override the load-average guard") + args = ap.parse_args() + + rc.RAW_DIR.mkdir(parents=True, exist_ok=True) + + if args.setup: + rc.setup_venv(HEADLINE_VENV, HEADLINE_PIN) + rc.preflight_venv(HEADLINE_VENV, HEADLINE_PIN) + rc.data_determinism_check() + return 0 + + # Preflight (cheap, always). + preflight_info = rc.preflight_venv(HEADLINE_VENV, HEADLINE_PIN) + rc.data_determinism_check() + rc.check_load(args.force, enforce=not args.smoke) + + estimators = args.only or list(BENCH_SPECS) + for est in estimators: + if est not in BENCH_SPECS: + ap.error(f"unknown estimator {est!r}") + + n_reps = 2 if args.smoke else args.reps + results_path = SMOKE_RESULTS_PATH if args.smoke else RESULTS_PATH + + # Merge-on-write so --only reruns update a single artifact. + payload: Dict[str, Any] = {"cells": {}} + if results_path.exists(): + with open(results_path) as f: + payload = json.load(f) + payload["run_metadata"] = rc.collect_run_metadata() + payload["headline_pin"] = HEADLINE_PIN + payload["smoke"] = args.smoke + fingerprint = rc.env_fingerprint( + payload["run_metadata"], HEADLINE_PIN, python_env=preflight_info + ) + payload["env_fingerprint"] = fingerprint + + t0 = time.time() + for est in estimators: + scales = ( + ["small"] + if args.smoke and "small" in BENCH_SPECS[est]["scales"] + else BENCH_SPECS[est]["scales"] + ) + for scale in scales: + cell = run_cell(est, scale, n_reps, args.smoke) + cell["env_fingerprint"] = fingerprint + payload["cells"][f"{est}/{scale}"] = cell + hard = rc.hard_flags(cell["flags"]) + if hard: + print(f" HARD GATE FAILURE {est}/{scale}: {hard}") + with open(results_path, "w") as f: + json.dump(payload, f, indent=2) + rc.cooldown(1 if args.smoke else rc.COOLDOWN_SECONDS) + + print(f"\nDone in {time.time() - t0:.0f}s -> {results_path}") + stale = sorted( + key for key, c in payload["cells"].items() if c.get("env_fingerprint") != fingerprint + ) + if stale: + print( + "\nARTIFACT NOT PUBLICATION-READY: cells captured under a " + f"different environment/protocol fingerprint: {stale}. " + "Re-run them under the current environment." + ) + failures = collect_hard_failures(payload) + if failures: + print( + "\nHARD-GATE FAILURES in the artifact (including any stale " + "merge-on-write cells; publication blocked until resolved):" + ) + for f_ in failures: + print(f" - {f_}") + return 1 + return 1 if stale else 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/benchmarks/refresh_2026_07/run_version_story.py b/benchmarks/refresh_2026_07/run_version_story.py new file mode 100644 index 00000000..3a52de52 --- /dev/null +++ b/benchmarks/refresh_2026_07/run_version_story.py @@ -0,0 +1,380 @@ +#!/usr/bin/env python3 +""" +Internal "optimization story" benchmark: released diff-diff 3.5.3 vs 3.7.0. + +Quantifies the June-July 2026 optimization arc (Rust demean_map kernel, +solve_ols marshalling slimming, CallawaySantAnna O(n_units) aggregation, +bootstrap memory tiling) by running the SAME benchmark scripts on the SAME +seed-42 data against two pinned PyPI wheels in isolated venvs, strictly +sequentially. v3.5.3 (2026-06-25) predates the entire 3.6.x wave. + +Both arms run at wheel defaults (backend auto -> Rust + Accelerate; no +DIFF_DIFF_* env knobs set) - i.e. exactly what `pip install diff-diff` users +got on each date. Note the opt-in solve_ols Cholesky fast path (#670) is in +NO released wheel and plays no role here. + +Methods: MultiPeriodDiD (absorbed-FE event study - exercises the demeaning +arc) and CallawaySantAnna (exercises the CS scaling arc) at the harness 20k +scale plus one larger story-scale cell each (~2M rows); BasicDiD as a cheap +"already fast, unchanged" context cell. SyntheticDiD is deliberately +excluded: its Frank-Wolfe rework shipped in v3.3.0, before both pins. + +This is a repo-internal artifact (results/version_story.{json,md}); it is +NOT part of the docs site. + +Usage: + python run_version_story.py --setup # venvs + preflight only + python run_version_story.py --smoke # small-scale compat check + python run_version_story.py # full run (idle machine!) + python run_version_story.py --report # regenerate md from json +""" + +import argparse +import json +import sys +import time +from pathlib import Path +from typing import Any, Dict, List + +import refresh_common as rc + +sys.path.insert(0, str(rc.REPO_ROOT)) + +from benchmarks.run_benchmarks import SCALE_CONFIGS # noqa: E402 + +PINS = {"dd353": "3.5.3", "dd370": "3.7.0"} +OLD_VENV, NEW_VENV = "dd353", "dd370" + +RESULTS_JSON = rc.RESULTS_DIR / "version_story.json" +RESULTS_MD = rc.RESULTS_DIR / "version_story.md" + +# Story-scale configs (~2M rows). Calibrate with --story-scale story50k +# first if unsure the old arm stays under ~10 min/rep. +STORY_CONFIGS: Dict[str, Dict[str, Dict[str, int]]] = { + "story2m": { + "multiperiod": {"n_units": 100_000, "n_pre": 10, "n_post": 10}, + "staggered": {"n_units": 100_000, "n_periods": 20, "n_cohorts": 8}, + }, + "story50k": { + "multiperiod": {"n_units": 50_000, "n_pre": 10, "n_post": 10}, + "staggered": {"n_units": 50_000, "n_periods": 20, "n_cohorts": 8}, + }, +} + +CELLS: List[Dict[str, Any]] = [ + { + "estimator": "basic", + "display": "BasicDiD (interaction OLS)", + "py_script": "benchmark_basic.py", + "py_estimator": "diff_diff.DifferenceInDifferences", + "dataset": "basic", + "scale": "20k", + "extra": ["--type", "basic"], + "reps": 6, + "timeout": 2400, + }, + { + "estimator": "multiperiod", + "display": "MultiPeriodDiD (absorbed FE)", + "py_script": "benchmark_multiperiod.py", + "py_estimator": "diff_diff.MultiPeriodDiD", + "dataset": "multiperiod", + "scale": "20k", + "reps": 6, + "timeout": 2400, + }, + { + "estimator": "callaway", + "display": "CallawaySantAnna", + "py_script": "benchmark_callaway.py", + "py_estimator": "diff_diff.CallawaySantAnna", + "dataset": "staggered", + "scale": "20k", + "reps": 6, + "timeout": 2400, + }, + { + "estimator": "multiperiod", + "display": "MultiPeriodDiD (absorbed FE)", + "py_script": "benchmark_multiperiod.py", + "py_estimator": "diff_diff.MultiPeriodDiD", + "dataset": "multiperiod", + "scale": "STORY", + "reps": 4, + "timeout": 3600, + }, + { + "estimator": "callaway", + "display": "CallawaySantAnna", + "py_script": "benchmark_callaway.py", + "py_estimator": "diff_diff.CallawaySantAnna", + "dataset": "staggered", + "scale": "STORY", + "reps": 4, + "timeout": 3600, + }, +] + + +def dataset_for(cell: Dict[str, Any], story_scale: str) -> Path: + scale = cell["scale"] + if scale == "STORY": + cfg = STORY_CONFIGS[story_scale][cell["dataset"]] + return rc.ensure_dataset(cell["dataset"], story_scale, cfg) + cfg = SCALE_CONFIGS[scale][cell["dataset"]] + return rc.ensure_dataset(cell["dataset"], scale, cfg) + + +def cell_extra(cell: Dict[str, Any], story_scale: str) -> List[str]: + extra = list(cell.get("extra", [])) + if cell["estimator"] == "multiperiod": + scale = cell["scale"] + cfg = ( + STORY_CONFIGS[story_scale]["multiperiod"] + if scale == "STORY" + else SCALE_CONFIGS[scale]["multiperiod"] + ) + extra += ["--n-pre", str(cfg["n_pre"]), "--n-post", str(cfg["n_post"])] + return extra + + +def run_story_cell( + cell: Dict[str, Any], + story_scale: str, + n_reps: int, + smoke: bool, +) -> Dict[str, Any]: + data_path = dataset_for(cell, story_scale) + extra = cell_extra(cell, story_scale) + scale_label = story_scale if cell["scale"] == "STORY" else cell["scale"] + out: Dict[str, Any] = { + "estimator": cell["display"], + "scale": scale_label, + "data_file": data_path.name, + "data_sha256": rc.sha256_file(data_path), + "flags": [], + } + print( + f"\n=== STORY {cell['display']} @ {scale_label} " + f"(reps={n_reps}, data={data_path.name}) ===" + ) + + for venv in (OLD_VENV, NEW_VENV): + pin = PINS[venv] + raw_out = rc.RAW_DIR / (f"story_{cell['estimator']}_{scale_label}_{venv}.json") + + def rep_fn(venv=venv, pin=pin, raw_out=raw_out): + res = rc.run_python_rep( + cell["py_script"], + data_path, + raw_out, + backend="auto", + venv=venv, + extra_args=extra, + timeout=cell["timeout"], + ) + # Wheel default = auto = rust backend on both pins. + rc.validate_python_provenance(res, pin, venv, "rust") + return res + + try: + arm = rc.run_arm( + f"{cell['estimator']}/{scale_label}/{pin}", + rep_fn, + n_reps, + allow_cv_rerun=not smoke, + ) + except Exception as exc: # old-wheel API drift: drop, never hack + print(f" [{pin}] CELL DROPPED: {exc}") + out["flags"].append(f"{pin}:dropped:{type(exc).__name__}") + out[pin] = {"dropped": True, "reason": str(exc)[:500]} + continue + rc.validate_estimator_field(arm["result"], cell["py_estimator"]) + arm["result"] = rc.slim_result(arm["result"]) + out[pin] = arm + out["flags"].extend(f"{pin}:{fl}" for fl in arm["flags"]) + rc.cooldown(1 if smoke else rc.COOLDOWN_SECONDS) + + # Cross-version invariants (balanced fixtures: estimates must agree). + old, new = out.get(PINS[OLD_VENV]), out.get(PINS[NEW_VENV]) + if old and new and not old.get("dropped") and not new.get("dropped"): + + def _att(entry): + res = entry["result"] + return float(res.get("overall_att", res.get("att"))) + + gap = abs(_att(old) - _att(new)) + out["att_gap_old_new"] = gap + if gap > 1e-6: + out["flags"].append(f"cross_version_att_shift:{gap:.2e}") + for lib in ("numpy_version", "pandas_version"): + vo = old["result"].get("provenance", {}).get(lib) + vn = new["result"].get("provenance", {}).get(lib) + if vo != vn: + out["flags"].append(f"venv_dep_mismatch:{lib}:{vo}!={vn}") + return out + + +def write_report(payload: Dict[str, Any]) -> None: + meta = payload["run_metadata"] + flagged = {k: c["flags"] for k, c in payload["cells"].items() if c.get("flags")} + lines = [ + "# diff-diff optimization story: released 3.5.3 vs 3.7.0", + "", + ] + if flagged: + lines += [ + "> **ATTENTION - flagged cells below; do not quote these numbers", + "> without resolving the flags:**", + ">", + ] + for k, fls in flagged.items(): + lines.append(f"> - `{k}`: {', '.join(fls)}") + lines.append("") + lines += [ + "Internal benchmark artifact (NOT part of the docs site). Same", + "benchmark scripts, same seed-42 data, two pinned PyPI wheels in", + "isolated venvs, wheel defaults (Rust + Accelerate backend, no", + "DIFF_DIFF_* env knobs), fresh subprocess per replication, strictly", + "sequential, untimed in-process warm-up, median of counted reps", + "(first rep excluded).", + "", + "v3.5.3 (2026-06-25) predates the June-July optimization arc;", + "v3.7.0 (2026-07-08) contains it (Rust demean_map FE-absorption", + "kernel, solve_ols marshalling slimming, CallawaySantAnna O(n_units)", + "aggregation + fused bootstrap, multiplier-bootstrap memory tiling).", + "The opt-in solve_ols Cholesky fast path (#670) is in NO released", + "wheel and plays no role in these numbers. SyntheticDiD is excluded:", + "its Frank-Wolfe rework shipped in v3.3.0, before both pins.", + "", + "| Estimator | Scale | Rows | 3.5.3 median (s) | 3.7.0 median (s) " "| Speedup |", + "|---|---|---|---|---|---|", + ] + for key, cell in payload["cells"].items(): + old = cell.get("3.5.3", {}) + new = cell.get("3.7.0", {}) + if old.get("dropped") or new.get("dropped") or not old or not new: + lines.append( + f"| {cell['estimator']} | {cell['scale']} | - | - | - | " + f"DROPPED ({'; '.join(cell['flags'])}) |" + ) + continue + t_old = old["timing"]["stats"]["median"] + t_new = new["timing"]["stats"]["median"] + rows = old["result"].get("metadata", {}).get("n_obs", "-") + speed = t_old / t_new if t_new else float("nan") + lines.append( + f"| {cell['estimator']} | {cell['scale']} | {rows:,} " + f"| {t_old:.3f} | {t_new:.3f} | **{speed:.2f}x** |" + if isinstance(rows, int) + else f"| {cell['estimator']} | {cell['scale']} | {rows} " + f"| {t_old:.3f} | {t_new:.3f} | **{speed:.2f}x** |" + ) + lines += [ + "", + "## Environment", + "", + f"- Captured: {meta['date_utc']}", + f"- Hardware: {meta['hardware']['cpu']}, " + f"{meta['hardware']['memory_gb']} GB, {meta['hardware']['os']} " + f"({meta['hardware']['arch']})", + f"- Wheels: diff-diff 3.5.3 and 3.7.0 (PyPI macosx-arm64, Rust " + f"backend + Apple Accelerate), python {meta['orchestrator_python']}", + f"- Protocol: {meta['protocol']}", + f"- Thread policy: {meta['thread_policy']}", + "", + "Flags (if any) per cell:", + "", + ] + for key, cell in payload["cells"].items(): + if cell.get("flags"): + lines.append(f"- `{key}`: {', '.join(cell['flags'])}") + if not any(c.get("flags") for c in payload["cells"].values()): + lines.append("- none") + lines.append("") + RESULTS_MD.write_text("\n".join(lines)) + print(f"wrote {RESULTS_MD}") + + +def main() -> int: + ap = argparse.ArgumentParser(description=__doc__) + ap.add_argument("--setup", action="store_true") + ap.add_argument( + "--smoke", action="store_true", help="small-scale old-wheel compat check (reps=2)" + ) + ap.add_argument( + "--report", action="store_true", help="regenerate version_story.md from existing json" + ) + ap.add_argument( + "--story-scale", + default="story2m", + choices=list(STORY_CONFIGS) + ["none"], + help="size of the large story cells (default story2m; " + "use story50k to calibrate, none to skip them)", + ) + ap.add_argument("--force", action="store_true") + args = ap.parse_args() + + rc.RAW_DIR.mkdir(parents=True, exist_ok=True) + + if args.report: + with open(RESULTS_JSON) as f: + write_report(json.load(f)) + return 0 + + if args.setup: + for venv, pin in PINS.items(): + rc.setup_venv(venv, pin) + rc.preflight_venv(venv, pin) + return 0 + + for venv, pin in PINS.items(): + rc.preflight_venv(venv, pin) + rc.check_load(args.force, enforce=not args.smoke) + + if args.smoke: + # Compat guard: every story script x the OLD wheel at small scale. + cells = [ + {**c, "scale": "small", "reps": 2, "timeout": 600} + for c in CELLS + if c["scale"] == "20k" # one per script, small data + ] + else: + cells = [c for c in CELLS if not (c["scale"] == "STORY" and args.story_scale == "none")] + + payload: Dict[str, Any] = {"cells": {}} + if RESULTS_JSON.exists() and not args.smoke: + with open(RESULTS_JSON) as f: + payload = json.load(f) + payload["run_metadata"] = rc.collect_run_metadata() + payload["pins"] = PINS + + t0 = time.time() + for cell in cells: + scale_label = args.story_scale if cell["scale"] == "STORY" else cell["scale"] + n_reps = 2 if args.smoke else cell["reps"] + out = run_story_cell(cell, args.story_scale, n_reps, args.smoke) + payload["cells"][f"{cell['estimator']}/{scale_label}"] = out + if not args.smoke: + with open(RESULTS_JSON, "w") as f: + json.dump(payload, f, indent=2) + rc.cooldown(1 if args.smoke else rc.COOLDOWN_SECONDS) + + print(f"\nStory run done in {time.time() - t0:.0f}s") + flagged = {k: c["flags"] for k, c in payload["cells"].items() if c.get("flags")} + if not args.smoke: + with open(RESULTS_JSON, "w") as f: + json.dump(payload, f, indent=2) + write_report(payload) + else: + print("(smoke mode: no artifacts written)") + if flagged: + print("\nFLAGGED CELLS (report banner added; exit nonzero):") + for k, fls in flagged.items(): + print(f" - {k}: {fls}") + return 1 + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/benchmarks/run_benchmarks.py b/benchmarks/run_benchmarks.py index bee27a0f..75894616 100644 --- a/benchmarks/run_benchmarks.py +++ b/benchmarks/run_benchmarks.py @@ -613,9 +613,9 @@ def run_basic_did_benchmark( n_replications: int = 1, backends: Optional[List[str]] = None, ) -> Dict[str, Any]: - """Run basic DiD / TWFE benchmarks (Python and R) with replications.""" + """Run BasicDiD (interaction OLS) benchmarks (Python and R) with replications.""" print(f"\n{'='*60}") - print(f"BASIC DID / TWFE BENCHMARK ({scale})") + print(f"BASIC DID BENCHMARK ({scale})") print(f"{'='*60}") if backends is None: @@ -650,7 +650,7 @@ def run_basic_did_benchmark( "benchmark_basic.py", data_path, py_output, - extra_args=["--type", "twfe"], + extra_args=["--type", "basic"], timeout=timeouts["python"], backend=backend, ) @@ -688,7 +688,7 @@ def run_basic_did_benchmark( "benchmark_fixest.R", data_path, r_output, - extra_args=["--type", "twfe"], + extra_args=["--type", "basic"], timeout=timeouts["r"], ) r_timings.append(r_result["timing"]["total_seconds"]) @@ -713,7 +713,7 @@ def run_basic_did_benchmark( comparison = compare_estimates( results["python"], results["r"], - "BasicDiD/TWFE", + "BasicDiD", scale=scale, python_pure_results=results.get("python_pure"), python_rust_results=results.get("python_rust"), diff --git a/docs/benchmarks.rst b/docs/benchmarks.rst index 3ad3b7ff..0388bd70 100644 --- a/docs/benchmarks.rst +++ b/docs/benchmarks.rst @@ -58,7 +58,12 @@ Tolerance Thresholds ~~~~~~~~~~~~~~~~~~~~ - **Point estimates (ATT)**: Absolute difference < 1e-4 or relative < 1% -- **Standard errors**: Relative difference < 10% +- **Standard errors**: Relative difference < 10%. Exception: SyntheticDiD + placebo SEs use a documented 35% Monte Carlo-bounded gate - both + implementations estimate the placebo variance by simulation and R's + placebo permutation is unseeded, so SEs agree in distribution rather + than draw-by-draw (see the SyntheticDiD methodology registry note; + the deterministic Frank-Wolfe ATT is gated at 1e-8) - **Confidence intervals**: Must overlap Benchmark Results @@ -67,6 +72,8 @@ Benchmark Results Summary Table ~~~~~~~~~~~~~ +.. refresh-table-start: summary + .. list-table:: :header-rows: 1 :widths: 25 20 20 15 20 @@ -97,11 +104,15 @@ Summary Table - Yes - **PASS** +.. refresh-table-end: summary + Basic DiD Results ~~~~~~~~~~~~~~~~~ **Data**: 100 units, 4 periods, true ATT = 5.0 (small scale) +.. refresh-table-start: accuracy_basic + .. list-table:: :header-rows: 1 @@ -126,6 +137,8 @@ Basic DiD Results - 0.041 - **22x faster** +.. refresh-table-end: accuracy_basic + **Validation**: PASS - Results are numerically identical across all implementations. MultiPeriodDiD Results @@ -133,6 +146,8 @@ MultiPeriodDiD Results **Data**: 200 units, 8 periods (4 pre, 4 post), true ATT = 3.0 (small scale) +.. refresh-table-start: accuracy_multiperiod + .. list-table:: :header-rows: 1 @@ -162,6 +177,8 @@ MultiPeriodDiD Results - 0.035 - **7x faster** (pure) +.. refresh-table-end: accuracy_multiperiod + **Validation**: PASS - Both average ATT and all period-level effects match R's ``fixest::feols(outcome ~ treated * time_f | unit)`` to machine precision. The regression includes unit fixed effects (absorbed via ``| unit`` in R, within- @@ -173,6 +190,8 @@ Synthetic DiD Results **Data**: 50 units (40 control, 10 treated), 20 periods, true ATT = 4.0 +.. refresh-table-start: accuracy_synthdid + .. list-table:: :header-rows: 1 @@ -197,11 +216,15 @@ Synthetic DiD Results - 8.19 - **2.4x faster** (pure) +.. refresh-table-end: accuracy_synthdid + **Validation**: PASS - ATT estimates are numerically identical across all implementations. Both diff-diff and R's synthdid use Frank-Wolfe optimization with two-pass sparsification and auto-computed regularization (``zeta_omega``, -``zeta_lambda``), producing identical unit and time weights. Both use -placebo-based variance estimation (Algorithm 4 from Arkhangelsky et al. 2021). +``zeta_lambda``), producing identical unit and time weights (reproduced at +< 1e-8 by the benchmark harness's id-aligned per-unit comparison; see the +SyntheticDiD methodology registry note). Both use placebo-based variance +estimation (Algorithm 4 from Arkhangelsky et al. 2021). The small SE difference (0.3% at small scale, up to ~7% at larger scales) is due to Monte Carlo variance in the placebo procedure, which randomly permutes @@ -213,6 +236,8 @@ Callaway-Sant'Anna Results **Data**: 200 units, 8 periods, 3 treatment cohorts, dynamic effects (small scale) +.. refresh-table-start: accuracy_callaway + .. list-table:: :header-rows: 1 @@ -237,6 +262,8 @@ Callaway-Sant'Anna Results - 0.070 ± 0.001 - **10x faster** +.. refresh-table-end: accuracy_callaway + **Validation**: PASS - Both point estimates and standard errors match R exactly. **Key findings from investigation:** @@ -271,6 +298,10 @@ implementations: - **Python (Pure)**: diff-diff with NumPy/SciPy only (no Rust backend) - **Python (Rust)**: diff-diff with optional Rust backend enabled +.. refresh-table-start: environment + +.. refresh-table-end: environment + .. note:: **v2.0.0 Rust Backend**: diff-diff v2.0.0 introduces an optional Rust backend @@ -291,6 +322,8 @@ Three-Way Performance Summary **BasicDiD/TWFE Results:** +.. refresh-table-start: perf_basic + .. list-table:: :header-rows: 1 :widths: 12 15 18 18 12 12 @@ -332,8 +365,12 @@ Three-Way Performance Summary - **2x** - 0.9x +.. refresh-table-end: perf_basic + **CallawaySantAnna Results:** +.. refresh-table-start: perf_callaway + .. list-table:: :header-rows: 1 :widths: 12 15 18 18 12 12 @@ -375,8 +412,12 @@ Three-Way Performance Summary - **4x** - 1.0x +.. refresh-table-end: perf_callaway + **SyntheticDiD Results:** +.. refresh-table-start: perf_synthdid + .. list-table:: :header-rows: 1 :widths: 12 15 18 18 12 12 @@ -406,6 +447,8 @@ Three-Way Performance Summary - **16.5x** - 0.1x +.. refresh-table-end: perf_synthdid + .. note:: **SyntheticDiD Performance**: diff-diff's pure Python backend achieves @@ -556,6 +599,8 @@ minimum wage policy changes: Results Comparison ~~~~~~~~~~~~~~~~~~ +.. refresh-table-start: mpdta + .. list-table:: :header-rows: 1 :widths: 25 25 25 25 @@ -577,6 +622,8 @@ Results Comparison - 0.039s ± 0.006s - **14.4x faster** +.. refresh-table-end: mpdta + **Key Findings:** 1. **Point estimates match exactly**: The overall ATT of -0.039951 is identical @@ -913,8 +960,9 @@ When to Trust Results ATT and all period-level effects match to machine precision. Use with confidence. - **SyntheticDiD**: Point estimates are numerically identical (< 1e-10 diff) and - standard errors match closely (0.3% diff at small scale). Both implementations - use Frank-Wolfe optimization with identical weights. Use + standard errors agree within placebo Monte Carlo dispersion. Both + implementations use Frank-Wolfe optimization with identical unit and time + weights (verified by id-aligned comparison in the benchmark harness). Use ``variance_method="placebo"`` (default) to match R's inference. Results are fully validated. @@ -932,8 +980,10 @@ Known Differences 2. **Aggregation Weights**: Overall ATT is a weighted average of ATT(g,t). Weighting schemes may differ between implementations. -3. **Placebo Variance**: SyntheticDiD SE estimates differ slightly (0.3-7%) - across implementations due to Monte Carlo variance in the placebo procedure. - Point estimates and unit/time weights are numerically identical since both - implementations use the same Frank-Wolfe optimizer. +3. **Placebo Variance**: SyntheticDiD SE estimates differ across + implementations due to Monte Carlo variance in the placebo procedure + (R's placebo permutation is unseeded). Point estimates and unit/time + weights are numerically identical since both implementations use the + same Frank-Wolfe optimizer (weights verified by id-aligned comparison + in the benchmark harness). diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md index dd45ced7..e89890c1 100644 --- a/docs/methodology/REGISTRY.md +++ b/docs/methodology/REGISTRY.md @@ -2342,6 +2342,31 @@ Convergence criterion: stop when objective decrease < min_decrease² (default mi **Validation:** (a) hand-computed 2-stratum FPC magnitude regression (`test_jackknife_full_design_fpc_reduces_se_magnitude` — asserts `SE_fpc == SE_nofpc · sqrt(1 - f)` at `rtol=1e-10`), (b) self-consistency between the returned SE and the stratum-aggregation formula applied to the returned LOO estimates, (c) single-PSU-stratum skip, (d) all-strata-skipped UserWarning + NaN, (e) unstratified single-PSU short-circuit, (f) deterministic-dispatch regression. - **Note:** P-value computation is variance-method dependent. Placebo (Algorithm 4) uses the empirical null formula `max(mean(|variance_effects| ≥ |att|), 1/(r+1))` because permuting control indices generates draws from the null distribution (centered on 0). Bootstrap (Algorithm 2) and jackknife (Algorithm 3) use the analytical p-value from `safe_inference(att, se)` (normal-theory): bootstrap draws are centered on `τ̂` (sampling distribution of the estimator) and jackknife pseudo-values are not null draws, so the empirical null formula is invalid for them. This matches R's `synthdid::vcov()` convention, where variance is returned and inference is normal-theory from the SE. +- **Note (benchmark SE gate is Monte Carlo-bounded, 2026-07):** The public + benchmark refresh (`benchmarks/refresh_2026_07/`) compares placebo SEs + against R `synthdid` across INDEPENDENT placebo draws: R's + `vcov(method="placebo")` does not seed its permutation sequence, so the two + implementations agree draw-by-draw only when the permutation sequence is + threaded explicitly (verified at < 1e-8 SE tolerance in + `tests/test_methodology_sdid.py::TestJackknifeSERParity::test_placebo_se_matches_r`); + across independent draws the SEs agree in distribution, not point-by-point. + The benchmark harness therefore hard-gates SyntheticDiD SEs at a 35% + relative bound (vs the standard 10% used for analytical-SE estimators), + records R's rep-to-rep SE values in the committed results artifact so the + observed Monte Carlo dispersion is auditable, and gates the deterministic + Frank-Wolfe ATT at 1e-8. The generated benchmark summary table carries the + same footnote. **Weight-vector gates:** unit (ω) and time (λ) weight + vectors both reproduce R at machine precision (observed ≤ 5e-13 on the + benchmark small fixture) and are hard-gated at 1e-8 max abs diff, aligned + by unit/period id. The id alignment is load-bearing: + `SyntheticDiDResults.get_unit_weights_df()` returns weights sorted by + DESCENDING WEIGHT while R's `synthdid` emits panel control-row order, so a + naive positional comparison of the two vectors fabricates apparent ~1e-2 + divergence out of pure ordering. The benchmark scripts emit id-keyed + vectors on both sides and the comparator aligns on ids (fail-closed on key + mismatch/duplicates). Per-run max-abs-diff metrics for both vectors are + committed with the results. + - **Note (coverage Monte Carlo calibration):** `benchmarks/data/sdid_coverage.json` carries empirical rejection rates across the three variance methods on 4 representative null-panel DGPs (500 seeds × B=200, regenerable via `benchmarks/python/coverage_sdid.py`). The fourth DGP (`stratified_survey`, added in PR #355) validates the survey-bootstrap calibration; jackknife is also reported with a documented anti-conservatism caveat; placebo is N/A on this DGP because its cohort packs into a single stratum with 0 never-treated units (stratified-permutation allocator is structurally infeasible — see `test_placebo_full_design_raises_on_zero_control_stratum` / `_undersupplied_stratum` for the enforced behavior). Under H0 the nominal rejection rate at each α equals α; rates substantially above α indicate anti-conservatism, rates below indicate over-coverage. | DGP | method | α=0.01 | α=0.05 | α=0.10 | mean SE / true SD | diff --git a/tests/test_benchmark_refresh_gates.py b/tests/test_benchmark_refresh_gates.py new file mode 100644 index 00000000..73a0ab2d --- /dev/null +++ b/tests/test_benchmark_refresh_gates.py @@ -0,0 +1,510 @@ +""" +Regression tests for the 2026-07 benchmark-refresh publication gates. + +These lock the fail-closed behavior of the refresh harness in +benchmarks/refresh_2026_07/: out-of-tolerance headline SE parity must gate +even when confidence intervals overlap, BOTH rendered Python arms must +independently satisfy the tolerances vs R, and the table generator must +refuse any results payload carrying hard correctness flags. + +No R, venvs, or subprocesses involved - pure unit tests on the gate logic. +""" + +import sys +from pathlib import Path + +import pytest + +# The refresh harness is repo tooling, not wheel content: CI's Python-test +# jobs copy tests/ to a temp directory and run against the installed wheel, +# where benchmarks/ does not exist. Skip the whole module there (same +# convention as the benchmarks/data golden-file skips). +REFRESH_DIR = Path(__file__).parent.parent / "benchmarks" / "refresh_2026_07" +if not (REFRESH_DIR / "refresh_common.py").exists(): + pytest.skip( + "benchmark refresh harness not present (installed-wheel test layout)", + allow_module_level=True, + ) +sys.path.insert(0, str(REFRESH_DIR)) + +import gen_benchmark_tables as gen # noqa: E402 +import refresh_common as rc # noqa: E402 + + +class TestHardFlagClassification: + def test_all_hard_tokens_detected(self): + flags = [ + "parity_fail", + "mpdta_known_answer_fail:r:-0.05", + "headline_att_gate_fail:python_pure:1e-2", + "headline_se_gate_fail:python_rust:0.500", + "ci_gate_fail", + "sdid_weights_gate_fail:unit_weights:python_rust:0.01", + "detail_keys_mismatch:event_study:python_pure", + "sdid_att_gate_fail:1e-6", + "pure_rust_att_gate_fail:2e-7", + ] + assert rc.hard_flags(flags) == flags + + def test_soft_flags_ignored(self): + soft = ["python_rust:cv_flag:0.150"] + assert rc.hard_flags(soft) == [] + + def test_rep_att_spread_is_hard(self): + # Non-deterministic point estimates across seeded replications must + # block publication (rep_att_gate_fail carries the hard token). + assert rc.hard_flags(["python_rust:rep_att_gate_fail:1.0e-09"]) + + +class TestHeadlineGates: + """Strict ATT/SE gates - CI overlap must NOT be an escape hatch.""" + + def test_se_mismatch_gates_despite_ci_overlap(self): + # ATT identical, SE 30% off with wide CIs that overlap heavily: + # compare_estimates() would report passed=True via ci_overlap; the + # headline gate must still fail the SE. + flags = rc.headline_gate_flags( + "python_rust", + py_att=2.5, + py_se=0.13, + r_att=2.5, + r_se=0.10, + att_atol=1e-4, + se_rtol=0.10, + ) + assert any(fl.startswith("headline_se_gate_fail:python_rust") for fl in flags) + assert rc.hard_flags(flags), "SE gate flag must be a hard flag" + + def test_att_mismatch_gates(self): + flags = rc.headline_gate_flags( + "python_pure", + py_att=2.60, + py_se=0.10, + r_att=2.50, + r_se=0.10, + att_atol=1e-4, + se_rtol=0.10, + ) + assert any(fl.startswith("headline_att_gate_fail:python_pure") for fl in flags) + assert rc.hard_flags(flags) + + def test_nonfinite_values_gate(self): + flags = rc.headline_gate_flags( + "python_rust", + py_att=float("nan"), + py_se=0.1, + r_att=2.5, + r_se=0.1, + att_atol=1e-4, + se_rtol=0.10, + ) + assert any("nonfinite" in fl for fl in flags) + assert rc.hard_flags(flags) + flags_se = rc.headline_gate_flags( + "python_rust", + py_att=2.5, + py_se=float("inf"), + r_att=2.5, + r_se=0.1, + att_atol=1e-4, + se_rtol=0.10, + ) + assert any(fl.startswith("headline_se_gate_fail") for fl in flags_se) + + def test_missing_values_gate(self): + flags = rc.headline_gate_flags( + "python_pure", + py_att=None, + py_se=None, + r_att=2.5, + r_se=0.1, + att_atol=1e-4, + se_rtol=0.10, + ) + assert flags and rc.hard_flags(flags) + + def test_ci_overlap_gated_per_arm(self): + # Pure-arm CI disjoint from R must gate even though ATT/SE checks + # are relative: ATT far apart with tiny SEs -> no overlap. + flags = rc.headline_gate_flags( + "python_pure", + py_att=10.0, + py_se=0.01, + r_att=2.5, + r_se=0.01, + att_atol=1e-4, + se_rtol=0.10, + ) + assert any(fl == "ci_gate_fail:python_pure" for fl in flags) + assert rc.hard_flags(["ci_gate_fail:python_pure"]) + + def test_overlapping_cis_produce_no_ci_flag(self): + flags = rc.headline_gate_flags( + "python_rust", + py_att=2.5, + py_se=0.13, + r_att=2.5, + r_se=0.10, + att_atol=1e-4, + se_rtol=0.10, + ) + assert not any(fl.startswith("ci_gate_fail") for fl in flags) + + def test_within_tolerance_produces_no_flags(self): + flags = rc.headline_gate_flags( + "python_rust", + py_att=2.5000000001, + py_se=0.1000001, + r_att=2.5, + r_se=0.10, + att_atol=1e-4, + se_rtol=0.10, + ) + assert flags == [] + + +class TestCompareEffectArraysFailClosed: + PY = [{"group": 3, "time": 4, "att": 1.0, "se": 0.1}] + R = [{"group": 3, "time": 4, "att": 1.0, "se": 0.1}] + + def test_clean_arrays_pass(self): + m = rc.compare_effect_arrays(self.PY, self.R, ["group", "time"], se_rtol=0.1) + assert m["keys_match"] and m["att_ok"] and m["se_ok"] + + def test_duplicate_join_keys_break_keys_match(self): + py = self.PY + [{"group": 3, "time": 4, "att": 1.0, "se": 0.1}] + m = rc.compare_effect_arrays(py, self.R, ["group", "time"], se_rtol=0.1) + assert not m["keys_match"] + assert m["n_dup_python"] == 1 + + def test_nonfinite_att_breaks_keys_match(self): + py = self.PY + [{"group": 3, "time": 5, "att": float("nan"), "se": 0.1}] + r = self.R + [{"group": 3, "time": 5, "att": float("nan"), "se": 0.1}] + m = rc.compare_effect_arrays(py, r, ["group", "time"], se_rtol=0.1) + assert not m["keys_match"] + assert m["n_dropped_python"] == 1 and m["n_dropped_r"] == 1 + + def test_missing_se_breaks_se_ok(self): + py = [{"group": 3, "time": 4, "att": 1.0, "se": None}] + m = rc.compare_effect_arrays(py, self.R, ["group", "time"], se_rtol=0.1) + assert m["att_ok"] and not m["se_ok"] + assert m["n_se_compared"] == 0 + + +class TestEstimatorFieldGuard: + """A silently ignored --type flag can never publish mislabeled numbers.""" + + def test_mismatch_raises(self): + with pytest.raises(RuntimeError, match="estimator mismatch"): + rc.validate_estimator_field( + {"estimator": "diff_diff.DifferenceInDifferences"}, + "diff_diff.TwoWayFixedEffects", + ) + + def test_missing_field_raises(self): + with pytest.raises(RuntimeError, match="estimator mismatch"): + rc.validate_estimator_field({}, "diff_diff.CallawaySantAnna") + + def test_match_passes(self): + rc.validate_estimator_field( + {"estimator": "fixest::feols (absorbed FE)"}, + "fixest::feols (absorbed FE)", + ) + + def test_refresh_specs_declare_distinct_wiring(self): + # Every headline spec must pin the expected estimator fields, and the + # basic/twfe split must be wired to the correct scripts. + import run_refresh + + for name, spec in run_refresh.BENCH_SPECS.items(): + assert spec.get("py_estimator"), f"{name}: missing py_estimator" + assert spec.get("r_estimator"), f"{name}: missing r_estimator" + assert spec.get("se_gate_rtol"), f"{name}: missing se_gate_rtol" + basic = run_refresh.BENCH_SPECS["basic"] + twfe = run_refresh.BENCH_SPECS["twfe"] + assert basic["py_script"] == "benchmark_basic.py" + assert basic["py_estimator"] == "diff_diff.DifferenceInDifferences" + assert twfe["py_script"] == "benchmark_twfe.py" + assert twfe["r_script"] == "benchmark_twfe.R" + assert twfe["py_estimator"] == "diff_diff.TwoWayFixedEffects" + assert twfe["r_estimator"] == "fixest::feols (absorbed FE)" + + +class TestSyntheticDiDSEGateDocumented: + """The wider SDID SE gate is a documented Monte Carlo bound. + + Both implementations estimate the placebo variance by Monte Carlo and + R's placebo permutation is unseeded, so SEs agree in distribution, not + draw-by-draw - see the REGISTRY.md SyntheticDiD note "benchmark SE gate + is Monte Carlo-bounded (2026-07)". + """ + + def test_gate_value_matches_registry_note(self): + import run_refresh + + assert run_refresh.BENCH_SPECS["synthdid"]["se_gate_rtol"] == 0.35 + + def test_within_mc_bound_passes_beyond_it_gates(self): + ok = rc.headline_gate_flags( + "python_rust", + py_att=3.84, + py_se=0.12, + r_att=3.84, + r_se=0.10, + att_atol=1e-4, + se_rtol=0.35, + ) + assert ok == [] # 20% is within the documented MC bound + bad = rc.headline_gate_flags( + "python_rust", + py_att=3.84, + py_se=0.145, + r_att=3.84, + r_se=0.10, + att_atol=1e-4, + se_rtol=0.35, + ) + assert any(fl.startswith("headline_se_gate_fail") for fl in bad) + + +class TestEnvironmentFingerprint: + META = { + "r_version": "R version 4.5.2", + "r_packages": {"did": "2.5.1"}, + "hardware": {"cpu": "Apple M4 Max", "os": "macOS 15"}, + "orchestrator_python": "3.14.4", + "protocol": "p", + "thread_policy": "t", + } + + def test_fingerprint_changes_with_environment(self): + fp1 = rc.env_fingerprint(self.META, "3.7.0") + meta2 = dict(self.META, r_packages={"did": "2.6.0"}) + assert fp1 != rc.env_fingerprint(meta2, "3.7.0") + assert fp1 != rc.env_fingerprint(self.META, "3.7.1") + assert fp1 == rc.env_fingerprint(dict(self.META), "3.7.0") + + def test_fingerprint_changes_with_python_arm_provenance(self): + py1 = {"numpy": "2.5.1", "pandas": "3.0.3", "python": "3.14.4"} + py2 = dict(py1, numpy="2.6.0") + fp1 = rc.env_fingerprint(self.META, "3.7.0", python_env=py1) + assert fp1 != rc.env_fingerprint(self.META, "3.7.0", python_env=py2) + assert fp1 == rc.env_fingerprint(self.META, "3.7.0", python_env=dict(py1)) + + def test_thread_env_vars_stripped_from_child_env(self, monkeypatch): + monkeypatch.setenv("OMP_NUM_THREADS", "1") + monkeypatch.setenv("RAYON_NUM_THREADS", "2") + env = rc._child_env() + assert "OMP_NUM_THREADS" not in env + assert "RAYON_NUM_THREADS" not in env + assert env["DIFF_DIFF_BENCH_USE_INSTALLED"] == "1" + + def test_slim_result_redacts_local_paths(self): + slim = rc.slim_result( + { + "provenance": { + "diff_diff_path": "/Users/someone/repo/venvs/dd370/lib/" + "python3.14/site-packages/diff_diff/__init__.py", + "python_executable": "/Users/someone/repo/venvs/dd370/bin/python", + } + } + ) + prov = slim["provenance"] + assert "someone" not in prov["diff_diff_path"] + assert prov["diff_diff_path"].startswith(".../") + assert "someone" not in prov["python_executable"] + + def test_generator_refuses_mixed_fingerprints(self): + payload = { + "env_fingerprint": "abc", + "cells": { + "basic/small": {"flags": [], "env_fingerprint": "abc"}, + "callaway/small": {"flags": [], "env_fingerprint": "OLD"}, + }, + } + with pytest.raises(SystemExit, match=r"different\s+environment"): + gen.assert_uniform_environment(payload) + + def test_generator_accepts_uniform_fingerprints(self): + payload = { + "env_fingerprint": "abc", + "cells": { + "basic/small": {"flags": [], "env_fingerprint": "abc"}, + }, + } + gen.assert_uniform_environment(payload) + + +class TestRunArmNonFiniteReplicationGate: + def test_nan_on_intermediate_rep_hard_flags(self): + # A NaN replication must gate even when the LAST rep (which feeds + # the headline gates) is finite. + results = iter( + [ + {"timing": {"total_seconds": 0.01}, "att": float("nan"), "se": 0.1}, + {"timing": {"total_seconds": 0.01}, "att": 2.5, "se": 0.1}, + ] + ) + out = rc.run_arm("test/nan", lambda: next(results), 2, allow_cv_rerun=False) + assert out["n_nonfinite_att"] == 1 + assert any(fl.startswith("rep_att_gate_fail:nonfinite") for fl in out["flags"]) + assert rc.hard_flags(out["flags"]) + + def test_nan_pass_with_high_cv_is_not_masked_by_rerun(self): + # First pass: NaN SE on rep 1 AND noisy timings (CV > 10%). The CV + # rerun must NOT discard the correctness failure by replacing the + # pass with a clean second run. + results = iter( + [ + {"timing": {"total_seconds": 0.01}, "att": 2.5, "se": float("nan")}, + {"timing": {"total_seconds": 0.01}, "att": 2.5, "se": 0.1}, + {"timing": {"total_seconds": 0.50}, "att": 2.5, "se": 0.1}, + # If a rerun were (wrongly) attempted, these would be consumed: + {"timing": {"total_seconds": 0.01}, "att": 2.5, "se": 0.1}, + {"timing": {"total_seconds": 0.01}, "att": 2.5, "se": 0.1}, + {"timing": {"total_seconds": 0.01}, "att": 2.5, "se": 0.1}, + ] + ) + consumed = {"n": 0} + + def rep_fn(): + consumed["n"] += 1 + return next(results) + + out = rc.run_arm("test/nan+cv", rep_fn, 3, allow_cv_rerun=True) + assert out["cv"] > rc.CV_FLAG # the rerun trigger was genuinely armed + assert consumed["n"] == 3, "rerun must be skipped when correctness flags exist" + assert any(fl.startswith("rep_se_gate_fail:nonfinite") for fl in out["flags"]) + assert rc.hard_flags(out["flags"]) + + def test_all_finite_reps_produce_no_nonfinite_flags(self): + results = iter( + [ + {"timing": {"total_seconds": 0.01}, "att": 2.5, "se": 0.1}, + {"timing": {"total_seconds": 0.01}, "att": 2.5, "se": 0.1}, + ] + ) + out = rc.run_arm("test/ok", lambda: next(results), 2, allow_cv_rerun=False) + assert out["n_nonfinite_att"] == 0 and out["n_nonfinite_se"] == 0 + assert not rc.hard_flags(out["flags"]) + + +class TestRunnerExitScansWholeArtifact: + def test_stale_flagged_cells_fail_the_run(self): + # merge-on-write can preserve older hard-flagged cells during --only + # reruns; the exit status must reflect the whole artifact. + import run_refresh + + payload = { + "cells": { + "basic/small": {"flags": []}, # fresh, clean + "synthdid/5k": {"flags": ["headline_se_gate_fail:python_pure:0.9"]}, + } + } + failures = run_refresh.collect_hard_failures(payload) + assert failures == ["synthdid/5k: ['headline_se_gate_fail:python_pure:0.9']"] + + def test_clean_artifact_no_failures(self): + import run_refresh + + assert run_refresh.collect_hard_failures({"cells": {"a/b": {"flags": []}}}) == [] + + +class TestWeightVectorComparison: + def test_identical_weights_pass(self): + m = rc.compare_weight_vectors([0.5, 0.5, 0.0], [0.5, 0.5, 0.0]) + assert m["ok"] and m["max_abs_diff"] == 0.0 + + def test_out_of_order_weights_align_by_ids(self): + # Python's get_unit_weights_df() sorts by DESCENDING WEIGHT while R + # emits panel order - identical weights must compare equal once + # aligned by unit id (a positional comparison would fabricate a + # 0.4 max diff here). + py_w, py_ids = [0.5, 0.3, 0.2], [7, 2, 5] # weight-sorted + r_w, r_ids = [0.3, 0.2, 0.5], [2, 5, 7] # panel order + m = rc.compare_weight_vectors(py_w, r_w, py_ids=py_ids, r_ids=r_ids) + assert m["aligned_by_ids"] and m["ok"] and m["max_abs_diff"] == 0.0 + + def test_id_alignment_catches_true_divergence(self): + m = rc.compare_weight_vectors( + [0.5, 0.3, 0.2], [0.3, 0.3, 0.4], py_ids=[7, 2, 5], r_ids=[2, 5, 7] + ) + assert m["aligned_by_ids"] and not m["ok"] + assert abs(m["max_abs_diff"] - 0.1) < 1e-15 + + def test_string_ids_from_r_align_with_int_ids(self): + # R rownames are character; Python ids are ints - keys normalize. + m = rc.compare_weight_vectors([0.6, 0.4], [0.4, 0.6], py_ids=[10, 3], r_ids=["3", "10"]) + assert m["aligned_by_ids"] and m["ok"] + + def test_missing_ids_fail_closed_when_required(self): + # The SDID publication gate documents id alignment; positionally + # identical weights must still FAIL if either side stops emitting + # ids (a script regression cannot silently weaken the contract). + m = rc.compare_weight_vectors([0.5, 0.5], [0.5, 0.5], require_ids=True) + assert not m["ok"] and not m["aligned_by_ids"] + m2 = rc.compare_weight_vectors( + [0.5, 0.5], [0.5, 0.5], py_ids=[1, 2], r_ids=None, require_ids=True + ) + assert not m2["ok"] + + def test_require_ids_passes_with_valid_ids(self): + m = rc.compare_weight_vectors( + [0.5, 0.5], [0.5, 0.5], py_ids=[1, 2], r_ids=[1, 2], require_ids=True + ) + assert m["ok"] and m["aligned_by_ids"] + + def test_positional_fallback_still_allowed_when_not_required(self): + m = rc.compare_weight_vectors([0.5, 0.5], [0.5, 0.5]) + assert m["ok"] and not m["aligned_by_ids"] + + def test_mismatched_id_sets_fail_closed(self): + m = rc.compare_weight_vectors([0.5, 0.5], [0.5, 0.5], py_ids=[1, 2], r_ids=[1, 3]) + assert not m["ok"] and m["max_abs_diff"] is None + + def test_duplicate_ids_fail_closed(self): + m = rc.compare_weight_vectors([0.5, 0.5], [0.5, 0.5], py_ids=[1, 1], r_ids=[1, 2]) + assert not m["ok"] + + def test_divergent_weights_fail(self): + m = rc.compare_weight_vectors([0.5, 0.5], [0.4, 0.6]) + assert not m["ok"] and abs(m["max_abs_diff"] - 0.1) < 1e-15 + + def test_length_mismatch_fails_closed(self): + m = rc.compare_weight_vectors([0.5, 0.5], [1.0]) + assert not m["ok"] and m["max_abs_diff"] is None + + def test_nonfinite_fails_closed(self): + m = rc.compare_weight_vectors([float("nan"), 0.5], [0.5, 0.5]) + assert not m["ok"] + + def test_empty_fails_closed(self): + assert not rc.compare_weight_vectors([], [])["ok"] + + +class TestGeneratorRefusesHardGatedPayloads: + @staticmethod + def _payload(flags): + return {"cells": {"callaway/small": {"flags": flags}}} + + def test_hard_flag_aborts_generation(self): + for flag in ( + "headline_se_gate_fail:python_pure:0.300", + "parity_fail", + "detail_keys_mismatch:event_study:python_rust", + "pure_rust_att_gate_fail:1e-6", + ): + with pytest.raises(SystemExit, match="REFUSING"): + gen.assert_no_hard_flags(self._payload([flag])) + + def test_rep_att_spread_aborts_generation(self): + with pytest.raises(SystemExit, match="REFUSING"): + gen.assert_no_hard_flags( + {"cells": {"basic/small": {"flags": ["python_rust:rep_att_gate_fail:1e-9"]}}} + ) + + def test_soft_flags_do_not_abort(self): + gen.assert_no_hard_flags(self._payload(["python_rust:cv_flag:0.200"])) + + def test_flag_free_payload_passes(self): + gen.assert_no_hard_flags(self._payload([]))