From 5c895253cb33eae965b83b3c17739ab97589d7a1 Mon Sep 17 00:00:00 2001
From: fabian-s <fabian.scheipl@googlemail.com>
Date: Tue, 21 Jul 2026 11:09:49 +0200
Subject: [PATCH 1/5] exact-CL2 default within the CL2 path (A18) + hardening
 test

Implemented by a codex subagent (gpt-5.6-terra, high), independently
re-verified: exact Bell-McCaffrey block B_g = I - 2 H_gg + (Ht^2)_gg with
(1-cap)^2 eigenvalue floor; cl2_adjustment = auto/exact/shortcut on pffr()
and coef.pffr() (auto: exact for G <= 100 and a dense-cost bound, visible
fallback to shortcut); diagnostics n_adjusted/min_block_eig/max_block_kappa
in fit metadata; docs + NEWS.

A18 hardening gate: 16 cells x {gaussian, poisson, binomial, Gamma} x
k {4,7} x unbalanced/influential-cluster designs - all SEs finite/positive,
3 influential Poisson/Gamma blocks hit the floor and match the uncapped
analogue to 8.8e-6, median SE ratio to shortcut 1.00-1.02
(hardening-exactcl2-results.csv). test-pffr-exactcl2.R: 19 PASS (re-run
independently); full test-pffr.R + sandwich/satterthwaite files passed.

Open: exact Bell-McCaffrey Satterthwaite df remains separate (documented).

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01WFfnErmWi9NCGCCJ8YwE95
---
 EXACTCL2-IMPLEMENTATION-NOTES.md    |  84 ++++++++++++
 NEWS.md                             |   7 +
 R/pffr-core.R                       | 196 +++++++++++++++++++++++-----
 R/pffr-methods.R                    |  15 ++-
 R/pffr.R                            |  16 ++-
 hardening-exactcl2-results.csv      |  17 +++
 hardening-exactcl2.R                | 141 ++++++++++++++++++++
 man/apply_sandwich_correction.Rd    |   6 +-
 man/coef.pffr.Rd                    |   7 +
 man/gam_sandwich_cluster_cl2.Rd     |  12 +-
 man/pffr.Rd                         |  13 +-
 man/pffr_build_metadata.Rd          |   3 +-
 man/pffr_compute_sandwich.Rd        |   6 +-
 man/pffr_vcov.Rd                    |   6 +-
 man/resolve_cl2_adjustment.Rd       |  35 +++++
 tests/testthat/test-pffr-exactcl2.R |  60 +++++++++
 tests/testthat/test-pffr.R          |  14 +-
 17 files changed, 591 insertions(+), 47 deletions(-)
 create mode 100644 EXACTCL2-IMPLEMENTATION-NOTES.md
 create mode 100644 hardening-exactcl2-results.csv
 create mode 100644 hardening-exactcl2.R
 create mode 100644 man/resolve_cl2_adjustment.Rd
 create mode 100644 tests/testthat/test-pffr-exactcl2.R

diff --git a/EXACTCL2-IMPLEMENTATION-NOTES.md b/EXACTCL2-IMPLEMENTATION-NOTES.md
new file mode 100644
index 00000000..940bdca4
--- /dev/null
+++ b/EXACTCL2-IMPLEMENTATION-NOTES.md
@@ -0,0 +1,84 @@
+# Exact-CL2 implementation notes
+
+Date: 2026-07-21
+
+## Change made
+
+`gam_sandwich_cluster_cl2()` now implements the exact Bell--McCaffrey
+leverage block
+
+```
+B_g = I - 2 H_gg + (H_t^2)_gg
+```
+
+where `(H_t^2)_gg` is evaluated as
+`(Xw_g Vp) (Xw' Xw) (Xw_g Vp)'`. Its eigendecomposition floors eigenvalues at
+`(1 - leverage_cap)^2`; the historical shortcut continues to cap the
+eigenvalues of `H_gg` and uses `(I - H_gg)^(-1/2)`.
+
+`pffr()` and `coef.pffr()` have a new `cl2_adjustment` argument:
+`"auto"` (default), `"exact"`, or `"shortcut"`. Fit metadata records the
+resolved adjustment plus `n_adjusted`, `min_block_eig`, and
+`max_block_kappa` for exact CL2. Legacy shortcut diagnostics remain available.
+
+## Default rule
+
+Within the existing CL2 path only, `"auto"` selects exact CL2 when:
+
+* `G <= 100`; and
+* `G * max(D_g) * p^2 <= 5e8`.
+
+The first bound confines exact CL2 to the small-to-moderate-cluster regime
+where its finite-sample correction is relevant. The second is a conservative
+dense-matrix proxy for the per-block multiplication
+`(Xw_g Vp) (Xw'Xw)`. Outside either bound the package visibly records and uses
+the shortcut. The package-level `sandwich = "auto"` policy is unchanged.
+
+## Tests run
+
+All R invocations used `OPENBLAS_NUM_THREADS=1 OMP_NUM_THREADS=1` and
+`nice -n 15`.
+
+* `Rscript hardening-exactcl2.R` — passed; wrote
+  `hardening-exactcl2-results.csv`.
+* Targeted test files run with `NOT_CRAN=true` — passed:
+  `test-pffr-exactcl2.R`, `test-pffr-s2-autopolicy.R`,
+  `test-pffr-sandwich-storage.R`, `test-pffr-sandwich-refit.R`,
+  `test-pffr-satterthwaite.R`, and `test-pffr.R`.
+* `devtools::load_all()` — passed.
+* `devtools::document()` — completed and regenerated the affected Rd files.
+  It printed pre-existing roxygen warnings about unrelated tags/method exports
+  in `pffr-core.R`, `poridge.R`, `pffr-methods.R`, and plot methods.
+
+## Hardening output summary
+
+The 16-cell table covers families Gaussian, Poisson, binomial, and Gamma;
+small/moderate response bases `k = 4, 7`; and unbalanced and
+influential-cluster designs. Each has five unbalanced clusters with curve
+counts `1, 1, 2, 3, 5`; the influential design gives one curve an extreme
+covariate value while retaining an ordinary-range outcome.
+
+Actual results:
+
+* All 16 exact covariances and coefficient SE vectors were finite and strictly
+  positive.
+* The median exact/shortcut SE ratio ranged from 1.00 to 1.02.
+* Three influential Poisson/Gamma blocks hit the shipped exact floor
+  (`n_adjusted = 1`); all other cells reported zero adjusted blocks.
+* The smallest observed pre-floor block eigenvalue was `3.20e-07`.
+* Exact versus the `cl2_exact_nocap` analogue (floor `tol = 1e-8`) differed by
+  at most `8.76e-06`; the nonzero differences occur in the floored cells.
+* The largest reported block condition proxy was `2.98e+06`.
+
+## Open questions
+
+The pointwise Satterthwaite degrees-of-freedom routine remains the established
+working-iid approximation based on the shortcut adjustment; this change
+updates the CL2 covariance block only. Exact Bell--McCaffrey df remains a
+separate methodological task, as already documented in the package.
+
+The requested local commits could not be created in this execution environment:
+Git returned `Unable to create .../.git/worktrees/refund-wt-s1/index.lock:
+Read-only file system`. The branch `exact-cl2-default` exists and all changes
+remain unstaged in this worktree; no unrelated `tests/testthat/Rplots.pdf`
+change was touched.
diff --git a/NEWS.md b/NEWS.md
index 72b2fbab..84073aa3 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -50,6 +50,13 @@
   report CL2 leverage diagnostics: the returned covariance carries
   `max_leverage` and `n_capped_clusters` attributes, and a warning is emitted
   when one or more clusters hit the leverage cap.
+* CL2 now uses the exact Bell--McCaffrey leverage block by default when the
+  finite-sample correction is relevant and the dense-block calculation is
+  affordable (at most 100 clusters and cost proxy `G * max(D_g) * p^2 <=
+  5e8`). At larger/prohibitive problems it uses the historical shortcut.
+  `cl2_adjustment = "exact"` or `"shortcut"` on `pffr()` and `coef.pffr()`
+  explicitly selects either variant. Exact-CL2 diagnostics include the number
+  of eigenvalue-floored blocks (`n_adjusted`).
 * AR(1) support improvements: `pffr()` now automatically switches to
   `algorithm = "bam"` and `method = "fREML"` when `rho` is supplied, and
   sets `discrete = TRUE` for non-Gaussian families.
diff --git a/R/pffr-core.R b/R/pffr-core.R
index 61c8567a..e225eba0 100644
--- a/R/pffr-core.R
+++ b/R/pffr-core.R
@@ -675,7 +675,8 @@ pffr_build_metadata <- function(
   ydata,
   sandwich,
   dof_correction = "none",
-  edf_type = "trace"
+  edf_type = "trace",
+  cl2_adjustment = "auto"
 ) {
   list(
     call = call,
@@ -698,6 +699,7 @@ pffr_build_metadata <- function(
     sandwich = sandwich,
     dof_correction = dof_correction,
     edf_type = edf_type,
+    cl2_adjustment = cl2_adjustment,
     # Covariance storage contract version: format 2 keeps $Vp/$Vc/$Ve
     # model-based ALWAYS; the robust covariance lives in $pffr$Vsandwich.
     cov_format = PFFR_COV_STORAGE_FORMAT,
@@ -1616,6 +1618,10 @@ gam_sandwich_cluster <- function(
 #'   If `FALSE` (default), use Bayesian sandwich (`B2 = Vp - Ve`).
 #' @param tol Eigenvalue floor for numerical stability.
 #' @param leverage_cap Cap for cluster leverage eigenvalues (< 1).
+#' @param cl2_adjustment CL2 leverage adjustment: `"auto"` (default) selects
+#'   the exact Bell--McCaffrey block where it is relevant and affordable,
+#'   `"exact"` forces that block, and `"shortcut"` uses the historical
+#'   per-cluster shortcut `(I - H_gg)^(-1/2)`.
 #' @param b2 Internal ablation switch (default `TRUE` = current behavior). When
 #'   `FALSE`, drop the additive Bayesian smoothing-bias term \eqn{B_2 = V_p -
 #'   V_e} (X5).
@@ -1623,8 +1629,10 @@ gam_sandwich_cluster <- function(
 #'   behavior). When `TRUE`, center the leverage-adjusted per-cluster
 #'   contributions \eqn{U_g^c = U_g - (\sum_g U_g)/G} before forming the meat
 #'   (X6).
-#' @returns A p x p covariance matrix with attributes `n_capped_clusters` and
-#'   `max_leverage` for the CL2 leverage diagnostic.
+#' @returns A p x p covariance matrix with leverage diagnostics. Exact CL2
+#'   returns `n_adjusted` (blocks floored at `(1 - leverage_cap)^2`),
+#'   `min_block_eig`, and `max_block_kappa`; the shortcut returns the legacy
+#'   `n_capped_clusters` and `max_leverage` attributes.
 #' @keywords internal
 gam_sandwich_cluster_cl2 <- function(
   b,
@@ -1632,6 +1640,7 @@ gam_sandwich_cluster_cl2 <- function(
   freq = FALSE,
   tol = 1e-8,
   leverage_cap = 0.999,
+  cl2_adjustment = c("auto", "exact", "shortcut"),
   b2 = TRUE,
   center_scores = FALSE
 ) {
@@ -1673,41 +1682,78 @@ gam_sandwich_cluster_cl2 <- function(
   G <- n_clusters_checked(cluster_id_work)
   groups <- unique(cluster_id_work)
 
-  Vp <- b$Vp
-  B2 <- if (freq) 0 else b$Vp - b$Ve
+  Vp_raw <- b$Vp
+  Vp <- 0.5 * (Vp_raw + t(Vp_raw))
   p <- ncol(Xw)
+  cl2_adjustment <- resolve_cl2_adjustment(
+    cl2_adjustment,
+    G = G,
+    maxDg = max(table(cluster_id_work)),
+    p = p
+  )
+  B2 <- if (freq) {
+    0
+  } else if (cl2_adjustment == "exact") {
+    Vp - b$Ve
+  } else {
+    Vp_raw - b$Ve
+  }
+  Cmat <- if (cl2_adjustment == "exact") crossprod(Xw) else NULL
   meat <- matrix(0, nrow = p, ncol = p)
   Usum <- numeric(p)
   n_capped_clusters <- 0L
   max_leverage <- NA_real_
+  n_adjusted <- 0L
+  min_block_eig <- Inf
+  max_block_kappa <- 0
 
   for (g in groups) {
     idx <- which(cluster_id_work == g)
     Xwg <- Xw[idx, , drop = FALSE]
     zg <- z[idx]
 
-    Hgg <- Xwg %*% Vp %*% t(Xwg)
+    Tg <- Xwg %*% Vp
+    Hgg <- Tg %*% t(Xwg)
     Hgg <- 0.5 * (Hgg + t(Hgg))
 
-    ee_H <- eigen(Hgg, symmetric = TRUE)
-    cluster_max_leverage <- max(ee_H$values, na.rm = TRUE)
-    if (is.finite(cluster_max_leverage)) {
-      max_leverage <- if (is.na(max_leverage)) {
-        cluster_max_leverage
-      } else {
-        max(max_leverage, cluster_max_leverage)
+    if (cl2_adjustment == "shortcut") {
+      ee_H <- eigen(Hgg, symmetric = TRUE)
+      cluster_max_leverage <- max(ee_H$values, na.rm = TRUE)
+      if (is.finite(cluster_max_leverage)) {
+        max_leverage <- if (is.na(max_leverage)) {
+          cluster_max_leverage
+        } else {
+          max(max_leverage, cluster_max_leverage)
+        }
       }
+      if (any(ee_H$values > leverage_cap, na.rm = TRUE)) {
+        n_capped_clusters <- n_capped_clusters + 1L
+        ee_H$values <- pmin(ee_H$values, leverage_cap)
+        Hgg <- ee_H$vectors %*%
+          diag(ee_H$values, nrow = length(ee_H$values)) %*%
+          t(ee_H$vectors)
+      }
+      Ag <- sym_inv_sqrt(diag(length(idx)) - Hgg, tol = tol)
+    } else {
+      # The exact Bell--McCaffrey block is B_g = I - 2 H_gg +
+      # (H_t^2)_gg, with (H_t^2)_gg = (Xw_g Vp) (Xw' Xw)
+      # (Xw_g Vp)'. This is a squared operator, so the gain-matched analogue
+      # of capping H_gg at `leverage_cap` floors B_g at (1 - cap)^2.
+      Bg <- diag(length(idx)) - 2 * Hgg + Tg %*% Cmat %*% t(Tg)
+      ee_B <- eigen(0.5 * (Bg + t(Bg)), symmetric = TRUE)
+      block_min <- min(ee_B$values)
+      min_block_eig <- min(min_block_eig, block_min)
+      max_block_kappa <- max(
+        max_block_kappa,
+        abs(max(ee_B$values)) / max(abs(block_min), 1e-300)
+      )
+      floor_val <- (1 - leverage_cap)^2
+      if (block_min < floor_val) n_adjusted <- n_adjusted + 1L
+      ee_B$values <- pmax(ee_B$values, floor_val)
+      Ag <- ee_B$vectors %*%
+        diag(1 / sqrt(ee_B$values), nrow = length(ee_B$values)) %*%
+        t(ee_B$vectors)
     }
-    if (any(ee_H$values > leverage_cap, na.rm = TRUE)) {
-      n_capped_clusters <- n_capped_clusters + 1L
-      ee_H$values <- pmin(ee_H$values, leverage_cap)
-      Hgg <- ee_H$vectors %*%
-        diag(ee_H$values, nrow = length(ee_H$values)) %*%
-        t(ee_H$vectors)
-    }
-
-    Mg <- diag(length(idx)) - Hgg
-    Ag <- sym_inv_sqrt(Mg, tol = tol)
     Ug <- crossprod(Xwg, Ag %*% zg)
     meat <- meat + Ug %*% t(Ug)
     Usum <- Usum + as.vector(Ug)
@@ -1723,9 +1769,14 @@ gam_sandwich_cluster_cl2 <- function(
   V <- hc1 * Vp %*% meat %*% Vp
   if (isTRUE(b2)) V <- V + B2
   V <- 0.5 * (V + t(V))
+  attr(V, "cl2_adjustment") <- cl2_adjustment
   attr(V, "n_capped_clusters") <- n_capped_clusters
   attr(V, "max_leverage") <- max_leverage
-  if (n_capped_clusters > 0) {
+  attr(V, "n_adjusted") <- n_adjusted
+  attr(V, "min_block_eig") <- if (is.finite(min_block_eig)) min_block_eig else
+    NA_real_
+  attr(V, "max_block_kappa") <- max_block_kappa
+  if (cl2_adjustment == "shortcut" && n_capped_clusters > 0) {
     max_leverage_label <- if (is.finite(max_leverage)) {
       sprintf("%.3f", max_leverage)
     } else {
@@ -1750,6 +1801,41 @@ gam_sandwich_cluster_cl2 <- function(
   V
 }
 
+#' Resolve the CL2 leverage adjustment
+#'
+#' Exact CL2 is the default when the finite-sample correction can matter
+#' (`G <= 100`) and the dense block multiplication is affordable. Its
+#' conservative cost proxy `G * maxDg * p^2` covers the dominant
+#' `(Xw_g Vp) (Xw'Xw)` multiplication; 5e8 scalar operations keeps this below
+#' the cost of a typical pffr fit on the supported auto-CL2 grids. Large-G
+#' CL2 uses the historical shortcut because the correction tends to zero and
+#' the exact block has no material finite-sample benefit there.
+#'
+#' @param cl2_adjustment One of `"auto"`, `"exact"`, or `"shortcut"`.
+#' @param G Number of clusters.
+#' @param maxDg Largest cluster size.
+#' @param p Coefficient-space dimension.
+#' @returns `"exact"` or `"shortcut"`.
+#' @keywords internal
+resolve_cl2_adjustment <- function(
+  cl2_adjustment = c("auto", "exact", "shortcut"),
+  G,
+  maxDg,
+  p
+) {
+  cl2_adjustment <- match.arg(cl2_adjustment)
+  if (cl2_adjustment != "auto") return(cl2_adjustment)
+
+  cost <- G * maxDg * p^2
+  if (
+    is.finite(G) && is.finite(maxDg) && is.finite(p) && G <= 100 && cost <= 5e8
+  ) {
+    "exact"
+  } else {
+    "shortcut"
+  }
+}
+
 #' Resolve a pffr fit's covariance matrices to a canonical representation
 #'
 #' The current storage contract (format 2) keeps `Vp`/`Vc`/`Ve` model-based
@@ -1899,6 +1985,8 @@ restore_model_cov <- function(object) {
 #' @param center_scores Internal ablation switch (default `FALSE`); center the
 #'   per-cluster score sums before the meat when `TRUE` (X6). Applies to
 #'   `"cluster"`/`"cl2"` only.
+#' @param cl2_adjustment CL2 leverage adjustment (`"auto"`, `"exact"`, or
+#'   `"shortcut"`), used only for `type = "cl2"`.
 #' @returns A covariance matrix (with CL2 leverage attributes for `type =
 #'   "cl2"`).
 #' @keywords internal
@@ -1910,7 +1998,8 @@ pffr_compute_sandwich <- function(
   dof_correction = "none",
   edf_type = "trace",
   b2 = TRUE,
-  center_scores = FALSE
+  center_scores = FALSE,
+  cl2_adjustment = "auto"
 ) {
   switch(
     type,
@@ -1928,7 +2017,8 @@ pffr_compute_sandwich <- function(
       cluster_id,
       freq = freq,
       b2 = b2,
-      center_scores = center_scores
+      center_scores = center_scores,
+      cl2_adjustment = cl2_adjustment
     ),
     hc = mgcv::vcov.gam(b, sandwich = TRUE, freq = freq),
     none = if (freq) b$Ve else (b$Vc %||% b$Vp),
@@ -1960,6 +2050,8 @@ pffr_compute_sandwich <- function(
 #' @param center_scores Internal ablation switch (default `FALSE` = current
 #'   behavior). When `TRUE`, center the per-cluster score sums before forming
 #'   the cluster/CL2 meat (X6). Same cache semantics as `b2`.
+#' @param cl2_adjustment CL2 leverage adjustment. `NULL` (default) inherits
+#'   the fitted choice; otherwise one of `"auto"`, `"exact"`, or `"shortcut"`.
 #' @returns A covariance matrix.
 #' @keywords internal
 pffr_vcov <- function(
@@ -1970,7 +2062,8 @@ pffr_vcov <- function(
   dof_correction = NULL,
   edf_type = NULL,
   b2 = TRUE,
-  center_scores = FALSE
+  center_scores = FALSE,
+  cl2_adjustment = NULL
 ) {
   # Ablation variants (X5/X6) bypass the fit-time and recompute caches entirely,
   # so they never overwrite or shadow the standard cached matrices.
@@ -1990,6 +2083,13 @@ pffr_vcov <- function(
 
   dof_correction <- dof_correction %||% (object$pffr$dof_correction %||% "none")
   edf_type <- edf_type %||% (object$pffr$edf_type %||% "trace")
+  cl2_adjustment_explicit <- !is.null(cl2_adjustment)
+  cl2_adjustment <- cl2_adjustment %||%
+    (object$pffr$cl2_adjustment %||% "auto")
+  cl2_adjustment <- match.arg(
+    cl2_adjustment,
+    c("auto", "exact", "shortcut")
+  )
 
   # Serve the cached fit-time robust matrix when the request matches it exactly.
   opts_match <- if (requested == "cluster") {
@@ -1999,11 +2099,20 @@ pffr_vcov <- function(
   } else {
     TRUE
   }
+  adjustment_match <- requested != "cl2" ||
+    !cl2_adjustment_explicit ||
+    identical(
+      cl2_adjustment,
+      object$pffr$sandwich_info$cl2_adjustment %||%
+        object$pffr$cl2_adjustment %||%
+        "shortcut"
+    )
   if (
     !ablation &&
       is.null(cluster) &&
       identical(requested, canon$fit_type) &&
       opts_match &&
+      adjustment_match &&
       !is.null(canon$Vsandwich)
   ) {
     return(
@@ -2024,6 +2133,17 @@ pffr_vcov <- function(
   key <- if (!ablation && is.null(cluster)) {
     if (requested == "cluster" && (freq || dof_correction != "none")) {
       paste(requested, freq, dof_correction, edf_type, sep = "|")
+    } else if (requested == "cl2") {
+      if (cl2_adjustment == "auto") {
+        if (freq) paste(requested, "freq", sep = "|") else requested
+      } else {
+        paste(
+          requested,
+          cl2_adjustment,
+          if (freq) "freq" else "bayes",
+          sep = "|"
+        )
+      }
     } else if (freq) {
       paste(requested, "freq", sep = "|")
     } else {
@@ -2054,7 +2174,8 @@ pffr_vcov <- function(
     dof_correction = dof_correction,
     edf_type = edf_type,
     b2 = b2,
-    center_scores = center_scores
+    center_scores = center_scores,
+    cl2_adjustment = cl2_adjustment
   )
   if (!is.null(cache) && !is.null(key)) {
     cache[[key]] <- V
@@ -2956,6 +3077,8 @@ pffr_upgrade_fit <- function(object) {
 #'   warning) for `"cl2"`, which already corrects per-cluster leverage.
 #' @param edf_type Which EDF the `"edf"` correction uses (`"trace"`/`"edf2"`/
 #'   `"basis"`).
+#' @param cl2_adjustment CL2 leverage adjustment (`"auto"` (default),
+#'   `"exact"`, or `"shortcut"`), used only for `type = "cl2"`.
 #' @returns Model with the robust covariance stored in `$pffr$Vsandwich`.
 #' @keywords internal
 apply_sandwich_correction <- function(
@@ -2963,7 +3086,8 @@ apply_sandwich_correction <- function(
   algorithm,
   type = "cluster",
   dof_correction = "none",
-  edf_type = "trace"
+  edf_type = "trace",
+  cl2_adjustment = "auto"
 ) {
   gam_obj <- if (as.character(algorithm) %in% c("gamm4", "gamm")) m$gam else m
 
@@ -2992,7 +3116,8 @@ apply_sandwich_correction <- function(
     cluster_id,
     freq = FALSE,
     dof_correction = dof_correction,
-    edf_type = edf_type
+    edf_type = edf_type,
+    cl2_adjustment = cl2_adjustment
   )
   Vsw_freq <- pffr_compute_sandwich(
     bread,
@@ -3000,11 +3125,14 @@ apply_sandwich_correction <- function(
     cluster_id,
     freq = TRUE,
     dof_correction = dof_correction,
-    edf_type = edf_type
+    edf_type = edf_type,
+    cl2_adjustment = cl2_adjustment
   )
 
   n_capped <- attr(Vsw, "n_capped_clusters") %||% 0L
   max_lev <- attr(Vsw, "max_leverage") %||% NA_real_
+  n_adjusted <- attr(Vsw, "n_adjusted") %||% 0L
+  resolved_adjustment <- attr(Vsw, "cl2_adjustment") %||% NA_character_
 
   gam_obj$pffr$Vsandwich <- Vsw
   gam_obj$pffr$Vsandwich_freq <- Vsw_freq
@@ -3014,6 +3142,10 @@ apply_sandwich_correction <- function(
     G = if (!is.null(cluster_id)) length(unique(cluster_id)) else NA_integer_,
     n_capped = n_capped,
     max_leverage = max_lev,
+    cl2_adjustment = resolved_adjustment,
+    n_adjusted = n_adjusted,
+    min_block_eig = attr(Vsw, "min_block_eig") %||% NA_real_,
+    max_block_kappa = attr(Vsw, "max_block_kappa") %||% NA_real_,
     dof_correction = if (type == "cluster") dof_correction else "none",
     edf_type = edf_type,
     version = as.character(utils::packageVersion("refund")),
@@ -3021,6 +3153,8 @@ apply_sandwich_correction <- function(
   )
   # Keep the legacy CL2 leverage-cap diagnostic slot populated.
   gam_obj$pffr$cl2_n_capped <- if (type == "cl2") n_capped else NULL
+  gam_obj$pffr$cl2_adjustment <- if (type == "cl2") resolved_adjustment else
+    NULL
   # Fresh cache for on-demand recomputation of other sandwich types
   # (fit$pffr$Vsandwich_cache[[type]]).
   gam_obj$pffr$Vsandwich_cache <- new.env(parent = emptyenv())
diff --git a/R/pffr-methods.R b/R/pffr-methods.R
index 05736df3..3324aecc 100755
--- a/R/pffr-methods.R
+++ b/R/pffr-methods.R
@@ -1291,6 +1291,11 @@ compute_pointwise_ci <- function(
 #' @param edf_type Which EDF the \code{"edf"} correction uses
 #'   (\code{"trace"}/\code{"edf2"}/\code{"basis"}; see \code{\link{pffr}}).
 #'   Defaults to \code{NULL} (inherit from the fit).
+#' @param cl2_adjustment Leverage adjustment within \code{sandwich = "cl2"}.
+#'   \code{NULL} (default) inherits the fit-time choice; \code{"auto"}
+#'   applies the documented feasibility rule, while \code{"exact"} and
+#'   \code{"shortcut"} force either CL2 variant. Supplying a value forces a
+#'   covariance recomputation.
 #' @param seWithMean logical, defaults to TRUE. Include uncertainty about the intercept/overall mean in  standard errors returned for smooth components?
 #' @param n1 see below
 #' @param n2 see below
@@ -1368,6 +1373,7 @@ coef.pffr <- function(
   cluster = NULL,
   dof_correction = NULL,
   edf_type = NULL,
+  cl2_adjustment = NULL,
   seWithMean = TRUE,
   n1 = 100,
   n2 = 40,
@@ -1398,6 +1404,12 @@ coef.pffr <- function(
   if (is.null(edf_type)) edf_type <- model_edf_type
   dof_correction <- match.arg(dof_correction, c("none", "edf"))
   edf_type <- match.arg(edf_type, c("trace", "edf2", "basis"))
+  if (!is.null(cl2_adjustment)) {
+    cl2_adjustment <- match.arg(
+      cl2_adjustment,
+      c("auto", "exact", "shortcut")
+    )
+  }
   # Only warn when the user *explicitly* asked for a dof correction on a
   # non-cluster sandwich; an inherited "edf" (from the fit) stays silent.
   if (dof_explicitly_set && dof_correction != "none" && sandwich != "cluster") {
@@ -1586,7 +1598,8 @@ coef.pffr <- function(
       dof_correction = dof_correction,
       edf_type = edf_type,
       b2 = b2,
-      center_scores = center_scores
+      center_scores = center_scores,
+      cl2_adjustment = cl2_adjustment
     )
 
     # Pointwise critical-value reference (S3). `crit` selects the pointwise
diff --git a/R/pffr.R b/R/pffr.R
index 0f6668da..a85fb987 100644
--- a/R/pffr.R
+++ b/R/pffr.R
@@ -180,7 +180,9 @@
 #'   handles both heteroskedasticity and within-curve autocorrelation — the
 #'   recommended choice for functional data when a fixed estimator is wanted.
 #'   \code{"cl2"}: leverage-adjusted cluster-robust sandwich (Bell-McCaffrey
-#'   style CL2), mainly relevant in smaller samples.
+#'   style CL2), mainly relevant in smaller samples. Within CL2, the exact
+#'   block calculation is used by default where feasible; see
+#'   \code{cl2_adjustment} below.
 #'   \code{"hc"}: observation-level HC sandwich via
 #'   \code{\link[mgcv]{vcov.gam}(sandwich = TRUE)}, which corrects for
 #'   heteroskedasticity but ignores within-curve correlation.
@@ -220,6 +222,13 @@
 #'   the penalized hat), \code{"edf2"} (mgcv's bias-corrected EDF), or
 #'   \code{"basis"} (the basis dimension). Only relevant when
 #'   \code{dof_correction = "edf"}.
+#' @param cl2_adjustment Leverage adjustment within \code{sandwich = "cl2"}:
+#'   \code{"auto"} (default) uses exact CL2 when \eqn{G \le 100} and the
+#'   dense-block cost proxy \eqn{G\max_g D_g p^2} is at most \eqn{5\times
+#'   10^8}; it otherwise falls back to the historical shortcut. The exact block
+#'   is \eqn{B_g = I - 2H_{gg} + (H_t^2)_{gg}}, with eigenvalues floored at
+#'   \eqn{(1 - 0.999)^2}. \code{"exact"} and \code{"shortcut"} force either
+#'   variant. Ignored unless the resolved sandwich is \code{"cl2"}.
 #' @param ... additional arguments that are valid for \code{\link[mgcv]{gam}},
 #'   \code{\link[mgcv]{bam}}, \code{'\link[gamm4]{gamm4}'} or
 #'   \code{'\link[mgcv]{jagam}'}. \code{subset} is not implemented.
@@ -320,6 +329,7 @@ pffr <- function(
   sandwich = c("auto", "cluster", "cl2", "hc", "none"),
   dof_correction = c("none", "edf"),
   edf_type = c("trace", "edf2", "basis"),
+  cl2_adjustment = c("auto", "exact", "shortcut"),
   ...
 ) {
   call <- match.call()
@@ -330,6 +340,7 @@ pffr <- function(
   sandwich <- match.arg(sandwich)
   dof_correction <- match.arg(dof_correction)
   edf_type <- match.arg(edf_type)
+  cl2_adjustment <- match.arg(cl2_adjustment)
   # "auto" is resolved after fitting (needs G / max D_g); defer the CR1 dof
   # compatibility check to then, so a legitimate auto -> "cluster" keeps its dof.
   if (dof_correction != "none" && !(sandwich %in% c("cluster", "auto"))) {
@@ -503,6 +514,7 @@ pffr <- function(
     prep$algorithm,
     type = sandwich,
     dof_correction = dof_correction,
-    edf_type = edf_type
+    edf_type = edf_type,
+    cl2_adjustment = cl2_adjustment
   )
 }
diff --git a/hardening-exactcl2-results.csv b/hardening-exactcl2-results.csv
new file mode 100644
index 00000000..39f55aab
--- /dev/null
+++ b/hardening-exactcl2-results.csv
@@ -0,0 +1,17 @@
+"family","k","design","G","max_cluster_size","n_adjusted","min_block_eig","max_block_kappa","max_abs_exact_nocap","median_se_ratio_to_shortcut","se_finite_positive","se_ratio_sane"
+"gaussian",4,"unbalanced",5,50,0,0.453590234808712,2.20463299969789,0,1.01263947332551,TRUE,TRUE
+"poisson",4,"unbalanced",5,50,0,0.348940004869249,2.82730371704937,0,1.0079082080914,TRUE,TRUE
+"binomial",4,"unbalanced",5,50,0,0.552153568687061,1.81109035005941,0,1.0010220332576,TRUE,TRUE
+"Gamma",4,"unbalanced",5,50,0,0.332486151450086,3.00764406468859,0,1.00903670061469,TRUE,TRUE
+"gaussian",7,"unbalanced",5,50,0,0.374043657671271,2.67348471091803,0,1.0189182605098,TRUE,TRUE
+"poisson",7,"unbalanced",5,50,0,0.301435120875162,3.31746346310504,0,1.00000000257308,TRUE,TRUE
+"binomial",7,"unbalanced",5,50,0,0.478396994803052,2.09031413420915,0,1.00000000004754,TRUE,TRUE
+"Gamma",7,"unbalanced",5,50,0,0.386030967291895,2.590465752049,0,1.00648632821506,TRUE,TRUE
+"gaussian",4,"influential",5,50,0,1.39826045642e-06,709818.235108635,0,1.00791990080075,TRUE,TRUE
+"poisson",4,"influential",5,50,1,3.20272329212834e-07,2982577.58727008,7.06659823389664e-06,1.01755206252338,TRUE,TRUE
+"binomial",4,"influential",5,50,0,1.58358505475409e-06,623867.247044085,0,1.01039988193779,TRUE,TRUE
+"Gamma",4,"influential",5,50,0,2.35791791958206e-06,423665.797661914,0,1.01294453124424,TRUE,TRUE
+"gaussian",7,"influential",5,50,0,1.12859513423658e-06,882176.398770039,0,1.01074965381931,TRUE,TRUE
+"poisson",7,"influential",5,50,1,7.21418514260108e-07,1385176.29210241,8.75780657148509e-07,1.01500324062661,TRUE,TRUE
+"binomial",7,"influential",5,50,0,1.35355708907881e-06,738794.066446026,0,1.00000000006042,TRUE,TRUE
+"Gamma",7,"influential",5,50,1,7.39633434387876e-07,1348528.03666834,4.54114080495374e-07,1.00490798665056,TRUE,TRUE
diff --git a/hardening-exactcl2.R b/hardening-exactcl2.R
new file mode 100644
index 00000000..96a10d1e
--- /dev/null
+++ b/hardening-exactcl2.R
@@ -0,0 +1,141 @@
+# Scoped exact-CL2 hardening gate (DESIGN amendment A18).
+#
+# Run from the package root with:
+# OPENBLAS_NUM_THREADS=1 OMP_NUM_THREADS=1 nice -n 15 Rscript hardening-exactcl2.R
+# The script writes hardening-exactcl2-results.csv in the package root.
+
+make_exactcl2_fixture <- function(family_name, k, design) {
+  set.seed(
+    match(family_name, c("gaussian", "poisson", "binomial", "Gamma")) *
+      100 +
+      k +
+      match(design, c("unbalanced", "influential"))
+  )
+  n_curve <- 12L
+  n_grid <- 10L
+  grid <- seq(0, 1, length.out = n_grid)
+  xlin <- stats::rnorm(n_curve)
+  eta <- outer(xlin, rep(0.35, n_grid)) +
+    matrix(rep(0.2 * sin(2 * pi * grid), each = n_curve), n_curve, n_grid) -
+    0.4
+  Y <- switch(
+    family_name,
+    gaussian = eta + matrix(stats::rnorm(n_curve * n_grid, sd = 0.35), n_curve),
+    poisson = matrix(
+      stats::rpois(n_curve * n_grid, lambda = exp(pmin(eta, 2))),
+      n_curve
+    ),
+    binomial = matrix(
+      stats::rbinom(n_curve * n_grid, 1, stats::plogis(eta)),
+      n_curve
+    ),
+    Gamma = matrix(
+      stats::rgamma(n_curve * n_grid, shape = 5, scale = exp(pmin(eta, 2)) / 5),
+      n_curve
+    )
+  )
+  # Keep the response in its ordinary range, then make one covariate curve
+  # extreme. This isolates high leverage from a separable/outcome-saturation
+  # failure and exercises the shipped eigenvalue floor (especially Gaussian,
+  # Poisson, and Gamma at k = 9).
+  if (design == "influential") xlin[1] <- 3000
+  # Five deliberately unbalanced independent clusters, with cluster 1 made
+  # influential above. Their curve counts 1, 1, 2, 3, 5 yield unequal blocks.
+  cluster <- rep(seq_len(5), times = c(1, 1, 2, 3, 5))
+  list(
+    data = data.frame(Y = I(Y), xlin = xlin),
+    yind = grid,
+    cluster = cluster,
+    family = switch(
+      family_name,
+      gaussian = stats::gaussian(),
+      poisson = stats::poisson(),
+      binomial = stats::binomial(),
+      Gamma = stats::Gamma(link = "log")
+    ),
+    k = k
+  )
+}
+
+run_exactcl2_hardening <- function() {
+  cases <- expand.grid(
+    family = c("gaussian", "poisson", "binomial", "Gamma"),
+    k = c(4L, 7L),
+    design = c("unbalanced", "influential"),
+    stringsAsFactors = FALSE
+  )
+  rows <- vector("list", nrow(cases))
+
+  for (i in seq_len(nrow(cases))) {
+    case <- cases[i, ]
+    fixture <- make_exactcl2_fixture(case$family, case$k, case$design)
+    fit <- suppressWarnings(suppressMessages(refund::pffr(
+      Y ~ xlin,
+      data = fixture$data,
+      yind = fixture$yind,
+      family = fixture$family,
+      bs.yindex = list(bs = "ps", k = fixture$k, m = c(2, 1)),
+      sandwich = "none"
+    )))
+
+    # `1 - sqrt(tol)` makes the exact floor equal tol: Study EX's
+    # cl2_exact_nocap comparator without duplicating the package implementation.
+    V_exact <- suppressWarnings(refund:::pffr_vcov(
+      fit,
+      sandwich = "cl2",
+      cluster = fixture$cluster,
+      cl2_adjustment = "exact"
+    ))
+    V_nocap <- suppressWarnings(refund:::gam_sandwich_cluster_cl2(
+      refund:::pffr_model_based_gam(fit),
+      refund:::build_cluster_id(fit$pffr, cluster = fixture$cluster),
+      cl2_adjustment = "exact",
+      leverage_cap = 1 - sqrt(1e-8)
+    ))
+    V_shortcut <- suppressWarnings(refund:::pffr_vcov(
+      fit,
+      sandwich = "cl2",
+      cluster = fixture$cluster,
+      cl2_adjustment = "shortcut"
+    ))
+
+    se_exact <- sqrt(pmax(diag(V_exact), 0))
+    se_shortcut <- sqrt(pmax(diag(V_shortcut), 0))
+    positive <- is.finite(se_exact) &
+      is.finite(se_shortcut) &
+      se_exact > 0 &
+      se_shortcut > 0
+    ratio <- se_exact[positive] / se_shortcut[positive]
+    rows[[i]] <- data.frame(
+      family = case$family,
+      k = case$k,
+      design = case$design,
+      G = length(unique(fixture$cluster)),
+      max_cluster_size = max(table(fixture$cluster)) * length(fixture$yind),
+      n_adjusted = attr(V_exact, "n_adjusted"),
+      min_block_eig = attr(V_exact, "min_block_eig"),
+      max_block_kappa = attr(V_exact, "max_block_kappa"),
+      max_abs_exact_nocap = max(abs(V_exact - V_nocap)),
+      median_se_ratio_to_shortcut = stats::median(ratio),
+      se_finite_positive = all(positive),
+      se_ratio_sane = length(ratio) > 0 &&
+        min(ratio) > 1e-3 &&
+        max(ratio) < 1e3,
+      stringsAsFactors = FALSE
+    )
+  }
+  do.call(rbind, rows)
+}
+
+if (sys.nframe() == 0L) {
+  devtools::load_all(quiet = TRUE)
+  results <- run_exactcl2_hardening()
+  utils::write.csv(results, "hardening-exactcl2-results.csv", row.names = FALSE)
+  print(results)
+  stopifnot(
+    all(results$se_finite_positive),
+    all(results$se_ratio_sane),
+    all(is.finite(results$n_adjusted)),
+    any(results$n_adjusted > 0)
+  )
+}
diff --git a/man/apply_sandwich_correction.Rd b/man/apply_sandwich_correction.Rd
index 2a497f12..904a494d 100644
--- a/man/apply_sandwich_correction.Rd
+++ b/man/apply_sandwich_correction.Rd
@@ -9,7 +9,8 @@ apply_sandwich_correction(
   algorithm,
   type = "cluster",
   dof_correction = "none",
-  edf_type = "trace"
+  edf_type = "trace",
+  cl2_adjustment = "auto"
 )
 }
 \arguments{
@@ -26,6 +27,9 @@ warning) for `"cl2"`, which already corrects per-cluster leverage.}
 
 \item{edf_type}{Which EDF the `"edf"` correction uses (`"trace"`/`"edf2"`/
 `"basis"`).}
+
+\item{cl2_adjustment}{CL2 leverage adjustment (`"auto"` (default),
+`"exact"`, or `"shortcut"`), used only for `type = "cl2"`.}
 }
 \value{
 Model with the robust covariance stored in `$pffr$Vsandwich`.
diff --git a/man/coef.pffr.Rd b/man/coef.pffr.Rd
index 11eceb3b..a6e8d4de 100644
--- a/man/coef.pffr.Rd
+++ b/man/coef.pffr.Rd
@@ -13,6 +13,7 @@
   cluster = NULL,
   dof_correction = NULL,
   edf_type = NULL,
+  cl2_adjustment = NULL,
   seWithMean = TRUE,
   n1 = 100,
   n2 = 40,
@@ -68,6 +69,12 @@ Ignored (with a warning) for \code{sandwich} other than \code{"cluster"}.}
 (\code{"trace"}/\code{"edf2"}/\code{"basis"}; see \code{\link{pffr}}).
 Defaults to \code{NULL} (inherit from the fit).}
 
+\item{cl2_adjustment}{Leverage adjustment within \code{sandwich = "cl2"}.
+\code{NULL} (default) inherits the fit-time choice; \code{"auto"}
+applies the documented feasibility rule, while \code{"exact"} and
+\code{"shortcut"} force either CL2 variant. Supplying a value forces a
+covariance recomputation.}
+
 \item{seWithMean}{logical, defaults to TRUE. Include uncertainty about the intercept/overall mean in  standard errors returned for smooth components?}
 
 \item{n1}{see below}
diff --git a/man/gam_sandwich_cluster_cl2.Rd b/man/gam_sandwich_cluster_cl2.Rd
index f5d1ac61..3bd0e45f 100644
--- a/man/gam_sandwich_cluster_cl2.Rd
+++ b/man/gam_sandwich_cluster_cl2.Rd
@@ -10,6 +10,7 @@ gam_sandwich_cluster_cl2(
   freq = FALSE,
   tol = 1e-08,
   leverage_cap = 0.999,
+  cl2_adjustment = c("auto", "exact", "shortcut"),
   b2 = TRUE,
   center_scores = FALSE
 )
@@ -27,6 +28,11 @@ If `FALSE` (default), use Bayesian sandwich (`B2 = Vp - Ve`).}
 
 \item{leverage_cap}{Cap for cluster leverage eigenvalues (< 1).}
 
+\item{cl2_adjustment}{CL2 leverage adjustment: `"auto"` (default) selects
+the exact Bell--McCaffrey block where it is relevant and affordable,
+`"exact"` forces that block, and `"shortcut"` uses the historical
+per-cluster shortcut `(I - H_gg)^(-1/2)`.}
+
 \item{b2}{Internal ablation switch (default `TRUE` = current behavior). When
 `FALSE`, drop the additive Bayesian smoothing-bias term \eqn{B_2 = V_p -
 V_e} (X5).}
@@ -37,8 +43,10 @@ contributions \eqn{U_g^c = U_g - (\sum_g U_g)/G} before forming the meat
 (X6).}
 }
 \value{
-A p x p covariance matrix with attributes `n_capped_clusters` and
-  `max_leverage` for the CL2 leverage diagnostic.
+A p x p covariance matrix with leverage diagnostics. Exact CL2
+  returns `n_adjusted` (blocks floored at `(1 - leverage_cap)^2`),
+  `min_block_eig`, and `max_block_kappa`; the shortcut returns the legacy
+  `n_capped_clusters` and `max_leverage` attributes.
 }
 \description{
 Computes a cluster-robust covariance matrix with a Bell-McCaffrey style
diff --git a/man/pffr.Rd b/man/pffr.Rd
index 4fe62cae..769aae24 100644
--- a/man/pffr.Rd
+++ b/man/pffr.Rd
@@ -17,6 +17,7 @@ pffr(
   sandwich = c("auto", "cluster", "cl2", "hc", "none"),
   dof_correction = c("none", "edf"),
   edf_type = c("trace", "edf2", "basis"),
+  cl2_adjustment = c("auto", "exact", "shortcut"),
   ...
 )
 }
@@ -81,7 +82,9 @@ penalty.}
   handles both heteroskedasticity and within-curve autocorrelation — the
   recommended choice for functional data when a fixed estimator is wanted.
   \code{"cl2"}: leverage-adjusted cluster-robust sandwich (Bell-McCaffrey
-  style CL2), mainly relevant in smaller samples.
+  style CL2), mainly relevant in smaller samples. Within CL2, the exact
+  block calculation is used by default where feasible; see
+  \code{cl2_adjustment} below.
   \code{"hc"}: observation-level HC sandwich via
   \code{\link[mgcv]{vcov.gam}(sandwich = TRUE)}, which corrects for
   heteroskedasticity but ignores within-curve correlation.
@@ -124,6 +127,14 @@ the penalized hat), \code{"edf2"} (mgcv's bias-corrected EDF), or
 \code{"basis"} (the basis dimension). Only relevant when
 \code{dof_correction = "edf"}.}
 
+\item{cl2_adjustment}{Leverage adjustment within \code{sandwich = "cl2"}:
+\code{"auto"} (default) uses exact CL2 when \eqn{G \le 100} and the
+dense-block cost proxy \eqn{G\max_g D_g p^2} is at most \eqn{5\times
+10^8}; it otherwise falls back to the historical shortcut. The exact block
+is \eqn{B_g = I - 2H_{gg} + (H_t^2)_{gg}}, with eigenvalues floored at
+\eqn{(1 - 0.999)^2}. \code{"exact"} and \code{"shortcut"} force either
+variant. Ignored unless the resolved sandwich is \code{"cl2"}.}
+
 \item{...}{additional arguments that are valid for \code{\link[mgcv]{gam}},
 \code{\link[mgcv]{bam}}, \code{'\link[gamm4]{gamm4}'} or
 \code{'\link[mgcv]{jagam}'}. \code{subset} is not implemented.}
diff --git a/man/pffr_build_metadata.Rd b/man/pffr_build_metadata.Rd
index 925ec667..9d183451 100644
--- a/man/pffr_build_metadata.Rd
+++ b/man/pffr_build_metadata.Rd
@@ -24,7 +24,8 @@ pffr_build_metadata(
   ydata,
   sandwich,
   dof_correction = "none",
-  edf_type = "trace"
+  edf_type = "trace",
+  cl2_adjustment = "auto"
 )
 }
 \description{
diff --git a/man/pffr_compute_sandwich.Rd b/man/pffr_compute_sandwich.Rd
index 5d60df7e..6069f230 100644
--- a/man/pffr_compute_sandwich.Rd
+++ b/man/pffr_compute_sandwich.Rd
@@ -12,7 +12,8 @@ pffr_compute_sandwich(
   dof_correction = "none",
   edf_type = "trace",
   b2 = TRUE,
-  center_scores = FALSE
+  center_scores = FALSE,
+  cl2_adjustment = "auto"
 )
 }
 \arguments{
@@ -34,6 +35,9 @@ only).}
 \item{center_scores}{Internal ablation switch (default `FALSE`); center the
 per-cluster score sums before the meat when `TRUE` (X6). Applies to
 `"cluster"`/`"cl2"` only.}
+
+\item{cl2_adjustment}{CL2 leverage adjustment (`"auto"`, `"exact"`, or
+`"shortcut"`), used only for `type = "cl2"`.}
 }
 \value{
 A covariance matrix (with CL2 leverage attributes for `type =
diff --git a/man/pffr_vcov.Rd b/man/pffr_vcov.Rd
index 051bc665..697cbe58 100644
--- a/man/pffr_vcov.Rd
+++ b/man/pffr_vcov.Rd
@@ -12,7 +12,8 @@ pffr_vcov(
   dof_correction = NULL,
   edf_type = NULL,
   b2 = TRUE,
-  center_scores = FALSE
+  center_scores = FALSE,
+  cl2_adjustment = NULL
 )
 }
 \arguments{
@@ -39,6 +40,9 @@ recomputation and are never served from or written to the cache.}
 \item{center_scores}{Internal ablation switch (default `FALSE` = current
 behavior). When `TRUE`, center the per-cluster score sums before forming
 the cluster/CL2 meat (X6). Same cache semantics as `b2`.}
+
+\item{cl2_adjustment}{CL2 leverage adjustment. `NULL` (default) inherits
+the fitted choice; otherwise one of `"auto"`, `"exact"`, or `"shortcut"`.}
 }
 \value{
 A covariance matrix.
diff --git a/man/resolve_cl2_adjustment.Rd b/man/resolve_cl2_adjustment.Rd
new file mode 100644
index 00000000..bdc3c084
--- /dev/null
+++ b/man/resolve_cl2_adjustment.Rd
@@ -0,0 +1,35 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/pffr-core.R
+\name{resolve_cl2_adjustment}
+\alias{resolve_cl2_adjustment}
+\title{Resolve the CL2 leverage adjustment}
+\usage{
+resolve_cl2_adjustment(
+  cl2_adjustment = c("auto", "exact", "shortcut"),
+  G,
+  maxDg,
+  p
+)
+}
+\arguments{
+\item{cl2_adjustment}{One of `"auto"`, `"exact"`, or `"shortcut"`.}
+
+\item{G}{Number of clusters.}
+
+\item{maxDg}{Largest cluster size.}
+
+\item{p}{Coefficient-space dimension.}
+}
+\value{
+`"exact"` or `"shortcut"`.
+}
+\description{
+Exact CL2 is the default when the finite-sample correction can matter
+(`G <= 100`) and the dense block multiplication is affordable. Its
+conservative cost proxy `G * maxDg * p^2` covers the dominant
+`(Xw_g Vp) (Xw'Xw)` multiplication; 5e8 scalar operations keeps this below
+the cost of a typical pffr fit on the supported auto-CL2 grids. Large-G
+CL2 uses the historical shortcut because the correction tends to zero and
+the exact block has no material finite-sample benefit there.
+}
+\keyword{internal}
diff --git a/tests/testthat/test-pffr-exactcl2.R b/tests/testthat/test-pffr-exactcl2.R
new file mode 100644
index 00000000..0b2db55f
--- /dev/null
+++ b/tests/testthat/test-pffr-exactcl2.R
@@ -0,0 +1,60 @@
+# A18 exact-CL2 hardening gate. The standalone script writes the same table
+# for release records; this fast test keeps its family/design coverage in CI.
+
+source(testthat::test_path("..", "..", "hardening-exactcl2.R"), local = TRUE)
+
+test_that("exact CL2 is robust on high-leverage family and basis fixtures", {
+  results <- run_exactcl2_hardening()
+
+  expect_equal(nrow(results), 16L)
+  expect_true(all(is.finite(results$n_adjusted)))
+  expect_gt(sum(results$n_adjusted), 0)
+  expect_true(all(is.finite(results$min_block_eig)))
+  expect_true(all(is.finite(results$max_block_kappa)))
+  expect_true(all(is.finite(results$max_abs_exact_nocap)))
+  expect_true(all(results$se_finite_positive))
+  expect_true(all(results$se_ratio_sane))
+})
+
+test_that("CL2 adjustment selection can be forced and is recorded", {
+  fixture <- make_exactcl2_fixture("gaussian", 4L, "influential")
+  fit_auto <- suppressWarnings(suppressMessages(pffr(
+    Y ~ xlin,
+    data = fixture$data,
+    yind = fixture$yind,
+    bs.yindex = list(bs = "ps", k = fixture$k, m = c(2, 1)),
+    sandwich = "cl2"
+  )))
+  fit_exact <- suppressWarnings(suppressMessages(pffr(
+    Y ~ xlin,
+    data = fixture$data,
+    yind = fixture$yind,
+    bs.yindex = list(bs = "ps", k = fixture$k, m = c(2, 1)),
+    sandwich = "cl2",
+    cl2_adjustment = "exact"
+  )))
+  fit_shortcut <- suppressWarnings(suppressMessages(pffr(
+    Y ~ xlin,
+    data = fixture$data,
+    yind = fixture$yind,
+    bs.yindex = list(bs = "ps", k = fixture$k, m = c(2, 1)),
+    sandwich = "cl2",
+    cl2_adjustment = "shortcut"
+  )))
+
+  expect_identical(fit_exact$pffr$sandwich_info$cl2_adjustment, "exact")
+  expect_identical(fit_shortcut$pffr$sandwich_info$cl2_adjustment, "shortcut")
+  expect_identical(fit_auto$pffr$sandwich_info$cl2_adjustment, "exact")
+  expect_equal(fit_auto$pffr$Vsandwich, fit_exact$pffr$Vsandwich)
+  expect_true(is.numeric(fit_exact$pffr$sandwich_info$n_adjusted))
+  expect_true(is.numeric(fit_shortcut$pffr$sandwich_info$n_adjusted))
+})
+
+test_that("automatic exact CL2 uses the documented relevance and cost rule", {
+  resolve <- refund:::resolve_cl2_adjustment
+  expect_identical(resolve("auto", G = 100, maxDg = 50, p = 20), "exact")
+  expect_identical(resolve("auto", G = 101, maxDg = 50, p = 20), "shortcut")
+  expect_identical(resolve("auto", G = 20, maxDg = 500, p = 300), "shortcut")
+  expect_identical(resolve("exact", G = 1000, maxDg = 5000, p = 500), "exact")
+  expect_identical(resolve("shortcut", G = 2, maxDg = 2, p = 2), "shortcut")
+})
diff --git a/tests/testthat/test-pffr.R b/tests/testthat/test-pffr.R
index c9dca72c..37ad8d3a 100644
--- a/tests/testthat/test-pffr.R
+++ b/tests/testthat/test-pffr.R
@@ -1855,17 +1855,18 @@ test_that("CL2 reports leverage diagnostics and warns when the cap is hit", {
     V_benign <- gam_sandwich_cluster_cl2(b, cluster_id, freq = FALSE),
     NA
   )
-  expect_identical(attr(V_benign, "n_capped_clusters"), 0L)
-  expect_true(is.finite(attr(V_benign, "max_leverage")))
-  expect_gt(attr(V_benign, "max_leverage"), 0)
-  expect_lt(attr(V_benign, "max_leverage"), 0.999)
+  expect_identical(attr(V_benign, "cl2_adjustment"), "exact")
+  expect_identical(attr(V_benign, "n_adjusted"), 0L)
+  expect_true(is.finite(attr(V_benign, "min_block_eig")))
+  expect_gt(attr(V_benign, "min_block_eig"), 0)
 
   expect_warning(
     V_capped <- gam_sandwich_cluster_cl2(
       b,
       cluster_id,
       freq = FALSE,
-      leverage_cap = 0.001
+      leverage_cap = 0.001,
+      cl2_adjustment = "shortcut"
     ),
     "CL2 leverage adjustment hit the leverage cap"
   )
@@ -1876,7 +1877,8 @@ test_that("CL2 reports leverage diagnostics and warns when the cap is hit", {
     b,
     cluster_id,
     freq = FALSE,
-    leverage_cap = 0.001
+    leverage_cap = 0.001,
+    cl2_adjustment = "shortcut"
   ))
   expect_equal(
     as.matrix(V_capped),

From 33d8fad840fcc481347efc3b809faa126305b65c Mon Sep 17 00:00:00 2001
From: fabian-s <fabian.scheipl@googlemail.com>
Date: Tue, 21 Jul 2026 11:55:34 +0200
Subject: [PATCH 2/5] auto-rule cost cap 5e8 -> 5e9, calibrated by measurement

Both adjustments pay the same per-cluster eigendecomposition, so the exact
block's marginal cost is only the extra multiplications: timed at ~0.03 s
extra at proxy 5e8 and ~0.7 s total at 5e9 (single-thread BLAS). The old
cap could deny a sub-second correctness upgrade in the saturated small-G
regime where Study EX measured the largest gain (+1.3pp at G=20). G <= 100
relevance bound unchanged (EX-validated range; gain < 0.15pp beyond).
Boundary tests updated; 20 PASS.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01WFfnErmWi9NCGCCJ8YwE95
---
 EXACTCL2-IMPLEMENTATION-NOTES.md    | 10 ++++++++++
 NEWS.md                             |  2 +-
 R/pffr-core.R                       | 19 ++++++++++++-------
 man/resolve_cl2_adjustment.Rd       | 17 +++++++++++------
 tests/testthat/test-pffr-exactcl2.R |  5 ++++-
 5 files changed, 38 insertions(+), 15 deletions(-)

diff --git a/EXACTCL2-IMPLEMENTATION-NOTES.md b/EXACTCL2-IMPLEMENTATION-NOTES.md
index 940bdca4..0ae36851 100644
--- a/EXACTCL2-IMPLEMENTATION-NOTES.md
+++ b/EXACTCL2-IMPLEMENTATION-NOTES.md
@@ -82,3 +82,13 @@ Git returned `Unable to create .../.git/worktrees/refund-wt-s1/index.lock:
 Read-only file system`. The branch `exact-cl2-default` exists and all changes
 remain unstaged in this worktree; no unrelated `tests/testthat/Rplots.pdf`
 change was touched.
+
+## Addendum 2026-07-21 (PI-side review): cost cap raised 5e8 -> 5e9
+
+Timing at the shipped cap showed the exact block's marginal cost over the
+shortcut is ~0.03 s (both paths pay the same per-cluster eigendecomposition;
+the marginal cost is only the extra multiplications), and ~0.7 s total at
+10x the cap. 5e8 could therefore deny a sub-second correctness upgrade in
+the saturated small-G regime where Study EX measured the largest gain
+(+1.3pp at G=20). Cap raised to 5e9; G <= 100 unchanged (EX-validated
+range; gain < 0.15pp beyond).
diff --git a/NEWS.md b/NEWS.md
index 84073aa3..a0e4f25f 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -53,7 +53,7 @@
 * CL2 now uses the exact Bell--McCaffrey leverage block by default when the
   finite-sample correction is relevant and the dense-block calculation is
   affordable (at most 100 clusters and cost proxy `G * max(D_g) * p^2 <=
-  5e8`). At larger/prohibitive problems it uses the historical shortcut.
+  5e9`, a measured sub-second marginal cost). At larger/prohibitive problems it uses the historical shortcut.
   `cl2_adjustment = "exact"` or `"shortcut"` on `pffr()` and `coef.pffr()`
   explicitly selects either variant. Exact-CL2 diagnostics include the number
   of eigenvalue-floored blocks (`n_adjusted`).
diff --git a/R/pffr-core.R b/R/pffr-core.R
index e225eba0..d3712dc5 100644
--- a/R/pffr-core.R
+++ b/R/pffr-core.R
@@ -1804,12 +1804,17 @@ gam_sandwich_cluster_cl2 <- function(
 #' Resolve the CL2 leverage adjustment
 #'
 #' Exact CL2 is the default when the finite-sample correction can matter
-#' (`G <= 100`) and the dense block multiplication is affordable. Its
-#' conservative cost proxy `G * maxDg * p^2` covers the dominant
-#' `(Xw_g Vp) (Xw'Xw)` multiplication; 5e8 scalar operations keeps this below
-#' the cost of a typical pffr fit on the supported auto-CL2 grids. Large-G
-#' CL2 uses the historical shortcut because the correction tends to zero and
-#' the exact block has no material finite-sample benefit there.
+#' (`G <= 100`; Study EX measured the exact-over-shortcut coverage gain
+#' decaying from ~+1pp at G = 20 to under +0.15pp at G = 100) and the dense
+#' block multiplication is affordable. The cost proxy `G * maxDg * p^2`
+#' covers the dominant `(Xw_g Vp) (Xw'Xw)` multiplication. The 5e9 bound is
+#' calibrated by measurement, not guessed: both adjustments pay the same
+#' per-cluster eigendecomposition, so the exact block's MARGINAL cost is only
+#' these multiplications -- timed at ~0.03 s extra at proxy 5e8 and ~0.7 s
+#' total at 5e9 (single-thread BLAS), i.e. the bound caps the correctness
+#' upgrade at well under a second and in particular never denies it in the
+#' saturated small-G regime where its benefit is largest. Large-G CL2 uses
+#' the historical shortcut because the correction tends to zero there.
 #'
 #' @param cl2_adjustment One of `"auto"`, `"exact"`, or `"shortcut"`.
 #' @param G Number of clusters.
@@ -1828,7 +1833,7 @@ resolve_cl2_adjustment <- function(
 
   cost <- G * maxDg * p^2
   if (
-    is.finite(G) && is.finite(maxDg) && is.finite(p) && G <= 100 && cost <= 5e8
+    is.finite(G) && is.finite(maxDg) && is.finite(p) && G <= 100 && cost <= 5e9
   ) {
     "exact"
   } else {
diff --git a/man/resolve_cl2_adjustment.Rd b/man/resolve_cl2_adjustment.Rd
index bdc3c084..009b4c2f 100644
--- a/man/resolve_cl2_adjustment.Rd
+++ b/man/resolve_cl2_adjustment.Rd
@@ -25,11 +25,16 @@ resolve_cl2_adjustment(
 }
 \description{
 Exact CL2 is the default when the finite-sample correction can matter
-(`G <= 100`) and the dense block multiplication is affordable. Its
-conservative cost proxy `G * maxDg * p^2` covers the dominant
-`(Xw_g Vp) (Xw'Xw)` multiplication; 5e8 scalar operations keeps this below
-the cost of a typical pffr fit on the supported auto-CL2 grids. Large-G
-CL2 uses the historical shortcut because the correction tends to zero and
-the exact block has no material finite-sample benefit there.
+(`G <= 100`; Study EX measured the exact-over-shortcut coverage gain
+decaying from ~+1pp at G = 20 to under +0.15pp at G = 100) and the dense
+block multiplication is affordable. The cost proxy `G * maxDg * p^2`
+covers the dominant `(Xw_g Vp) (Xw'Xw)` multiplication. The 5e9 bound is
+calibrated by measurement, not guessed: both adjustments pay the same
+per-cluster eigendecomposition, so the exact block's MARGINAL cost is only
+these multiplications -- timed at ~0.03 s extra at proxy 5e8 and ~0.7 s
+total at 5e9 (single-thread BLAS), i.e. the bound caps the correctness
+upgrade at well under a second and in particular never denies it in the
+saturated small-G regime where its benefit is largest. Large-G CL2 uses
+the historical shortcut because the correction tends to zero there.
 }
 \keyword{internal}
diff --git a/tests/testthat/test-pffr-exactcl2.R b/tests/testthat/test-pffr-exactcl2.R
index 0b2db55f..d4faf7af 100644
--- a/tests/testthat/test-pffr-exactcl2.R
+++ b/tests/testthat/test-pffr-exactcl2.R
@@ -54,7 +54,10 @@ test_that("automatic exact CL2 uses the documented relevance and cost rule", {
   resolve <- refund:::resolve_cl2_adjustment
   expect_identical(resolve("auto", G = 100, maxDg = 50, p = 20), "exact")
   expect_identical(resolve("auto", G = 101, maxDg = 50, p = 20), "shortcut")
-  expect_identical(resolve("auto", G = 20, maxDg = 500, p = 300), "shortcut")
+  # 9e8 ops: below the 5e9 cost cap (raised from 5e8, see notes 2026-07-21)
+  expect_identical(resolve("auto", G = 20, maxDg = 500, p = 300), "exact")
+  # 6.4e10 ops: above the cap even at eligible G
+  expect_identical(resolve("auto", G = 100, maxDg = 1000, p = 800), "shortcut")
   expect_identical(resolve("exact", G = 1000, maxDg = 5000, p = 500), "exact")
   expect_identical(resolve("shortcut", G = 2, maxDg = 2, p = 2), "shortcut")
 })

From 9a4e3c8b3f250b2dd1258dd631667a98d9a62fef Mon Sep 17 00:00:00 2001
From: fabian-s <fabian.scheipl@googlemail.com>
Date: Tue, 21 Jul 2026 12:21:05 +0200
Subject: [PATCH 3/5] review fixes: stale 5e8 docs -> 5e9; wire cl2_adjustment
 into metadata call

Independent review (+ Copilot comments 1-3, all valid): the cap raise in
33d8fad8 missed the pffr.R roxygen (and its generated Rd) and the notes
file's Default-rule body. All now 5e9. Also pass cl2_adjustment through to
pffr_build_metadata() - previously the parameter existed but the sole
caller never supplied it (the field was authoritative only after
apply_sandwich_correction overwrote it).

Copilot comment 4 (use max(abs(eigenvalues)) for kappa) NOT actioned:
B_g is a principal submatrix of the PSD matrix (I-H)^2, hence PSD
(verified numerically, min block eigenvalue +0.23), so abs(max(v)) ==
max(abs(v)) and the diagnostic is already the true condition number.

Review also verified the exact-block identity against a dense (I-H)^2
reference to 4.4e-16 incl. cross-cluster terms. Tests: 20/20.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01WFfnErmWi9NCGCCJ8YwE95
---
 EXACTCL2-IMPLEMENTATION-NOTES.md | 2 +-
 R/pffr.R                         | 6 ++++--
 man/pffr.Rd                      | 3 ++-
 3 files changed, 7 insertions(+), 4 deletions(-)

diff --git a/EXACTCL2-IMPLEMENTATION-NOTES.md b/EXACTCL2-IMPLEMENTATION-NOTES.md
index 0ae36851..84710e00 100644
--- a/EXACTCL2-IMPLEMENTATION-NOTES.md
+++ b/EXACTCL2-IMPLEMENTATION-NOTES.md
@@ -26,7 +26,7 @@ resolved adjustment plus `n_adjusted`, `min_block_eig`, and
 Within the existing CL2 path only, `"auto"` selects exact CL2 when:
 
 * `G <= 100`; and
-* `G * max(D_g) * p^2 <= 5e8`.
+* `G * max(D_g) * p^2 <= 5e9` (raised from 5e8, see Addendum below).
 
 The first bound confines exact CL2 to the small-to-moderate-cluster regime
 where its finite-sample correction is relevant. The second is a conservative
diff --git a/R/pffr.R b/R/pffr.R
index a85fb987..8f71ce73 100644
--- a/R/pffr.R
+++ b/R/pffr.R
@@ -225,7 +225,8 @@
 #' @param cl2_adjustment Leverage adjustment within \code{sandwich = "cl2"}:
 #'   \code{"auto"} (default) uses exact CL2 when \eqn{G \le 100} and the
 #'   dense-block cost proxy \eqn{G\max_g D_g p^2} is at most \eqn{5\times
-#'   10^8}; it otherwise falls back to the historical shortcut. The exact block
+#'   10^9} (a measured sub-second marginal cost); it otherwise falls back to
+#'   the historical shortcut. The exact block
 #'   is \eqn{B_g = I - 2H_{gg} + (H_t^2)_{gg}}, with eigenvalues floored at
 #'   \eqn{(1 - 0.999)^2}. \code{"exact"} and \code{"shortcut"} force either
 #'   variant. Ignored unless the resolved sandwich is \code{"cl2"}.
@@ -496,7 +497,8 @@ pffr <- function(
     ydata = prep$ydata,
     sandwich = sandwich,
     dof_correction = dof_correction,
-    edf_type = edf_type
+    edf_type = edf_type,
+    cl2_adjustment = cl2_adjustment
   )
 
   m <- pffr_attach_metadata(m, prep$algorithm, ret)
diff --git a/man/pffr.Rd b/man/pffr.Rd
index 769aae24..1521b53f 100644
--- a/man/pffr.Rd
+++ b/man/pffr.Rd
@@ -130,7 +130,8 @@ the penalized hat), \code{"edf2"} (mgcv's bias-corrected EDF), or
 \item{cl2_adjustment}{Leverage adjustment within \code{sandwich = "cl2"}:
 \code{"auto"} (default) uses exact CL2 when \eqn{G \le 100} and the
 dense-block cost proxy \eqn{G\max_g D_g p^2} is at most \eqn{5\times
-10^8}; it otherwise falls back to the historical shortcut. The exact block
+10^9} (a measured sub-second marginal cost); it otherwise falls back to
+the historical shortcut. The exact block
 is \eqn{B_g = I - 2H_{gg} + (H_t^2)_{gg}}, with eigenvalues floored at
 \eqn{(1 - 0.999)^2}. \code{"exact"} and \code{"shortcut"} force either
 variant. Ignored unless the resolved sandwich is \code{"cl2"}.}

From 39b0a356c97cb974337a5baedb9a2a65463969ca Mon Sep 17 00:00:00 2001
From: fabian-s <fabian.scheipl@googlemail.com>
Date: Tue, 21 Jul 2026 12:30:43 +0200
Subject: [PATCH 4/5] move implementation notes + hardening results out of the
 package tree

Both are PR documentation, not package files (non-standard top-level files
would also trip R CMD check); their full content is attached to PR #121 as
a comment. The hardening *script* stays on the branch so the results
remain regenerable.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01WFfnErmWi9NCGCCJ8YwE95
---
 EXACTCL2-IMPLEMENTATION-NOTES.md | 94 --------------------------------
 hardening-exactcl2-results.csv   | 17 ------
 2 files changed, 111 deletions(-)
 delete mode 100644 EXACTCL2-IMPLEMENTATION-NOTES.md
 delete mode 100644 hardening-exactcl2-results.csv

diff --git a/EXACTCL2-IMPLEMENTATION-NOTES.md b/EXACTCL2-IMPLEMENTATION-NOTES.md
deleted file mode 100644
index 84710e00..00000000
--- a/EXACTCL2-IMPLEMENTATION-NOTES.md
+++ /dev/null
@@ -1,94 +0,0 @@
-# Exact-CL2 implementation notes
-
-Date: 2026-07-21
-
-## Change made
-
-`gam_sandwich_cluster_cl2()` now implements the exact Bell--McCaffrey
-leverage block
-
-```
-B_g = I - 2 H_gg + (H_t^2)_gg
-```
-
-where `(H_t^2)_gg` is evaluated as
-`(Xw_g Vp) (Xw' Xw) (Xw_g Vp)'`. Its eigendecomposition floors eigenvalues at
-`(1 - leverage_cap)^2`; the historical shortcut continues to cap the
-eigenvalues of `H_gg` and uses `(I - H_gg)^(-1/2)`.
-
-`pffr()` and `coef.pffr()` have a new `cl2_adjustment` argument:
-`"auto"` (default), `"exact"`, or `"shortcut"`. Fit metadata records the
-resolved adjustment plus `n_adjusted`, `min_block_eig`, and
-`max_block_kappa` for exact CL2. Legacy shortcut diagnostics remain available.
-
-## Default rule
-
-Within the existing CL2 path only, `"auto"` selects exact CL2 when:
-
-* `G <= 100`; and
-* `G * max(D_g) * p^2 <= 5e9` (raised from 5e8, see Addendum below).
-
-The first bound confines exact CL2 to the small-to-moderate-cluster regime
-where its finite-sample correction is relevant. The second is a conservative
-dense-matrix proxy for the per-block multiplication
-`(Xw_g Vp) (Xw'Xw)`. Outside either bound the package visibly records and uses
-the shortcut. The package-level `sandwich = "auto"` policy is unchanged.
-
-## Tests run
-
-All R invocations used `OPENBLAS_NUM_THREADS=1 OMP_NUM_THREADS=1` and
-`nice -n 15`.
-
-* `Rscript hardening-exactcl2.R` — passed; wrote
-  `hardening-exactcl2-results.csv`.
-* Targeted test files run with `NOT_CRAN=true` — passed:
-  `test-pffr-exactcl2.R`, `test-pffr-s2-autopolicy.R`,
-  `test-pffr-sandwich-storage.R`, `test-pffr-sandwich-refit.R`,
-  `test-pffr-satterthwaite.R`, and `test-pffr.R`.
-* `devtools::load_all()` — passed.
-* `devtools::document()` — completed and regenerated the affected Rd files.
-  It printed pre-existing roxygen warnings about unrelated tags/method exports
-  in `pffr-core.R`, `poridge.R`, `pffr-methods.R`, and plot methods.
-
-## Hardening output summary
-
-The 16-cell table covers families Gaussian, Poisson, binomial, and Gamma;
-small/moderate response bases `k = 4, 7`; and unbalanced and
-influential-cluster designs. Each has five unbalanced clusters with curve
-counts `1, 1, 2, 3, 5`; the influential design gives one curve an extreme
-covariate value while retaining an ordinary-range outcome.
-
-Actual results:
-
-* All 16 exact covariances and coefficient SE vectors were finite and strictly
-  positive.
-* The median exact/shortcut SE ratio ranged from 1.00 to 1.02.
-* Three influential Poisson/Gamma blocks hit the shipped exact floor
-  (`n_adjusted = 1`); all other cells reported zero adjusted blocks.
-* The smallest observed pre-floor block eigenvalue was `3.20e-07`.
-* Exact versus the `cl2_exact_nocap` analogue (floor `tol = 1e-8`) differed by
-  at most `8.76e-06`; the nonzero differences occur in the floored cells.
-* The largest reported block condition proxy was `2.98e+06`.
-
-## Open questions
-
-The pointwise Satterthwaite degrees-of-freedom routine remains the established
-working-iid approximation based on the shortcut adjustment; this change
-updates the CL2 covariance block only. Exact Bell--McCaffrey df remains a
-separate methodological task, as already documented in the package.
-
-The requested local commits could not be created in this execution environment:
-Git returned `Unable to create .../.git/worktrees/refund-wt-s1/index.lock:
-Read-only file system`. The branch `exact-cl2-default` exists and all changes
-remain unstaged in this worktree; no unrelated `tests/testthat/Rplots.pdf`
-change was touched.
-
-## Addendum 2026-07-21 (PI-side review): cost cap raised 5e8 -> 5e9
-
-Timing at the shipped cap showed the exact block's marginal cost over the
-shortcut is ~0.03 s (both paths pay the same per-cluster eigendecomposition;
-the marginal cost is only the extra multiplications), and ~0.7 s total at
-10x the cap. 5e8 could therefore deny a sub-second correctness upgrade in
-the saturated small-G regime where Study EX measured the largest gain
-(+1.3pp at G=20). Cap raised to 5e9; G <= 100 unchanged (EX-validated
-range; gain < 0.15pp beyond).
diff --git a/hardening-exactcl2-results.csv b/hardening-exactcl2-results.csv
deleted file mode 100644
index 39f55aab..00000000
--- a/hardening-exactcl2-results.csv
+++ /dev/null
@@ -1,17 +0,0 @@
-"family","k","design","G","max_cluster_size","n_adjusted","min_block_eig","max_block_kappa","max_abs_exact_nocap","median_se_ratio_to_shortcut","se_finite_positive","se_ratio_sane"
-"gaussian",4,"unbalanced",5,50,0,0.453590234808712,2.20463299969789,0,1.01263947332551,TRUE,TRUE
-"poisson",4,"unbalanced",5,50,0,0.348940004869249,2.82730371704937,0,1.0079082080914,TRUE,TRUE
-"binomial",4,"unbalanced",5,50,0,0.552153568687061,1.81109035005941,0,1.0010220332576,TRUE,TRUE
-"Gamma",4,"unbalanced",5,50,0,0.332486151450086,3.00764406468859,0,1.00903670061469,TRUE,TRUE
-"gaussian",7,"unbalanced",5,50,0,0.374043657671271,2.67348471091803,0,1.0189182605098,TRUE,TRUE
-"poisson",7,"unbalanced",5,50,0,0.301435120875162,3.31746346310504,0,1.00000000257308,TRUE,TRUE
-"binomial",7,"unbalanced",5,50,0,0.478396994803052,2.09031413420915,0,1.00000000004754,TRUE,TRUE
-"Gamma",7,"unbalanced",5,50,0,0.386030967291895,2.590465752049,0,1.00648632821506,TRUE,TRUE
-"gaussian",4,"influential",5,50,0,1.39826045642e-06,709818.235108635,0,1.00791990080075,TRUE,TRUE
-"poisson",4,"influential",5,50,1,3.20272329212834e-07,2982577.58727008,7.06659823389664e-06,1.01755206252338,TRUE,TRUE
-"binomial",4,"influential",5,50,0,1.58358505475409e-06,623867.247044085,0,1.01039988193779,TRUE,TRUE
-"Gamma",4,"influential",5,50,0,2.35791791958206e-06,423665.797661914,0,1.01294453124424,TRUE,TRUE
-"gaussian",7,"influential",5,50,0,1.12859513423658e-06,882176.398770039,0,1.01074965381931,TRUE,TRUE
-"poisson",7,"influential",5,50,1,7.21418514260108e-07,1385176.29210241,8.75780657148509e-07,1.01500324062661,TRUE,TRUE
-"binomial",7,"influential",5,50,0,1.35355708907881e-06,738794.066446026,0,1.00000000006042,TRUE,TRUE
-"Gamma",7,"influential",5,50,1,7.39633434387876e-07,1348528.03666834,4.54114080495374e-07,1.00490798665056,TRUE,TRUE

From 276ba7179e5db2702cea3de9110b3d99ec0d9484 Mon Sep 17 00:00:00 2001
From: fabian-s <fabian.scheipl@googlemail.com>
Date: Tue, 21 Jul 2026 12:52:40 +0200
Subject: [PATCH 5/5] remove the hardening script from the package tree as well

Non-standard top-level file; full script preserved as PR #121
documentation alongside the notes and results it generates.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01WFfnErmWi9NCGCCJ8YwE95
---
 hardening-exactcl2.R | 141 -------------------------------------------
 1 file changed, 141 deletions(-)
 delete mode 100644 hardening-exactcl2.R

diff --git a/hardening-exactcl2.R b/hardening-exactcl2.R
deleted file mode 100644
index 96a10d1e..00000000
--- a/hardening-exactcl2.R
+++ /dev/null
@@ -1,141 +0,0 @@
-# Scoped exact-CL2 hardening gate (DESIGN amendment A18).
-#
-# Run from the package root with:
-# OPENBLAS_NUM_THREADS=1 OMP_NUM_THREADS=1 nice -n 15 Rscript hardening-exactcl2.R
-# The script writes hardening-exactcl2-results.csv in the package root.
-
-make_exactcl2_fixture <- function(family_name, k, design) {
-  set.seed(
-    match(family_name, c("gaussian", "poisson", "binomial", "Gamma")) *
-      100 +
-      k +
-      match(design, c("unbalanced", "influential"))
-  )
-  n_curve <- 12L
-  n_grid <- 10L
-  grid <- seq(0, 1, length.out = n_grid)
-  xlin <- stats::rnorm(n_curve)
-  eta <- outer(xlin, rep(0.35, n_grid)) +
-    matrix(rep(0.2 * sin(2 * pi * grid), each = n_curve), n_curve, n_grid) -
-    0.4
-  Y <- switch(
-    family_name,
-    gaussian = eta + matrix(stats::rnorm(n_curve * n_grid, sd = 0.35), n_curve),
-    poisson = matrix(
-      stats::rpois(n_curve * n_grid, lambda = exp(pmin(eta, 2))),
-      n_curve
-    ),
-    binomial = matrix(
-      stats::rbinom(n_curve * n_grid, 1, stats::plogis(eta)),
-      n_curve
-    ),
-    Gamma = matrix(
-      stats::rgamma(n_curve * n_grid, shape = 5, scale = exp(pmin(eta, 2)) / 5),
-      n_curve
-    )
-  )
-  # Keep the response in its ordinary range, then make one covariate curve
-  # extreme. This isolates high leverage from a separable/outcome-saturation
-  # failure and exercises the shipped eigenvalue floor (especially Gaussian,
-  # Poisson, and Gamma at k = 9).
-  if (design == "influential") xlin[1] <- 3000
-  # Five deliberately unbalanced independent clusters, with cluster 1 made
-  # influential above. Their curve counts 1, 1, 2, 3, 5 yield unequal blocks.
-  cluster <- rep(seq_len(5), times = c(1, 1, 2, 3, 5))
-  list(
-    data = data.frame(Y = I(Y), xlin = xlin),
-    yind = grid,
-    cluster = cluster,
-    family = switch(
-      family_name,
-      gaussian = stats::gaussian(),
-      poisson = stats::poisson(),
-      binomial = stats::binomial(),
-      Gamma = stats::Gamma(link = "log")
-    ),
-    k = k
-  )
-}
-
-run_exactcl2_hardening <- function() {
-  cases <- expand.grid(
-    family = c("gaussian", "poisson", "binomial", "Gamma"),
-    k = c(4L, 7L),
-    design = c("unbalanced", "influential"),
-    stringsAsFactors = FALSE
-  )
-  rows <- vector("list", nrow(cases))
-
-  for (i in seq_len(nrow(cases))) {
-    case <- cases[i, ]
-    fixture <- make_exactcl2_fixture(case$family, case$k, case$design)
-    fit <- suppressWarnings(suppressMessages(refund::pffr(
-      Y ~ xlin,
-      data = fixture$data,
-      yind = fixture$yind,
-      family = fixture$family,
-      bs.yindex = list(bs = "ps", k = fixture$k, m = c(2, 1)),
-      sandwich = "none"
-    )))
-
-    # `1 - sqrt(tol)` makes the exact floor equal tol: Study EX's
-    # cl2_exact_nocap comparator without duplicating the package implementation.
-    V_exact <- suppressWarnings(refund:::pffr_vcov(
-      fit,
-      sandwich = "cl2",
-      cluster = fixture$cluster,
-      cl2_adjustment = "exact"
-    ))
-    V_nocap <- suppressWarnings(refund:::gam_sandwich_cluster_cl2(
-      refund:::pffr_model_based_gam(fit),
-      refund:::build_cluster_id(fit$pffr, cluster = fixture$cluster),
-      cl2_adjustment = "exact",
-      leverage_cap = 1 - sqrt(1e-8)
-    ))
-    V_shortcut <- suppressWarnings(refund:::pffr_vcov(
-      fit,
-      sandwich = "cl2",
-      cluster = fixture$cluster,
-      cl2_adjustment = "shortcut"
-    ))
-
-    se_exact <- sqrt(pmax(diag(V_exact), 0))
-    se_shortcut <- sqrt(pmax(diag(V_shortcut), 0))
-    positive <- is.finite(se_exact) &
-      is.finite(se_shortcut) &
-      se_exact > 0 &
-      se_shortcut > 0
-    ratio <- se_exact[positive] / se_shortcut[positive]
-    rows[[i]] <- data.frame(
-      family = case$family,
-      k = case$k,
-      design = case$design,
-      G = length(unique(fixture$cluster)),
-      max_cluster_size = max(table(fixture$cluster)) * length(fixture$yind),
-      n_adjusted = attr(V_exact, "n_adjusted"),
-      min_block_eig = attr(V_exact, "min_block_eig"),
-      max_block_kappa = attr(V_exact, "max_block_kappa"),
-      max_abs_exact_nocap = max(abs(V_exact - V_nocap)),
-      median_se_ratio_to_shortcut = stats::median(ratio),
-      se_finite_positive = all(positive),
-      se_ratio_sane = length(ratio) > 0 &&
-        min(ratio) > 1e-3 &&
-        max(ratio) < 1e3,
-      stringsAsFactors = FALSE
-    )
-  }
-  do.call(rbind, rows)
-}
-
-if (sys.nframe() == 0L) {
-  devtools::load_all(quiet = TRUE)
-  results <- run_exactcl2_hardening()
-  utils::write.csv(results, "hardening-exactcl2-results.csv", row.names = FALSE)
-  print(results)
-  stopifnot(
-    all(results$se_finite_positive),
-    all(results$se_ratio_sane),
-    all(is.finite(results$n_adjusted)),
-    any(results$n_adjusted > 0)
-  )
-}