Differentiable SCETlib prediction as a rabbit parameter model

[lucalavezzo]

What

A new package wremnants/postprocessing/scetlib_ad/ providing a rabbit ParamModel
in which every theory parameter SCETlib exposes is a continuous fit parameter with
exact derivatives, instead of a discrete template morph whose joint response with
the others is an outer product.

parameter group	count (Z)	cache/runcard flag it needs
alphaS	1	has_as for the PDF-consistent version (see below)
nonperturbative λ (Collins–Soper + TMD)	8	—
theory nuisance parameters	10	a non-off [TNPs] block
resummation kappa_R	1	diff_scales
matching transition points x1..x3	3	diff_scales — sign-inverted, see Validation
kappa_F	1	diff_scales and has_muf (inert without it)
PDF eigenvector coefficients	n_eig	n_eig > 0

Which of these exist is a property of the cache, not of this code: the model
reads gradient_param_names() and registers what it finds.

The prediction comes from the SCETlib autodiff-sigmaul branch, whose
ScetlibCachedXsecTF replays a prepared cache (compressed bin rules for the
resummed piece + a frozen fixed-order grid for the nonsingular) and returns exact
first and second derivatives from clad.

How it works

1  SCETlib cached rule replay  ->  sigma(p; g)       boson level, gen grid
2  fold through the response R ->  sigma_reco(p; b)  gen -> reco
3  ratio to the reference      ->  rnorm(b, proc)    handed to rabbit

There is no surrogate: autodiff differentiates the real prediction.
ScetlibCachedXsecTF is an ordinary TF-differentiable function whose backward pass
is itself a custom_gradient contracting Hessian-vector products, so nested
GradientTapes work and TF drives every C++ call. The model just calls it inside
the graph, exactly as examples/matched_ad/tf_gradients.py does. An earlier
revision injected a local quadratic behind stop_gradient instead; that path was
removed (62fb5881, "one differentiation path, no fallback") — there is now one
path, so there is nothing to keep in sync.

That imposes one requirement which is easy to break by accident: map rabbit's fit
vector into SCETlib's layout with a constant 0/1 matrix multiply, never
tensor_scatter_nd_update. rabbit's vector holds only the fitted parameters, POIs
first, while SCETlib's holds every registered parameter in registry order, so some
mapping is unavoidable. A scatter's backward pass contains a gather, whose gradient
TF represents as tf.IndexedSlices, and the bridge's second-order py_function
payloads call .numpy() on the incoming cotangent and fail on it — so anything past
first order breaks. The matmul is bit-identical (entries are exactly 0 and 1) and
free at these sizes (at most ~25 × 25).

Other notes:

• One rabbit job. Second derivatives come from the bridge's own HVP
  contraction, so the composite Hessian is exact and fit + postfit covariance run
  together. Confirmed on the Asimov fits below: one job, edmval ~ 1e-27.
• --jitCompile off is required (XLA cannot compile a PyFunc) and is
  enforced at construction with an actionable message.
• GenFold sums cache bins onto the card's gen grid — handling a different
  nesting order, a signed-Y cache folded onto |Y|, and a cache finer than the fit
  binning — and verifies every gen bin is exactly tiled rather than assuming it.
  Note it does not normalise the Y convention out of the values it returns: a
  positive-side-only cache yields half the |Y| cross section. That cancels in the
  ratio the fit uses, and does not cancel if you compare absolutely.

Profile scales as unit nuisances

set_diff_scales(1) makes kappa_R and the transition points differentiable. But
their templates encode multiplicative (κ_R: 0.5 / 1 / 2) or asymmetric
(x2: 0.35 / 0.6 / 0.75) steps, which a symmetric Gaussian prior on the physical
value cannot express. So params.REPARAM maps them to unit nuisances —
exp(θ·ln2) for the scales, a quadratic for the transition point — and TF
differentiates the map, so the chain rule keeps gradient and Hessian exact.

The map itself is verified bit-exact (0.000e+00): θ = −1/0/+1 land on the
physical values, and the analytic Jacobian column times dκ/dθ equals TF's AD
gradient per bin, with dκ/dθ = 0.6931471806 = ln2 and 0.2 respectively. (The map
being right is not the same as the underlying scale response being right — see the
transition points below.)

Methodology note for anyone re-testing this: do not validate these
derivatives with a central difference across the anchor. The value surrogate has a
knot there (c_val forces exactness at the anchor), so a symmetric difference
averages two different slopes and the error is flat in h over four decades —
which looks like a real failure and is not. Compare against the analytic Jacobian.

Validation

Against a validated SCETlib production run and the histmaker, on the analysis
binning (Z, ptll × yll, 210 gen bins, cache_aspair):

test	result
cached rules vs a live SCETlib evaluation (small cache)	6e-15 value, 1.5e-15 gradient, 3.4e-14 Hessian
σ_gen vs the production driver, full space	7.4e-06
σ_reco vs the histmaker's corrected reco, shape	0.128% yield-weighted
σ_reco vs indata.norm, absolute (no rescaling)	0.149%, total 0.998847
injection closure (small cache)	exact, every other parameter unmoved, EDM 6e-22

Per-variation response against the Corr templates, all 37 labels the reference
carries (validate_variations.py), as max|dev| / mean|dev|:

block	agreement	verdict
20 TNP variations	2.2e-16 … 7.4e-04 / 2e-18 … 3.3e-05	excellent
8 NP λ variations	4.6e-04 … 4.9e-03 / 9e-06 … 2.1e-04	good
mufdown / mufup	2.9e-03 / 1.4e-02	usable, worth understanding
kappaFO2.-kappaf0.5	4.5e-03	good
kappaFO0.5-kappaf2.	4.0e-02	the κ_R down direction is 10× worse than up
3 transition_points*	1.1e-01 … 2.0e-01	WRONG SIGN — see below

The transition-point directions disagree in sign with their templates. All three
move the prediction the opposite way from the reference:

label                          model rng        ref rng
transition_points0.2_0.35_1.0  [1.0000,1.1593]  [0.9602,1.0000]   model UP,   ref DOWN
transition_points0.2_0.75_1.0  [0.8957,1.0000]  [1.0000,1.0207]   model DOWN, ref UP
transition_points0.3_0.6_0.9   [0.8971,1.0103]  [0.9990,1.0135]   model DOWN, ref UP

The third moves x1/x3 rather than x2, so this is not a mapping slip in the
validation table, and the reparametrisation map is separately verified bit-exact.
The central values match the templates' (x = 0.2, 0.6, 1.0). Leading suspect is a
convention difference between set_diff_scales and the production
transition_points setting. Until it is resolved, resumTransition* must not be
floated for a physics result — and note that the σ(α_s) in fit B below did float
resumTransition2, so its number and its ρ are provisional.

Also worth flagging: the κ_R down direction agreeing only to 4% matters more than it
looks, because κ_R is the direction that dominates σ(α_s).

Asimov reco fits

-t -1, real card, login-node CPU, one job each:

floating	edmval	saturated 2ΔNLL / ndof	wall	σ(α_s)
6 (α_s + 5 λ)	9.7e-28	0.0 / 774	46 min	6.16e-04
8 (+ κ_R, x2 — x2 provisional)	9.7e-28	0.0 / 772	58 min	1.81e-03

Every parameter returns exactly at truth, and the Hessian is finite and sensible.

ρ(α_s, resumScaleMuR) = +0.927. κ_R and α_s both set the strength of the
resummed logs, so they trade off almost freely in the qT shape — floating κ_R as a
continuous nuisance costs a factor 2.9 on σ(α_s). The κ_R treatment is therefore
the dominant choice for the α_s uncertainty, more than any NP λ. (Also
ρ(λ2, λ2_ν) = −0.97, the two low-qT damping knobs being near-degenerate as
expected.) Neither σ(α_s) is quotable: these caches carry no PDF eigenvectors.

Two bugs this validation caught

The AD path initially disagreed with a validated production run by 1–3%. Every
earlier test compared the AD path against itself (exact to 1e-15), so none of them
could see it. Running the same configuration through SCETlib's own production driver
split it into two independent bugs whose product reproduced the discrepancy to 1e-5
in every bin:

1. Ours — the cache builder configured the calculation without
   configure_ew_parameters / configure_fiducial_volumes /
   configure_calculation, so it silently used SCETlib's default EW inputs rather
   than the runcard's. Flat +1.61%. Fixed here.
2. Upstream — DrellYan::operator()'s ad::Node_shared was hoisted out of the
   node loop, so every node reused node 0's shared sub-expressions. qT-dependent,
   ±1.4%. Reported and cherry-picked upstream as b919b61 (which also turned up
   a second occurrence).

After both: A/driver = 0.000e+00, B_cacheON/driver = 3.028e-06.

Cache contents are load-bearing — read this before quoting anything

The model registers what the cache has, and a cache built with --no-pdf is missing
directions silently unless the card still carries the templates:

• has_as = 0 → α_s is a derivative at fixed PDF. The α_s pair rides on the
  existing alphas slot so one parameter moves the calculation and the PDF;
  without it, do not quote α_s. --no-pdf now says so out loud.
• has_muf = 0 → resumScaleMuF has an identically zero derivative. It is
  refused at fit time rather than silently doing nothing.
• n_eig = 0 → no PDF uncertainty from the model, so the card's pdf* templates
  must be kept.

Measured build cost at 210 bins: rules 9.0 min, FO warm 20.6 min, and 19.4 min per
PDF member for the fixed-order variations. The α_s + μF cache (4 members) is
~1.9 h; the full 29 eigenvector pairs (62 members) is ~20 h. The member loop is
the only serial axis left — the bin loop inside it is already _parallel_run — so
that is what sharding would have to split.

Not yet done

• The transition-point sign inversion above. Highest priority; it gates
  floating resumTransition* at all.
• No PDF eigenvector cache yet (the ~20 h build, plus beamfunc grids for 58
  members). Until then the card must keep pdf*.
• TNPs are not floated by default — fit_params excludes resumTNP_*. They
  work when asked for; the default is deliberate, not an oversight.
• Below qT ≈ 2 GeV the model's nonsingular cutoff (0.1 GeV) differs from the
  reference's (--qtCutoff 1.0). It largely cancels in the ratio the fit uses, and
  is being handled separately.
• Real-data fits, impacts, and a comparison of the continuous κ_R treatment against
  the discrete resumFOScale templates it would replace.

Changes to the shared Z dilepton histmaker

Most of this PR is a self-contained package, but five files are shared and one of
them changes a default. Flagging it explicitly so it can be objected to.

--responseGenBinning now defaults to theoryCorr. It adds two histograms to a
Z dilepton unfolding run: the reco x gen response on the theory correction's own gen
grid (770 gen bins for the Z — 70 qT bins to 100 GeV x 11 |Y| bins, |Y| truncated
at the reco yll edge 2.5) and the gen total on that same grid. The grid is read from
the correction file at runtime; no bin count is hardcoded.

Who pays for it, and how much:

• Nothing changes for a default run. The block is inside if args.unfolding:, and
  --unfolding is store_true/False, as is --poiAsNoi. A mz_dilepton run that
  passes neither never reaches the flag.
• The only runs affected are Z dilepton poi-as-noi unfolding histmakers — which are
  the runs that want a response matrix in the first place.
• mz_dilepton is Z-only, so no W histmaker is touched (those use
  mw_with_mu_eta_pt.py).
• Cost: 2.49 MB of an 86.2 MB output (2.9%), measured on
  260826_Z_histmaker_respgrid by summing the two histograms' stored datasets
  (nominal_prefsr_yieldsResponse 2.47 MB, prefsr_response 0.02 MB). They are
  23.5 MB in memory and compress ~9.5x because the reco x gen matrix is very sparse.
• --responseGenBinning none opts out in one flag.

The two preconditions are no longer hard errors when the flag is on only by default:
--poiAsNoi off means there is no reco x gen histogram to put on a finer grid, so the
response is skipped with a warning rather than failing a run that never asked for it.
Passing --responseGenBinning theoryCorr explicitly still raises, so an explicit
request is never silently downgraded.

--responseGenPtVExtend (the qT > 100 extension, 803 bins) stays PROVISIONAL and
off by default: the corrections do not exist on the wider grid yet, so the MC would be
uncorrected there while the model predicts a corrected cross section.

Notes for review

• The package is standalone: no code dependency on the scetlib_np package of
  SCETlib NP continuous parameter model #701. The two scetlib_np strings in response.py are datacard conventions
  (the auxiliary group name and the metadata key existing cards carry) and are
  overridable via response_group=. One optional script,
  compare_to_np_model.py, imports scetlib_np lazily — that script exists to
  cross-check the two models, so it is the one place the dependency is intended.
• Guards for the failure modes that are otherwise silent: cache/card binning
  mismatch, cache anchor vs the card's recorded nonperturbative values, template
  nuisances that would double-count a fitted parameter (regex-based, so ^pdf\d+
  catches the eigenvectors without catching pdfAlphaS), TNPs fitted without
  priors, and any parameter whose Jacobian column is identically zero.
• The diff is 22 files. 17 are the package itself; the other five are shared and
  are called out in "Changes to the shared Z dilepton histmaker" below.
• Linting: the three repo-wide commands CI actually runs (isort --profile black --line-length 88, flake8 --max-line-length 88 with the F-code select list, and
  black --check) are clean in the CI container. pylint is not installed locally,
  so the pre-commit hook's final step could not be run — CI's lint job does not run
  pylint either.

🤖 Generated with Claude Code

https://claude.ai/code/session_01CnJ9YKK8c1q1sCDouM6y1c