Unified constraint structure (cw/x0); ParamModel Gaussian priors#139
Merged
davidwalter2 merged 17 commits intoJun 23, 2026
Conversation
--blindingGroup: regex-defined groups of parameters that share a single deterministic blinding offset (seeded from the group regex string), so relative pulls and differences between matched parameters stay meaningful while absolute values remain blinded. Parameters matching both --unblind and --blindingGroup are a configuration error and abort the fit. match_regexp_params now returns the union of exact and regex matches (deduplicated, order-preserving) instead of short-circuiting on the first exact match, so mixed exact/regex expression lists behave consistently. --paramModelPriors: opt-in Gaussian priors on ParamModel parameters. The Fitter reads optional prior_sigmas / prior_means attributes from the ParamModel and adds the corresponding 0.5*(x-mu)^2/sigma^2 penalty to the NLL constraint term; the applied priors are stored in the output metadata. Off by default; without the flag all ParamModel params float free as before. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Per review on WMass#133 (parsing.py:252): there is no rabbit-side CLA anymore. If a ParamModel declares prior_sigmas, the priors are applied; whether and how to enable them (e.g. a token in the --paramModel spec) is the model's own decision. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Per review on WMass#133: param_prior_sigmas / param_prior_means are stored once as TF constants holding the priors as declared (NaN where no prior); the compute-safe masked forms are derived inside _compute_lc, and the output metadata reads the arrays via .numpy(). Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
CompositeParamModel summed npoi but concatenated params in model order, so a POU-carrying model placed before a POI-carrying one leaked its POUs into the composite's POI slice. Every npoi-sliced consumer then misbehaved: get_poi() squared the POUs (allowNegativeParam=False default) and, on real data, blinded them as "signal strength modifiers" -- visible with --computeSaturatedProjectionTests + a custom npou>0 param model, where the saturated refit ran with squared+blinded theory parameters (frozen ones unable to compensate), invalidating the saturated p-value. Fix: assemble the composite per-block, [poi(m1), poi(m2), ... | pou(m1), pou(m2), ...], and reassemble each submodel's native [poi | pou] vector in compute(). Legacy-valid orderings produce identical layouts. Additionally: - allowNegativeParam is derived from the POI-carrying submodels only (mixed flags raise: the fitter applies a single squaring transform to the POI block). load_models' any() over all models had the same trap. - submodel prior_sigmas/prior_means propagate through the permutation (NaN = no prior), so a composite refit keeps the same prior penalty. - param_impact_groups are merged (name-based, permutation-safe). - is_linear no longer claims linearity for products of >1 parameter-dependent factor or sqrt-stored POIs. - blinding log message: "signal strength modifier" -> "parameter" (review request). Tested in tests/test_composite_param_model.py (layout, compute slicing vs manual evaluation, gradient flow, legacy-ordering invariance, flag derivation/conflicts, prior+group propagation). Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
ParamModel priors now use the same constraint structure as the nuisance constraints: param_cw holds the per-parameter constraint weight (1/sigma^2; 0 = free) and param_x0 the constraint center (the prior mean) — the ParamModel-block analogues of indata.constraintweights and theta0. Consequences: - _compute_lc penalizes cw * 0.5 * (param - x0)^2, same form as the nuisance term. - prefit variances are 1/cw uniformly, so priored params get proper prefit uncertainties (pulls/constraints work without special cases). - Toys treat the prior like any auxiliary constraint: frequentist toys fluctuate param_x0 around its default with width sqrt(1/cw) exactly like theta0, and bayesian toys sample the priored param values from their priors. Previously priors were static in toys, silently underestimating the toy spread of anything correlated with a priored parameter. This also fixes the long-standing FIXMEs in bayesassign / frequentistassign: constraint centers now fluctuate around their defaults (not zero) with width sqrt(1/cw), and unconstrained entries are no longer randomized with unit width. - Gaussian priors on POIs require allowNegativeParam=True; with the squared storage the penalty would silently apply to sqrt(poi), so the Fitter raises instead. Verified on the test tensor: identical postfit optimum as the previous implementation, prior width visible in the prefit variances, per-toy fluctuation of the prior centers, and unchanged behavior for models without priors. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
The prior centers (param_x0) are auxiliary observables like theta0, so they now appear as global-impact sources, one column per priored param (named <param>_prior, appended after the systs / syst groups in the output axes): - likelihood-based global impacts: _compute_global_impacts_x0 already differentiates the full constraint term; the ParamModel-block rows (nonzero where priors are declared) are now kept instead of discarded, in per-1-prefit-sigma units automatically via sqrt(d2lc/dx2) = 1/sigma. Their variance is subtracted from the residual data-stat term, which previously absorbed it. - gaussian global impacts: _dxdvars additionally differentiates with respect to param_x0 and the impacts are dx/dparam_x0 * sigma. - observable (histogram) impacts: both variants extended through _dndvars with the chain-rule term pdndx @ dxdparam_x0. - profiled chi2: the residual covariance in _residuals_profiled gains the prior-center variance contribution. Verified on the test tensor: the likelihood and Gaussian paths agree exactly on the prior impact, and the variance closure sum(sources^2) + stat^2 + binByBinStat^2 = sigma_tot^2 holds to all printed digits. Output schemas are unchanged for fits without priors. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
_compute_lc is now one expression: cw * 0.5 * (x - x0)^2 summed over the
full effective parameter vector [poi, model_nui, theta], with
cw = concat(param_cw, constraintweights) and x0 = concat(param_x0,
theta0). No special-casing of the ParamModel block — a model without
priors simply contributes cw = 0 entries. The full-NLL normalization is
generalized to 0.5*log(2 pi / cw) for constrained entries (identical to
the previous 0.9189*cw for the standard cw in {0, 1}).
var_param_x0 is promoted to an attribute (the analogue of var_theta0)
and the toy randomization of param_x0 is unconditional, mirroring
theta0.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
lucalavezzo
force-pushed
the
split/param-priors-blinding
branch
from
664d055 to
c674aae
Compare
The constraint structure is now held in two full-length objects,
index-aligned with x:
cw constraint weights (1/sigma^2 for priored ParamModel params,
indata.constraintweights for systs; 0 = unconstrained)
x0 constraint centers (prior means / theta0), a tf.Variable
var_x0 prefit center variances (1/cw; 0 where free)
fitter.theta0, param_x0, param_cw, var_theta0 and var_param_x0 are
gone. Every consumer is migrated, in most cases with less code since
the +- nparams index offsets cancel:
- _compute_lc is one term: cw * 0.5 * (x - x0)^2.
- prefit_variance is a single tf.where over cw.
- frequentistassign / bayesassign are single full-vector expressions.
- _dxdvars / _dndvars differentiate w.r.t. x0 once; the dxdtheta0 /
dxdparam0 pair and its plumbing through the gaussian global impacts
collapse into a dxdx0 split helper.
- _residuals_profiled has a single x0-variance term.
- nonprofiled_impacts and global_asym_impacts operate on x0 directly
(x0[idx] instead of theta0[idx - nparams]).
- the saturated-projection re-init preserves the toy-randomized centers
via the x0 syst-block slice.
No behavior change: postfit results, impacts, closures and the
global-asym test suite are identical to the previous two-block
implementation. External code that accessed fitter.theta0 must use
fitter.x0[fitter.param_model.nparams:] instead.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
An expression that names one parameter exactly can no longer also match parameters whose names merely extend it (prefix matches could silently unblind more than intended); families are matched with an explicit pattern, e.g. alphaS.*. The parameters an --unblind expression resolves to are now reported at INFO level. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
The variation magnitude used 1/cw (the variance) instead of 1/sqrt(cw) (the sigma); identical for the 0/1 tensor constraint weights, but wrong for ParamModel priors where cw = 1/sigma^2 is genuinely non-unit. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
The saturated-projection re-init dropped --blindingGroup (falling back to per-name blinding seeds) and reset the ParamModel block of x0 to its defaults, losing toy-fluctuated prior centers. Both are now carried through the composite [POIs | POUs] permutation. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
… default - prior_sigmas/means and the cw/x0 numpy intermediates are built in self.indata.dtype (via as_numpy_dtype) rather than a hardcoded np.float64, matching how the rest of init_fit_parms respects the input dtype (davidwalter2 review on WMass#139). - free (unpriored) entries of the ParamModel x0 block now default to the parameter own default instead of 0. cw = 0 makes this irrelevant to the likelihood, but nonprofiled_impacts reads x0 as a parameter natural center, so a frozen free param is now varied around its default rather than around 0. With no priors declared x0default now equals xdefault. - fix a stale theta0 reference in the global_asym_impacts docstring. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…ored params in asym impacts Item 5 (exact + cheaper): x0 enters the NLL only through cw*(x-x0)^2, so the loss derivative w.r.t. any center with cw = 0 is identically zero. _dxdvars now perturbs only the constraint sources (the nuisance block plus priored params) via a scattered tangent and returns the response split into (dxdtheta0, dxdparam0), shrinking the Jacobian from [npar, npar] to [npar, nsyst (+ n_prior)]; the gaussian global-impacts helpers take the pre-split pieces (dropping _split_x0_sources). _dndvars likewise drops its x0 watch since the forward model has no x0 dependence (verified d(yields)/dx0 = 0 exactly). Validated bit-identical to the previous code on the test tensor: parms gaussian/likelihood impacts, observable impacts, and the profiled-chi2 residual covariance all unchanged. Item 6b: global_asym_impacts now scans priored ParamModel params as sources too (full-x indexed, labelled <param>_prior), matching the source set of gaussian_global_impacts_parms, with ParamModel impact groups folded into the grouped envelopes. The per-source center shift is scaled by the source prefit width sqrt(var_x0) (1 for unit-constrained nuisances, the prior sigma for priored params), so a unit-sigma asym impact stays in the same units as the gaussian one. On the test tensor the sig_prior asym impact (-1.10e-2, +1.11e-2) reproduces the gaussian value (1.106e-2) to <1%, the small spread being genuine non-Gaussianity; nuisance impacts are unchanged. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
lucalavezzo
changed the title
davidwalter2
left a comment
davidwalter2
left a comment
Collaborator
There was a problem hiding this comment.
The params and systs are not really unified in this PR, there are lot's of places where the two are split and/or merged together with specific treatment which isn't getting clear to my why it is needed. Is it possible to unify them better and work with one collection in most cases?
Can the modifications on blinding groups and the bugfix on composite param models each one be moved into separate PRs?
| ) | ||
| self.param_prior_active = True | ||
| self.param_prior_mask = tf.constant(mask, dtype=tf.bool) | ||
| self.param_prior_sigmas = tf.constant( |
Collaborator
There was a problem hiding this comment.
These are duplicated in ParamModel and fitter class which calls for trouble, we should probably only define them in only one of them.
| ) | ||
| t_lws = time.perf_counter() | ||
| dxdtheta0_tf, _, _ = fitter._dxdvars() | ||
| dxdtheta0_tf, dxdparam0_tf, _, _ = fitter._dxdvars() |
Collaborator
There was a problem hiding this comment.
I think theta0 and param0 should be treated on equal footing, so here instead of returning the two dxdtheta0_tf and dxdparam0_tf collection return one dxdx0_tf
| # _x0_constraint_source_idxs. The resulting columns are ordered | ||
| # [theta block | priored params] and split below. | ||
| src_idxs = self._x0_constraint_source_idxs() | ||
| n_src = self.indata.nsyst + ( |
Collaborator
There was a problem hiding this comment.
This looks very error prone, what do we gain from sticking together two collections, I think it is much cleaner to work with a single vector self.x0 which replaces the previous self.theta0 and contains all params and systs.
| # ParamModel prior penalties (rows in the ParamModel block, nonzero where | ||
| # the model declares priors). Per-1-sigma units come out automatically | ||
| # since sc = sqrt(d2lc/dx2) = sqrt(cw) = 1/sigma. | ||
| impacts_x0 = _compute_global_impacts_x0(x, compute_lc_fn, cov_dexpdx) |
Collaborator
There was a problem hiding this comment.
When params and systs are treated on equal footing I think it should be possible to work with impacts_x0 without splitting them into the two groups, or am I missing something?
Per davidwalter2 review: treat the ParamModel params and the systs on equal footing instead of carrying / splitting two collections. - Fitter builds a single x0_source_idxs (the constrained centers, ordered [nuisance block | priored params]) and var_x0_sources, derived from cw. - _dxdvars / _dndvars return ONE dxdx0 / dndx0 over those sources instead of the (dxdtheta0, dxdparam0) pair; the helper methods _x0_constraint_source_idxs and _x0_source_vars are gone. - the likelihood global impacts gather impacts_x0 over source_idxs in one step (no theta-block slice + separate prior gather); the gaussian impacts take one dxdx0 / varx0 and compute the per-source impact uniformly as dxdx0 * sqrt(var) (systs have var = 1, priored params carry their sigma). - the only remaining nuisance-vs-param boundary is a local slice when assembling the grouped axis in its historical order [syst groups | stat | bbb | priors], driven by n_param_sources. - global_asym_impacts warm-start maps each source to its single dxdx0 column. - drop the duplicated param_prior_sigmas / _means / _mask Fitter attributes; priors live only on the model, and the output metadata reads them there. param_prior_active is now derived from cw. Validated bit-identical to the previous code on the test tensor: likelihood and gaussian global impacts (parms + observables), and the asym sig_prior impact still matches the gaussian one to <1%. Source axis unchanged. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
This was referenced Jun 22, 2026
Integrates the merged sibling PRs into the priors/unification branch: WMass#140 fit-ops (minimizer flags, externalPostfit), WMass#142 scalePoiScan, WMass#143 CompositeParamModel layout fix, WMass#144 blinding groups. Conflicts (both trivial): - fitter.py: blinding_group=options.blindingGroup (main's direct form; the SimpleNamespace test fixtures already carry blindingGroup=[]). - helpers.py: kept main's cleaned-up version (dropped the explanatory comment removed on main). Also dropped tests/test_composite_param_model.py, which the merge re-introduced from this branch but was deliberately removed on main. Validated: lint clean; pytest (15) + the five CI script-tests pass; the ParamModel-prior closure is unchanged (gaussian sig_prior 1.105769e-02, asym (-1.097e-2,+1.114e-2)); parsing.py is identical to main so all of WMass#140/WMass#142/WMass#144's CLAs are present. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Collaborator
Author
|
Thanks @davidwalter2 I tried to improve a bit the unification after discussing. Everything should be more folded into a single cw / x0 / var_x0. Can you take a look to see if this is closer to what had in mind? |
lucalavezzo
changed the title
davidwalter2
left a comment
davidwalter2
left a comment
Collaborator
There was a problem hiding this comment.
A few more comments
| # what was applied without parsing the rabbit log. | ||
| if getattr(ifitter, "param_prior_active", False): | ||
| pm = ifitter.param_model | ||
| sigmas = np.asarray(pm.prior_sigmas, dtype=np.float64) |
Collaborator
There was a problem hiding this comment.
Should this not rather be tf objects instead of np? And use ifitter.indata.dtype here and below?
| x_nom + dxdtheta0[:, i] * shift, the Gaussian-approximation new | ||
| minimum for the shifted theta0. Should drastically reduce the | ||
| number of optimizer iterations on near-Gaussian nuisances. | ||
| x_nom + dxds[:, source] * shift, the Gaussian-approximation new |
Collaborator
There was a problem hiding this comment.
should it not be dxdx0 rather than dxds?
| t_lws = time.perf_counter() | ||
| dxdtheta0_tf, _, _ = fitter._dxdvars() | ||
| dxdtheta0_np = dxdtheta0_tf.numpy() | ||
| dxdx0_tf, _, _ = fitter._dxdvars() |
Collaborator
There was a problem hiding this comment.
just call it dxdx0 since tf is the default used elsewhere
| bin_by_bin_stat_mode, | ||
| systgroupidxs, | ||
| impacts_theta0_sq, | ||
| impacts_sources_sq, |
Collaborator
There was a problem hiding this comment.
should that not be impacts_x0_sq?
| else: | ||
| self.x.assign( | ||
| tf.concat([self.param_model.xparamdefault, self.theta0], axis=0) | ||
| tf.concat( |
Collaborator
There was a problem hiding this comment.
Is there now a self.x0default attribute defined above that let's you get rid of this if-else block?
| ], | ||
| axis=0, | ||
| ) | ||
| defaults = tf.concat( |
Collaborator
There was a problem hiding this comment.
Same comment as above, can self.x0default be used instead?
| ], | ||
| axis=0, | ||
| ) | ||
| self.var_x0_sources = tf.gather(self.var_x0, self.x0_source_idxs) |
Collaborator
There was a problem hiding this comment.
I think we should get rid of this attribute and use self.var_x0 everywhere instead, even if we carry some zeros with us, they should not matter
| return pd2ldbeta2 | ||
|
|
||
| def _dxdvars(self): | ||
| # Response of the postfit minimum to a unit shift of each constraint |
Collaborator
There was a problem hiding this comment.
Instead of working with this reduced set of source attributes I suggest to replace the self.theta0 from the previous implementation with self.x0 even if we have parameters without priors that we carry for no reason, I think it is better to have code simplicity versus overly specific logic
Per davidwalter2's latest comments, simplify to a single full-length x0 instead of carrying a reduced 'source' set: - _dxdvars / _dndvars differentiate w.r.t. the whole self.x0 (no scatter / perturbation, no x0= loss argument); columns for unconstrained centers (cw = 0) are exactly zero and carried along. - the global impacts work with the full impacts_x0 directly: impacts_x0_sq = square(impacts_x0), var_x0 = reduce_sum, no gather; the per-source axis is now the whole parameter list (matching the traditional impacts). gaussian per-source = dxdx0 * sqrt(var_x0) (one expression). - grouping is one operation over full-x column-index lists (syst groups shifted by nmodel_params, ParamModel impact groups in the param block) so param and syst sources combine uniformly; the grouped global axis now uses the param-impact groups (like the traditional grouped impacts). - drop x0_source_idxs / var_x0_sources; param_prior_idxs kept only for the asym scan. xdefaultassign / bayesassign use self.x0default. - nits: drop the param-model-specific prior comment; meta priors use indata.dtype; rename dxds -> dxdx0 in global_asym_impacts. Validated: lint clean, pytest (15) + the five CI script-tests pass, and the ParamModel-prior impacts are unchanged (likelihood = gaussian source[sig] = 1.105769e-02, asym (-1.097e-2, +1.114e-2)). NOTE for review: the global-impact source axis is now the full parameter list (was systs + <param>_prior); priored params appear under their own name with unconstrained params as exact-zero columns. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Based on #140 #143 #144 to be merged after those.