From fbd4228517b1a44620f3a8d263844d7d73cec532 Mon Sep 17 00:00:00 2001 From: Luca Lavezzo Date: Thu, 30 Apr 2026 09:14:14 -0400 Subject: [PATCH 1/7] asymmetric traditoinal impacts; unblinding groups --- bin/rabbit_fit.py | 114 ++++++++++++ rabbit/fitter.py | 306 +++++++++++++++++++++++++++------ rabbit/impacts/asym_impacts.py | 168 ++++++++++++++++++ rabbit/io_tools.py | 8 +- rabbit/parsing.py | 14 ++ 5 files changed, 558 insertions(+), 52 deletions(-) create mode 100644 rabbit/impacts/asym_impacts.py diff --git a/bin/rabbit_fit.py b/bin/rabbit_fit.py index dcf8b1a..9eb253e 100755 --- a/bin/rabbit_fit.py +++ b/bin/rabbit_fit.py @@ -212,6 +212,56 @@ def make_parser(): action="store_true", help="compute impacts of frozen (non-profiled) systematics", ) + parser.add_argument( + "--asymImpacts", + default=False, + action="store_true", + help="Compute traditional asymmetric impacts on POIs by running a " + "Delta(2NLL)=1 contour scan per nuisance. Skips structurally " + "symmetric and unconstrained nuisances by default.", + ) + parser.add_argument( + "--asymImpactsInclude", + default=None, + nargs="+", + help="Regex(es) restricting which nuisances are scanned for --asymImpacts.", + ) + parser.add_argument( + "--asymImpactsExclude", + default=None, + nargs="+", + help="Regex(es) excluding nuisances from --asymImpacts.", + ) + parser.add_argument( + "--asymImpactsAll", + default=False, + action="store_true", + help="Disable the structural-symmetric skip in --asymImpacts and scan " + "every constrained nuisance regardless of template content.", + ) + parser.add_argument( + "--asymImpactsHess", + default="exact", + choices=["exact", "hvp", "frozen", "bfgs", "sr1"], + help="Constraint-Hessian mode for the contour-scan in --asymImpacts. " + "'exact' rebuilds the full N x N NLL Hessian each iteration (slow, " + "reference). 'hvp' uses Hessian-vector products via a LinearOperator " + "(exact, typically fastest for large N). 'frozen' uses cov^-1 from " + "the postfit as a constant Hessian (cheapest but produces silent " + "failures on non-Gaussian profiles -- speed reference only). " + "'bfgs'/'sr1' use a quasi-Newton estimate built up from the gradient " + "sequence (no extra TF calls per iteration; may need more iterations).", + ) + parser.add_argument( + "--asymImpactsTol", + default=1e-6, + type=float, + help="trust-constr xtol/gtol for the --asymImpacts contour-scan. " + "Looser values (e.g. 1e-3) terminate before the iterate is on the " + "Delta(2NLL)=q contour, producing silent constraint violations. " + "Tighter values (1e-5, 1e-6) are slower with no benefit unless your " + "fit has nuisances whose profile is far from quadratic." + ) parser.add_argument( "--lCurveScan", default=False, @@ -328,6 +378,37 @@ def save_hists(args, mappings, fitter, ws, prefit=True, profile=False): logger.info(f" 2*deltaNLL: {round(chi2val, 2)}") logger.info(rf" p-value: {round(p_val * 100, 2)}%") + # DEBUG: decompose chi2_sat = 2*(NLL_nom - NLL_sat) + # = 2*(ln_nom-ln_sat) [data-vs-prediction misfit] + # + 2*(lc_nom-lc_sat) [constrained-nuisance pull cost] + # + 2*(lbeta_nom-lbeta_sat) [BBB pull cost] + ln_sat_t, lc_sat_t, lbeta_sat_t, _, _ = ( + fitter_saturated._compute_nll_components(profile=True) + ) + ln_sat = float(ln_sat_t.numpy()) + lc_sat = float(lc_sat_t.numpy()) + lbeta_sat = ( + float(lbeta_sat_t.numpy()) if lbeta_sat_t is not None else 0.0 + ) + d_ln = 2.0 * (ws.results["ln_nom"] - ln_sat) + d_lc = 2.0 * (ws.results["lc_nom"] - lc_sat) + d_lbeta = 2.0 * (ws.results["lbeta_nom"] - lbeta_sat) + d_sum = d_ln + d_lc + d_lbeta + logger.info(f"Saturated chi2 decomposition for {mapping.key}:") + logger.info(f" 2*Delta ln (data-vs-pred) : {d_ln:.4f}") + logger.info(f" 2*Delta lc (nuis. pulls) : {d_lc:.4f}") + logger.info(f" 2*Delta lbeta (BBB pulls) : {d_lbeta:.4f}") + logger.info( + f" sum : {d_sum:.4f}" + f" (closure vs chi2val={chi2val:.4f}: {d_sum - chi2val:+.4f})" + ) + if chi2val > 0: + logger.info( + f" fractions: data {d_ln / chi2val * 100:5.1f}% " + f"nuis {d_lc / chi2val * 100:5.1f}% " + f"BBB {d_lbeta / chi2val * 100:5.1f}%" + ) + ws.add_chi2(chi2val, ndf, prefit, mapping, saturated=True) if args.saveHistsPerProcess and not mapping.skip_per_process: @@ -448,6 +529,22 @@ def fit(args, fitter, ws, dofit=True): nllvalreduced = fitter.reduced_nll().numpy() + # DEBUG: decompose nominal-fit NLL into (Poisson data, nuisance constraints, BBB) + # so the per-projection saturated chi2 below can be split into: + # 2*Delta ln -> data-vs-prediction misfit + # 2*Delta lc -> constrained-nuisance pull cost + # 2*Delta lb -> BBB pull cost + ln_nom_t, lc_nom_t, lbeta_nom_t, _, _ = fitter._compute_nll_components(profile=True) + ln_nom = float(ln_nom_t.numpy()) + lc_nom = float(lc_nom_t.numpy()) + lbeta_nom = float(lbeta_nom_t.numpy()) if lbeta_nom_t is not None else 0.0 + logger.info("NLL decomposition (nominal postfit):") + logger.info(f" ln (Poisson data) : {ln_nom:.4f}") + logger.info(f" lc (nuis. constr.): {lc_nom:.4f}") + logger.info(f" lbeta (BBB) : {lbeta_nom:.4f}") + logger.info(f" sum : {ln_nom + lc_nom + lbeta_nom:.4f}") + logger.info(f" nllvalreduced : {nllvalreduced:.4f}") + ndfsat = ( tf.size(fitter.nobs) - fitter.poi_model.npoi - fitter.indata.nsystnoconstraint ).numpy() @@ -463,6 +560,9 @@ def fit(args, fitter, ws, dofit=True): ws.results.update( { "nllvalreduced": nllvalreduced, + "ln_nom": ln_nom, + "lc_nom": lc_nom, + "lbeta_nom": lbeta_nom, "ndfsat": ndfsat, "edmval": edmval, "postfit_profile": not args.noPostfitProfileBB, @@ -486,6 +586,20 @@ def fit(args, fitter, ws, dofit=True): *fitter.nonprofiled_impacts_parms(), base_name="nonprofiled_impacts_asym" ) + if args.asymImpacts: + ws.add_impacts_asym_hist( + *fitter.asym_impacts_parms( + nll_min=fitter.reduced_nll().numpy(), + include=args.asymImpactsInclude, + exclude=args.asymImpactsExclude, + skip_symmetric=not args.asymImpactsAll, + hess_mode=args.asymImpactsHess, + contour_xtol=args.asymImpactsTol, + contour_gtol=args.asymImpactsTol, + ), + base_name="impacts_asym", + ) + # Likelihood scans if args.scan is not None: parms = np.array(fitter.parms).astype(str) if len(args.scan) == 0 else args.scan diff --git a/rabbit/fitter.py b/rabbit/fitter.py index 5e3c29f..68304dc 100644 --- a/rabbit/fitter.py +++ b/rabbit/fitter.py @@ -11,7 +11,12 @@ from rabbit import io_tools from rabbit import tfhelpers as tfh -from rabbit.impacts import global_impacts, nonprofiled_impacts, traditional_impacts +from rabbit.impacts import ( + asym_impacts, + global_impacts, + nonprofiled_impacts, + traditional_impacts, +) from rabbit.tfhelpers import edmval_cov logger = logging.child_logger(__name__) @@ -140,6 +145,7 @@ def __init__( poi_model, options.setConstraintMinimum, unblind=options.unblind, + blinding_group=getattr(options, "blindingGroup", []), freeze_parameters=options.freezeParameters, ) @@ -261,6 +267,7 @@ def init_fit_parms( poi_model, set_constraint_minimum=[], unblind=False, + blinding_group=[], freeze_parameters=[], ): self.poi_model = poi_model @@ -276,7 +283,7 @@ def init_fit_parms( trainable=False, name="offset_theta", ) - self.init_blinding_values(unblind) + self.init_blinding_values(unblind, blinding_group) self.parms = np.concatenate([self.poi_model.pois, self.indata.systs]) @@ -441,14 +448,16 @@ def defreeze_params(self, unfrozen_parmeter_expressions): ] self.update_frozen_params() - def init_blinding_values(self, unblind_parameter_expressions=[]): + def init_blinding_values( + self, unblind_parameter_expressions=[], blinding_group_expressions=[] + ): logger.debug(f"Unblind parameters with {unblind_parameter_expressions}") + all_param_names = [ + *self.poi_model.pois, + *[self.indata.systs[i] for i in self.indata.noiidxs], + ] unblind_parameters = match_regexp_params( - unblind_parameter_expressions, - [ - *self.poi_model.pois, - *[self.indata.systs[i] for i in self.indata.noiidxs], - ], + unblind_parameter_expressions, all_param_names ) # check if dataset is an integer (i.e. if it is real data or not) and use this to choose the random seed @@ -473,14 +482,47 @@ def deterministic_random_from_string(s, mean=0.0, std=5.0): value = rng.normal(loc=mean, scale=std) return value + # Build a map: param name -> seed string. Default is the param's own name; for + # parameters matching a blinding group the seed is the group regex string, so all + # members of the group share an identical deterministic offset (preserving relative + # differences while blinding absolute values). + if isinstance(blinding_group_expressions, str): + blinding_group_expressions = [blinding_group_expressions] + param_to_seed = {} + param_to_group = {} + for expr in blinding_group_expressions: + matched = match_regexp_params(expr, all_param_names) + for p in matched: + if p in param_to_seed: + continue # first matching group wins + param_to_seed[p] = expr + param_to_group[p] = expr + + # Refuse to run if --unblind and --blindingGroup match the same parameter: + # the user's intent is ambiguous (unblind it, or share a group offset?), and + # silently picking one risks accidentally unblinding a sensitive parameter. + overlap = [p for p in unblind_parameters if p in param_to_group] + if overlap: + details = ", ".join( + f"{p.decode() if isinstance(p, bytes) else p}" + f" (group '{param_to_group[p]}')" + for p in overlap + ) + raise RuntimeError( + "The following parameters match both --unblind and --blindingGroup; " + "refusing to proceed to avoid an ambiguous (un)blinding. " + f"Tighten the regexes to make the intent explicit: {details}" + ) + # multiply offset to nois self._blinding_values_theta = np.zeros(self.indata.nsyst, dtype=np.float64) for i in self.indata.noiidxs: param = self.indata.systs[i] if param in unblind_parameters: continue - logger.debug(f"Blind parameter {param}") - value = deterministic_random_from_string(param) + seed = param_to_seed.get(param, param) + logger.debug(f"Blind parameter {param} (seed='{seed}')") + value = deterministic_random_from_string(seed) self._blinding_values_theta[i] = value # add offset to pois @@ -489,8 +531,9 @@ def deterministic_random_from_string(s, mean=0.0, std=5.0): param = self.poi_model.pois[i] if param in unblind_parameters: continue - logger.debug(f"Blind signal strength modifier for {param}") - value = deterministic_random_from_string(param) + seed = param_to_seed.get(param, param) + logger.debug(f"Blind signal strength modifier for {param} (seed='{seed}')") + value = deterministic_random_from_string(seed) self._blinding_values_poi[i] = np.exp(value) def set_blinding_offsets(self, blind=True): @@ -903,6 +946,78 @@ def gaussian_global_impacts_parms(self): return impacts, impacts_grouped + def asymmetric_nuisance_mask(self, atol=0.0): + """Boolean mask of length nsyst, True where the nuisance has nonzero + asymmetric (logkhalfdiff) tensor content.""" + return asym_impacts.asymmetric_nuisance_mask(self.indata, atol=atol) + + def asym_impacts_parms( + self, + nll_min=None, + q=1, + include=None, + exclude=None, + skip_symmetric=True, + skip_unconstrained=True, + contour_xtol=1e-6, + contour_gtol=1e-6, + contour_maxiter=200, + hess_mode="exact", + ): + """Traditional asymmetric impacts via per-nuisance contour scan. + + Args: + nll_min: postfit reduced NLL. Computed from the current fit state + if None. + q: contour level (q=1 -> 1 sigma). + include: optional regex(es) restricting which nuisances to scan. + exclude: optional regex(es) excluding nuisances from the scan. + skip_symmetric: skip nuisances whose template content is structurally + symmetric (their asymmetric impact equals the Gaussian impact). + skip_unconstrained: skip nuisances with constraintweight=0 (no + finite Delta(2NLL)=q contour). + """ + if nll_min is None: + nll_min = float(self.reduced_nll().numpy()) + + nsyst = self.indata.nsyst + cw = self.indata.constraintweights.numpy() + syst_names = np.array(self.indata.systs).astype(bytes) + + selected = np.ones(nsyst, dtype=bool) + if skip_unconstrained: + selected &= cw > 0 + if skip_symmetric: + selected &= self.asymmetric_nuisance_mask() + if include is not None: + keep = match_regexp_params(include, syst_names) + keep_set = set(keep) + selected &= np.array([n in keep_set for n in syst_names]) + if exclude is not None: + drop = match_regexp_params(exclude, syst_names) + drop_set = set(drop) + selected &= np.array([n not in drop_set for n in syst_names]) + + selected_idxs = np.where(selected)[0] + selected_names = syst_names[selected_idxs] + + logger.info( + f"asym_impacts_parms: selected {len(selected_idxs)}/{nsyst} nuisances " + f"(skip_symmetric={skip_symmetric}, skip_unconstrained={skip_unconstrained})" + ) + + return asym_impacts.asym_impacts_parms( + self, + nll_min, + selected_idxs, + selected_names, + q=q, + contour_xtol=contour_xtol, + contour_gtol=contour_gtol, + contour_maxiter=contour_maxiter, + hess_mode=hess_mode, + ) + def nonprofiled_impacts_parms(self, unconstrained_err=1.0): return nonprofiled_impacts.nonprofiled_impacts_parms( self.x, @@ -2395,23 +2510,111 @@ def nll_scan2D(self, param_tuple, scan_range, scan_points, use_prefit=False): return x_scans, y_scans, dnlls - def contour_scan(self, param, nll_min, q=1, signs=[-1, 1], fun=None): - # TODO this is basically traditional asymmetric impacts - def scipy_loss(x): - self.x.assign(x) - val = self.loss_val() - loss = val.numpy() - nll_min - 0.5 * q - return loss[None,] - - def scipy_grad(x): + def contour_scan( + self, + param, + nll_min, + q=1, + signs=[-1, 1], + fun=None, + xtol=1e-6, + gtol=1e-6, + maxiter=200, + hess_mode="exact", + ): + # Layered cache: trust-constr calls scipy_loss many times during line + # search (only val needed), and scipy_grad / scipy_hess on accepted + # steps. The cache is keyed by x content so repeated requests at the + # same point are free. + lg_cache = {"x": None, "val": None, "grad": None} + + def _ensure_loss_grad(x): + if lg_cache["x"] is not None and np.array_equal(lg_cache["x"], x): + return self.x.assign(x) val, grad = self.loss_val_grad() - return grad.numpy()[None,] + lg_cache["x"] = np.array(x, copy=True) + lg_cache["val"] = float(val.numpy()) - nll_min - 0.5 * q + lg_cache["grad"] = grad.numpy() + + # Constraint Hessian. Modes: + # "exact": recompute the full NLL Hessian at every accepted iteration + # (~25 s/eval for thousands of params; dominant cost). Reference. + # "hvp": LinearOperator whose matvec computes one Hessian-vector + # product via a nested GradientTape (~2x the cost of a gradient). + # Exact (no approximation); avoids materializing the N x N matrix. + # trust-constr only multiplies H against trial directions in its + # inner CG, so HVP is typically much faster than "exact". + # "frozen": constant Hessian = postfit precision matrix (cov^-1), + # computed once. Cheapest, but the Lagrangian model is wrong off + # the postfit, so trust-constr's KKT/optimality criterion can be + # satisfied while the constraint violation is large -- producing + # silent failures on non-Gaussian profiles. Useful only as a + # speed reference, not as a production default. + # "bfgs" / "sr1": quasi-Newton Hessian estimate built up by + # trust-constr from the gradient sequence (no extra TF calls). + # Cheapest per iteration; may need more iterations to converge. + # SR1 is more robust for non-convex local geometry than BFGS. + if hess_mode == "hvp": + def scipy_hess(x, v): + _ensure_loss_grad(x) + n = len(x) + scale = float(v[0]) + + def _matvec(p): + self.x.assign(x) + p_tf = tf.convert_to_tensor(p, dtype=self.indata.dtype) + _, _, hp = self.loss_val_grad_hessp(p_tf) + return scale * hp.numpy() - def scipy_hess(x, v): - self.x.assign(x) - val, grad, hess = self.loss_val_grad_hess() - return v[0] * hess.numpy() + return scipy.sparse.linalg.LinearOperator( + shape=(n, n), matvec=_matvec, dtype=np.float64 + ) + + elif hess_mode == "frozen": + postfit_hess = np.linalg.inv(self.cov.numpy()) + + def scipy_hess(x, v): + return v[0] * postfit_hess + + elif hess_mode == "bfgs": + scipy_hess = scipy.optimize.BFGS() + + elif hess_mode == "sr1": + scipy_hess = scipy.optimize.SR1() + + elif hess_mode == "exact": + h_cache = {"x": None, "hess": None} + + def _ensure_hess(x): + if h_cache["x"] is not None and np.array_equal(h_cache["x"], x): + return + self.x.assign(x) + val, grad, hess = self.loss_val_grad_hess() + h_cache["x"] = np.array(x, copy=True) + h_cache["hess"] = hess.numpy() + # opportunistically refresh the loss/grad cache. + lg_cache["x"] = h_cache["x"] + lg_cache["val"] = float(val.numpy()) - nll_min - 0.5 * q + lg_cache["grad"] = grad.numpy() + + def scipy_hess(x, v): + _ensure_hess(x) + return v[0] * h_cache["hess"] + + else: + raise ValueError( + f"contour_scan: unknown hess_mode={hess_mode!r}; " + "expected one of 'exact', 'hvp', 'frozen', 'bfgs', 'sr1'." + ) + + def scipy_loss(x): + _ensure_loss_grad(x) + return np.array([lg_cache["val"]]) + + def scipy_grad(x): + _ensure_loss_grad(x) + return lg_cache["grad"][None, :] nlc = scipy.optimize.NonlinearConstraint( fun=scipy_loss, @@ -2425,30 +2628,24 @@ def scipy_hess(x, v): params_values = np.full((len(signs), len(self.parms)), np.nan) xval = tf.identity(self.x) - xval_init = xval.numpy() + xval_np = xval.numpy() idx = np.where(self.parms.astype(str) == param)[0][0] - x0 = xval[idx] - - # initial guess from covariance - initial_fit = False - - xup = xval[idx] + (self.cov[idx, idx] * q) ** 0.5 - xdn = xval[idx] - (self.cov[idx, idx] * q) ** 0.5 + x0 = xval_np[idx] + + # Gaussian-optimal warm start on the Delta(2NLL)=q contour: + # maximizing +/- dx[idx] subject to dx^T H dx = q (with H = cov^{-1}) + # gives dx = sign * sqrt(q) * cov[:,idx] / sqrt(cov[idx,idx]). + # This lands all parameters on the contour in the Gaussian limit and + # is exact for nuisances with a near-quadratic likelihood, so the + # constrained minimization typically converges in just a few steps. + cov_col = self.cov[:, idx].numpy() + sigma_idx = float(self.cov[idx, idx].numpy()) ** 0.5 + gauss_dx = (q**0.5) * cov_col / sigma_idx for i, sign in enumerate(signs): - # Objective function and its derivatives - if sign == -1: - xval_init[idx] = xdn - else: - xval_init[idx] = xup - - if initial_fit: - # perform initial fit where contour is expected - self.x.assign(xval_init) - self.freeze_params(param) - self.minimize() - self.defreeze_params(param) + xval_init = xval_np + sign * gauss_dx + t_side0 = time.perf_counter() opt = {} if fun is None: @@ -2511,20 +2708,27 @@ def objective_hessp(x, pval): jac=True, constraints=[nlc], options={ - "maxiter": 50000, - "xtol": 1e-10, - "gtol": 1e-10, - # "barrier_tol": 1e-10, + "maxiter": maxiter, + "xtol": xtol, + "gtol": gtol, }, **opt, ) - logger.info(f"Success: {res.success}") + t_side = time.perf_counter() - t_side0 + logger.info( + f"Success: {res.success} sign={sign} time={t_side:.2f}s " + f"(nit={getattr(res, 'nit', '?')}, " + f"nfev={getattr(res, 'nfev', '?')}, " + f"njev={getattr(res, 'njev', '?')}, " + f"nhev={getattr(res, 'nhev', '?')})" + ) logger.debug(f"Status: {res.status}") if not res.success: logger.warning(f"Message: {res.message}") logger.warning(f"Optimality (gtol): {res.optimality}") logger.warning(f"Constraint Violation: {res.constr_violation}") + self.x.assign(xval) continue params_values[i] = res["x"] - xval diff --git a/rabbit/impacts/asym_impacts.py b/rabbit/impacts/asym_impacts.py new file mode 100644 index 0000000..3b76929 --- /dev/null +++ b/rabbit/impacts/asym_impacts.py @@ -0,0 +1,168 @@ +""" +Traditional asymmetric impacts. + +For each selected nuisance, find the asymmetric +/- 1 sigma points on the +Delta(2NLL)=q likelihood contour via constrained minimization (contour_scan). +The shifts of every fitter parameter at those points are the asymmetric +impacts. Group impacts are obtained by quadrature envelope of the contained +nuisances, separately for the down and up sides. + +Nuisances that are structurally symmetric (logkhalfdiff identically zero) +or unconstrained (constraintweight = 0) are skipped by default: the first +case has asymmetric impact equal to the Gaussian impact (already produced by +traditional_impacts.impacts_parms), and the second has no finite Delta(2NLL) +contour. +""" + +import time + +import numpy as np +from wums import logging + +logger = logging.child_logger(__name__) + + +def asymmetric_nuisance_mask(indata, atol=0.0): + """Boolean mask of length nsyst, True for nuisances with nonzero asymmetric + (logkhalfdiff) tensor content. False if the entire tensor is symmetric or + if a particular nuisance has zero halfdiff.""" + if indata.symmetric_tensor: + return np.zeros(indata.nsyst, dtype=bool) + + nsyst = indata.nsyst + if indata.sparse: + idx = indata.logk.indices.numpy() + vals = indata.logk.values.numpy() + syst_col = idx[:, -1] + halfdiff_entries = (syst_col >= nsyst) & (np.abs(vals) > atol) + nz_systs = np.unique(syst_col[halfdiff_entries]) - nsyst + mask = np.zeros(nsyst, dtype=bool) + mask[nz_systs] = True + return mask + + # Dense layout: logk has shape [nbinsfull, nproc, 2, nsyst] when asymmetric; + # axis -2 is [logkavg, logkhalfdiff]. Reduce over (bin, proc) for halfdiff. + halfdiff = indata.logk[..., 1, :].numpy() + axes = tuple(range(halfdiff.ndim - 1)) + return np.any(np.abs(halfdiff) > atol, axis=axes) + + +def _envelope(values): + """Quadrature envelope of asymmetric impacts within a group. + + Args: + values: numpy array of shape (n_in_group, 2, n_total_params), where + axis 1 is [down, up]. + + Returns: + Array of shape (2, n_total_params) with the envelope's lower (negative) + and upper (positive) sides. + """ + zeros = np.zeros((values.shape[0], values.shape[-1]), dtype=values.dtype) + vmin = np.min(values, axis=1) + vmax = np.max(values, axis=1) + lower = -np.sqrt(np.sum(np.minimum(zeros, vmin) ** 2, axis=0)) + upper = np.sqrt(np.sum(np.maximum(zeros, vmax) ** 2, axis=0)) + return np.stack([lower, upper]) + + +def asym_impacts_parms( + fitter, + nll_min, + selected_idxs, + selected_names, + q=1, + contour_xtol=1e-4, + contour_gtol=1e-4, + contour_maxiter=200, + hess_mode="exact", +): + """Run a per-nuisance contour scan and assemble the asymmetric-impact tensor. + + Args: + fitter: the Fitter instance (used for contour_scan and indata). + nll_min: postfit reduced NLL. + selected_idxs: indices into the syst axis (0..nsyst-1) of nuisances to scan. + selected_names: names of those nuisances (bytes), used as impact-axis labels. + q: contour level (q=1 -> 1 sigma, q=4 -> 2 sigma). + + Returns: + parms: np.ndarray of bytes, the impact-axis labels. + impacts: np.ndarray of shape (n_scanned, 2, n_total_params). + Axis 1 is [down, up] matching axis_downUpVar. + group_names: np.ndarray of bytes for groups containing scanned nuisances + (plus a trailing "Total"). + impacts_grouped: np.ndarray of shape (n_groups, 2, n_total_params). + """ + n_scanned = len(selected_idxs) + n_total = len(fitter.parms) + impacts = np.zeros((n_scanned, 2, n_total)) + + logger.info(f"asym_impacts: scanning {n_scanned} nuisances") + + t_per = np.zeros(n_scanned, dtype=np.float64) + t_total0 = time.perf_counter() + + for i, name in enumerate(selected_names): + name_str = name.decode() if isinstance(name, bytes) else name + logger.info(f" [{i + 1}/{n_scanned}] contour scan for {name_str}") + t0 = time.perf_counter() + _, params_values = fitter.contour_scan( + name_str, + nll_min, + q=q, + signs=[-1, 1], + xtol=contour_xtol, + gtol=contour_gtol, + maxiter=contour_maxiter, + hess_mode=hess_mode, + ) + t_per[i] = time.perf_counter() - t0 + logger.info(f" took {t_per[i]:.2f}s") + # signs=[-1, +1] -> bin 0 = down, bin 1 = up (matches axis_downUpVar). + # params_values rows are NaN where convergence failed; leave as zero impact. + if not np.any(np.isnan(params_values)): + impacts[i] = params_values + else: + valid = ~np.any(np.isnan(params_values), axis=1) + impacts[i, valid] = params_values[valid] + + if n_scanned > 0: + t_total = time.perf_counter() - t_total0 + logger.info( + f"asym_impacts: total {t_total:.1f}s " + f"(mean {t_per.mean():.2f}s, min {t_per.min():.2f}s, " + f"max {t_per.max():.2f}s per nuisance)" + ) + + # Build grouped impacts via quadrature envelope, separately for down/up. + syst_names = np.array(fitter.indata.systs).astype(bytes) + selected_set = set(selected_idxs.tolist()) + pos_in_scanned = {int(idx): k for k, idx in enumerate(selected_idxs)} + + group_names = [] + group_impacts = [] + for gname, gidxs in zip(fitter.indata.systgroups, fitter.indata.systgroupidxs): + gidxs = np.asarray(gidxs).astype(int) + in_scanned = [pos_in_scanned[i] for i in gidxs if int(i) in selected_set] + if not in_scanned: + continue + sub = impacts[in_scanned] + group_names.append(gname) + group_impacts.append(_envelope(sub)) + + if n_scanned > 0: + group_names.append(b"Total") + group_impacts.append(_envelope(impacts)) + + if group_impacts: + impacts_grouped = np.stack(group_impacts) + else: + impacts_grouped = np.zeros((0, 2, n_total)) + + return ( + np.asarray(selected_names), + impacts, + np.asarray(group_names), + impacts_grouped, + ) diff --git a/rabbit/io_tools.py b/rabbit/io_tools.py index ad3f89a..fabbd2c 100644 --- a/rabbit/io_tools.py +++ b/rabbit/io_tools.py @@ -47,7 +47,13 @@ def read_impacts_poi( ): # read impacts of a single POI - if asym and impact_type == "traditional": + if ( + asym + and impact_type == "traditional" + and "impacts_asym" not in fitresult.keys() + ): + # Fallback: read asymmetric traditional impacts from a generic + # contour scan output if --asymImpacts wasn't run. h_impacts = fitresult["contour_scans"].get()[{"confidence_level": "1.0"}] else: impact_name = "impacts" diff --git a/rabbit/parsing.py b/rabbit/parsing.py index 1945ece..3bc9043 100644 --- a/rabbit/parsing.py +++ b/rabbit/parsing.py @@ -237,6 +237,20 @@ def common_parser(): E.g. use '--unblind ^signal$' to unblind a parameter named signal or '--unblind' to unblind all. """, ) + parser.add_argument( + "--blindingGroup", + type=str, + default=[], + nargs="*", + help=""" + Specify list of regex defining groups of parameters that share a single deterministic + blinding offset (seeded from the regex string itself). Useful to keep relative pulls / + differences between matched parameters meaningful while still blinding their absolute + values. E.g. '--blindingGroup ^alphaS_y\\d+$' applies the same offset to all alphaS + rapidity-bin parameters. Parameters not matching any group keep per-name blinding. + Overlap with --unblind is treated as a configuration error and aborts the fit. + """, + ) parser.add_argument( "--setConstraintMinimum", default=[], From b49d09e4bf8a2b0c2bb8282508d080b669e10969 Mon Sep 17 00:00:00 2001 From: Luca Lavezzo <38077849+lucalavezzo@users.noreply.github.com> Date: Thu, 30 Apr 2026 09:42:16 -0400 Subject: [PATCH 2/7] fix --- rabbit/fitter.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rabbit/fitter.py b/rabbit/fitter.py index 8a03f15..86f04e4 100644 --- a/rabbit/fitter.py +++ b/rabbit/fitter.py @@ -538,7 +538,7 @@ def init_blinding_values( ): logger.debug(f"Unblind parameters with {unblind_parameter_expressions}") all_param_names = [ - *self.param_model.params[: self.param_model.npoi, + *self.param_model.params[: self.param_model.npoi], *[self.indata.systs[i] for i in self.indata.noiidxs], ] unblind_parameters = match_regexp_params( From 467a7454a2cd4ebdd36b4cf6033c167829354353 Mon Sep 17 00:00:00 2001 From: Luca Lavezzo Date: Mon, 4 May 2026 09:35:07 -0400 Subject: [PATCH 3/7] global asym impacts --- bin/rabbit_fit.py | 65 ++++- rabbit/fitter.py | 88 +++++-- rabbit/impacts/global_asym_impacts.py | 209 ++++++++++++++++ tests/test_global_asym_impacts.py | 338 ++++++++++++++++++++++++++ 4 files changed, 685 insertions(+), 15 deletions(-) create mode 100644 rabbit/impacts/global_asym_impacts.py create mode 100644 tests/test_global_asym_impacts.py diff --git a/bin/rabbit_fit.py b/bin/rabbit_fit.py index f535244..e30c145 100755 --- a/bin/rabbit_fit.py +++ b/bin/rabbit_fit.py @@ -266,6 +266,15 @@ def make_parser(): "'bfgs'/'sr1' use a quasi-Newton estimate built up from the gradient " "sequence (no extra TF calls per iteration; may need more iterations).", ) + parser.add_argument( + "--asymImpactsMaxiter", + default=200, + type=int, + help="trust-constr maxiter for the --asymImpacts contour-scan. Lower " + "values cap the cost of nuisances that don't converge (useful for " + "sanity benchmarks); higher values give genuinely-non-Gaussian " + "nuisances more chances to converge.", + ) parser.add_argument( "--asymImpactsTol", default=1e-6, @@ -274,7 +283,49 @@ def make_parser(): "Looser values (e.g. 1e-3) terminate before the iterate is on the " "Delta(2NLL)=q contour, producing silent constraint violations. " "Tighter values (1e-5, 1e-6) are slower with no benefit unless your " - "fit has nuisances whose profile is far from quadratic." + "fit has nuisances whose profile is far from quadratic.", + ) + parser.add_argument( + "--globalAsymImpacts", + default=False, + action="store_true", + help="Compute fully likelihood-based asymmetric global impacts on POIs " + "by shifting each constrained nuisance's theta0 by +/- 1 prefit sigma " + "and re-running the fit. In the Gaussian limit this reproduces " + "--gaussianGlobalImpacts; deviations measure non-Gaussianity of the " + "joint profile. Cost is N_selected x 2 full minimizations -- gate with " + "--globalAsymImpactsInclude in practice.", + ) + parser.add_argument( + "--globalAsymImpactsInclude", + default=None, + nargs="+", + help="Regex(es) restricting which nuisances are scanned for " + "--globalAsymImpacts.", + ) + parser.add_argument( + "--globalAsymImpactsExclude", + default=None, + nargs="+", + help="Regex(es) excluding nuisances from --globalAsymImpacts.", + ) + parser.add_argument( + "--globalAsymImpactsSigma", + default=1.0, + type=float, + help="theta0 shift magnitude for --globalAsymImpacts, in units of the " + "prefit constraint width (1.0 = 1 prefit sigma).", + ) + parser.add_argument( + "--globalAsymImpactsLinearWarmstart", + default=False, + action="store_true", + help="EXPERIMENTAL: warm-start each --globalAsymImpacts refit at the " + "Gaussian-approximation new minimum x_nom + dxdtheta0[:, i] * shift " + "(same Jacobian as --gaussianGlobalImpacts). On near-Gaussian " + "nuisances this should reduce per-nuisance refit cost by 10-50x. " + "Adds one --gaussianGlobalImpacts-equivalent precompute up front. " + "Off by default until validated on real tensors.", ) parser.add_argument( "--lCurveScan", @@ -665,10 +716,22 @@ def fit(args, fitter, ws, dofit=True): hess_mode=args.asymImpactsHess, contour_xtol=args.asymImpactsTol, contour_gtol=args.asymImpactsTol, + contour_maxiter=args.asymImpactsMaxiter, ), base_name="impacts_asym", ) + if args.globalAsymImpacts: + ws.add_impacts_asym_hist( + *fitter.global_asym_impacts_parms( + include=args.globalAsymImpactsInclude, + exclude=args.globalAsymImpactsExclude, + sigma=args.globalAsymImpactsSigma, + linear_warmstart=args.globalAsymImpactsLinearWarmstart, + ), + base_name="global_impacts_asym", + ) + # Likelihood scans if args.scan is not None: parms = np.array(fitter.parms).astype(str) if len(args.scan) == 0 else args.scan diff --git a/rabbit/fitter.py b/rabbit/fitter.py index 8a03f15..8c8bb75 100644 --- a/rabbit/fitter.py +++ b/rabbit/fitter.py @@ -14,6 +14,7 @@ from rabbit import tfhelpers as tfh from rabbit.impacts import ( asym_impacts, + global_asym_impacts, global_impacts, nonprofiled_impacts, traditional_impacts, @@ -28,23 +29,26 @@ def solve_quad_eq(a, b, c): def match_regexp_params(regular_expressions, parameter_names): + # Match parameter names against a list of expressions where each entry may + # be either an exact parameter name or a regex (re.match, anchored to the + # start). Returns the union of exact and regex matches, preserving the + # parameter_names order and de-duplicating. Mixing exact and regex entries + # in the same call is supported. if isinstance(regular_expressions, str): regular_expressions = [regular_expressions] - # Check for exact matches first - exact_matches = [ - s for expr in regular_expressions for s in parameter_names if s.decode() == expr - ] - if exact_matches: - return exact_matches - - # Fall back to regex matching + exact_lookup = set(regular_expressions) compiled_expressions = [re.compile(expr) for expr in regular_expressions] - return [ - s - for s in parameter_names - if any(regex.match(s.decode()) for regex in compiled_expressions) - ] + + matched = [] + seen = set() + for s in parameter_names: + decoded = s.decode() if hasattr(s, "decode") else s + if decoded in exact_lookup or any(r.match(decoded) for r in compiled_expressions): + if decoded not in seen: + seen.add(decoded) + matched.append(s) + return matched class FitterCallback: @@ -538,7 +542,7 @@ def init_blinding_values( ): logger.debug(f"Unblind parameters with {unblind_parameter_expressions}") all_param_names = [ - *self.param_model.params[: self.param_model.npoi, + *self.param_model.params[: self.param_model.npoi], *[self.indata.systs[i] for i in self.indata.noiidxs], ] unblind_parameters = match_regexp_params( @@ -1209,6 +1213,61 @@ def asym_impacts_parms( hess_mode=hess_mode, ) + def global_asym_impacts_parms( + self, + include=None, + exclude=None, + skip_unconstrained=True, + sigma=1.0, + linear_warmstart=False, + ): + """Fully likelihood-based asymmetric global impacts. + + For each selected nuisance i, shift theta0[i] by +/- sigma (in units of + the prefit constraint width) and re-run the full fit. POI shifts at + each sign are the asymmetric global impacts. + + Args: + include: optional regex(es) restricting which nuisances to scan. + exclude: optional regex(es) excluding nuisances from the scan. + skip_unconstrained: skip nuisances with constraintweight=0; their + "1 prefit sigma" is undefined. + sigma: shift magnitude in prefit-sigma units. + linear_warmstart: experimental, see + global_asym_impacts.global_asym_impacts_parms. + """ + nsyst = self.indata.nsyst + cw = self.indata.constraintweights.numpy() + syst_names = np.array(self.indata.systs).astype(bytes) + + selected = np.ones(nsyst, dtype=bool) + if skip_unconstrained: + selected &= cw > 0 + if include is not None: + keep = match_regexp_params(include, syst_names) + keep_set = set(keep) + selected &= np.array([n in keep_set for n in syst_names]) + if exclude is not None: + drop = match_regexp_params(exclude, syst_names) + drop_set = set(drop) + selected &= np.array([n not in drop_set for n in syst_names]) + + selected_idxs = np.where(selected)[0] + selected_names = syst_names[selected_idxs] + + logger.info( + f"global_asym_impacts_parms: selected {len(selected_idxs)}/{nsyst} " + f"nuisances (skip_unconstrained={skip_unconstrained})" + ) + + return global_asym_impacts.global_asym_impacts_parms( + self, + selected_idxs, + selected_names, + sigma=sigma, + linear_warmstart=linear_warmstart, + ) + def nonprofiled_impacts_parms(self, unconstrained_err=1.0): return nonprofiled_impacts.nonprofiled_impacts_parms( self.x, @@ -2869,6 +2928,7 @@ def _ensure_loss_grad(x): # Cheapest per iteration; may need more iterations to converge. # SR1 is more robust for non-convex local geometry than BFGS. if hess_mode == "hvp": + def scipy_hess(x, v): _ensure_loss_grad(x) n = len(x) diff --git a/rabbit/impacts/global_asym_impacts.py b/rabbit/impacts/global_asym_impacts.py new file mode 100644 index 0000000..f928d3a --- /dev/null +++ b/rabbit/impacts/global_asym_impacts.py @@ -0,0 +1,209 @@ +""" +Fully likelihood-based asymmetric global impacts. + +For each selected nuisance i, shift the auxiliary (global) observable theta0[i] +by +/- 1 prefit sigma and re-run the full fit. The resulting POI shifts are the +asymmetric global impacts. This is option (c) of the three definitions listed +in global_impacts.py and complements the Gaussian variants there: + + - global_impacts_parms / gaussian_global_impacts_parms: analytic, single-sided. + - global_asym_impacts_parms (this module): fully likelihood, two-sided, exact. + +In the Gaussian limit this reproduces gaussian_global_impacts_parms; deviations +measure the non-Gaussianity of the joint profile around the postfit minimum. + +Differs from nonprofiled_impacts in two ways: + - nuisance i is *profiled* (not frozen), so it can equilibrate at the new + constraint center together with all correlated nuisances; + - beta (BBB) parameters are profiled too, regardless of --noPostfitProfileBB. +""" + +import time + +import numpy as np +import tensorflow as tf +from wums import logging + +logger = logging.child_logger(__name__) + + +def _envelope(values): + """Quadrature envelope of asymmetric impacts within a group, separately for + the down and up sides. + + Args: + values: numpy array of shape (n_in_group, 2, n_total_params), where + axis 1 is [down, up]. + + Returns: + Array of shape (2, n_total_params). + """ + zeros = np.zeros((values.shape[0], values.shape[-1]), dtype=values.dtype) + vmin = np.min(values, axis=1) + vmax = np.max(values, axis=1) + lower = -np.sqrt(np.sum(np.minimum(zeros, vmin) ** 2, axis=0)) + upper = np.sqrt(np.sum(np.maximum(zeros, vmax) ** 2, axis=0)) + return np.stack([lower, upper]) + + +def global_asym_impacts_parms( + fitter, + selected_idxs, + selected_names, + sigma=1.0, + signs=(-1, 1), + linear_warmstart=False, +): + """Run a per-nuisance theta0-shift + re-fit and assemble the asymmetric + global impact tensor. + + Args: + fitter: the Fitter instance (used for x, theta0 and minimize). + selected_idxs: indices into the syst axis (0..nsyst-1) of nuisances to + scan. + selected_names: names of those nuisances (bytes), used as impact-axis + labels. + sigma: shift magnitude in units of the prefit constraint width + (constraints are unit-sigma in rabbit, so 1.0 = 1 prefit sigma). + signs: sequence (down, up). Bin 0 of axis_downUpVar -> first sign. + linear_warmstart: experimental. If True, warm-start each refit at + 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. + Requires fitter.cov to exist (same prerequisite as + --gaussianGlobalImpacts). + + Returns: + parms: np.ndarray of bytes, the impact-axis labels. + impacts: np.ndarray of shape (n_scanned, 2, n_total_params). + Axis 1 is [down, up] matching axis_downUpVar. + group_names: np.ndarray of bytes for groups containing scanned + nuisances (plus a trailing "Total"). + impacts_grouped: np.ndarray of shape (n_groups, 2, n_total_params). + """ + n_scanned = len(selected_idxs) + n_total = len(fitter.parms) + impacts = np.zeros((n_scanned, 2, n_total)) + + nparams = fitter.param_model.nparams + + # Snapshot postfit nominal state to restore between iterations. + x_nom = tf.identity(fitter.x.value()) + theta0_nom = tf.identity(fitter.theta0.value()) + theta0_nom_np = theta0_nom.numpy() + x_nom_np = x_nom.numpy() + + logger.info( + f"global_asym_impacts: shifting theta0 by +/- {sigma} sigma and " + f"re-fitting for {n_scanned} nuisances" + + (" (linear warm-start enabled)" if linear_warmstart else "") + ) + + # Optional Gaussian-approximation warm-start. + # dxdtheta0 has shape [npar, nsyst]; column i gives the linearised + # response of the postfit minimum to a unit shift of theta0[i]. Computing + # it once before the loop is the same cost as one --gaussianGlobalImpacts + # call. + dxdtheta0_np = None + if linear_warmstart: + if fitter.cov is None: + raise RuntimeError( + "global_asym_impacts: linear_warmstart requires fitter.cov " + "(incompatible with --noHessian)." + ) + t_lws = time.perf_counter() + dxdtheta0_tf, _, _ = fitter._dxdvars() + dxdtheta0_np = dxdtheta0_tf.numpy() + logger.info( + f"global_asym_impacts: dxdtheta0 prepared in " + f"{time.perf_counter() - t_lws:.2f}s" + ) + + t_per = np.zeros(n_scanned, dtype=np.float64) + t_total0 = time.perf_counter() + + for k, i in enumerate(selected_idxs): + i = int(i) + name = selected_names[k] + name_str = name.decode() if isinstance(name, bytes) else name + logger.info(f" [{k + 1}/{n_scanned}] theta0-shift refit for {name_str}") + + t0 = time.perf_counter() + for j, sign in enumerate(signs): + shift = float(sign) * float(sigma) + + # Always shift the constraint center for nuisance i by `shift`. + theta0_shifted = theta0_nom_np.copy() + theta0_shifted[i] += shift + + # Warm-start x. With linear_warmstart, use the Gaussian-approx new + # minimum x_nom + dxdtheta0[:, i] * shift -- on near-Gaussian + # nuisances this lands at the new minimum to within roundoff. + # Without it, just shift x[nparams+i] by `shift` so the nuisance + # itself starts at the new constraint center. + if linear_warmstart: + x_shifted = x_nom_np + dxdtheta0_np[:, i] * shift + else: + x_shifted = x_nom_np.copy() + x_shifted[nparams + i] += shift + + fitter.theta0.assign(theta0_shifted) + fitter.x.assign(x_shifted) + + try: + fitter.minimize() + if fitter.binByBinStat: + fitter._profile_beta() + impacts[k, j] = (fitter.x.value() - x_nom).numpy() + except Exception as e: + logger.warning( + f" refit for {name_str} sign={sign:+d} failed: {e}; " + "leaving impact at zero" + ) + + t_per[k] = time.perf_counter() - t0 + logger.info(f" took {t_per[k]:.2f}s") + + # Restore the fit state so downstream postfit computations see the nominal. + fitter.theta0.assign(theta0_nom) + fitter.x.assign(x_nom) + if fitter.binByBinStat: + fitter._profile_beta() + + if n_scanned > 0: + t_total = time.perf_counter() - t_total0 + logger.info( + f"global_asym_impacts: total {t_total:.1f}s " + f"(mean {t_per.mean():.2f}s, min {t_per.min():.2f}s, " + f"max {t_per.max():.2f}s per nuisance)" + ) + + # Grouped impacts via quadrature envelope, separately for down/up. + selected_set = set(int(idx) for idx in selected_idxs) + pos_in_scanned = {int(idx): k for k, idx in enumerate(selected_idxs)} + + group_names = [] + group_impacts = [] + for gname, gidxs in zip(fitter.indata.systgroups, fitter.indata.systgroupidxs): + gidxs = np.asarray(gidxs).astype(int) + in_scanned = [pos_in_scanned[i] for i in gidxs if int(i) in selected_set] + if not in_scanned: + continue + group_names.append(gname) + group_impacts.append(_envelope(impacts[in_scanned])) + + if n_scanned > 0: + group_names.append(b"Total") + group_impacts.append(_envelope(impacts)) + + if group_impacts: + impacts_grouped = np.stack(group_impacts) + else: + impacts_grouped = np.zeros((0, 2, n_total)) + + return ( + np.asarray(selected_names), + impacts, + np.asarray(group_names), + impacts_grouped, + ) diff --git a/tests/test_global_asym_impacts.py b/tests/test_global_asym_impacts.py new file mode 100644 index 0000000..4d4807f --- /dev/null +++ b/tests/test_global_asym_impacts.py @@ -0,0 +1,338 @@ +"""End-to-end correctness tests for --globalAsymImpacts. + +Builds a tensor with both symmetric and asymmetric nuisances via +tests/make_tensor.py, then runs three fits and checks: + + T1 Gaussian-limit closure: at sigma=0.01, globalAsym up/sigma matches + gaussianGlobalImpacts to ~1e-3 relative. + T2 Symmetry at small sigma: up = -down to ~1e-3 relative. + T3 Asymmetry detection: at sigma=1.0, the structurally asymmetric nuisance + slope_2_signal_ch1 has |up+down| above a threshold; structurally + symmetric nuisances stay close to zero. + T4 State restoration: postfit `parms` histogram values are bit-identical + between a baseline fit and the --globalAsymImpacts fit (same seed). + T5 Output schema: global_impacts_asym and *_grouped exist with axes + (impacts, downUpVar, parms). + T6 Unconstrained nuisances skipped: slope_background / slope_signal not + in the impacts axis. + +Run: source /opt/env.sh && source setup.sh + python tests/test_global_asym_impacts.py +""" + +from __future__ import annotations + +import os +import shutil +import subprocess +import sys +from pathlib import Path + +import numpy as np + +from rabbit import io_tools + +OUTDIR = Path("/tmp/test_global_asym_impacts") +TENSOR = OUTDIR / "test_tensor.hdf5" + +# Names taken from tests/make_tensor.py +# +# slope_2_signal_ch1 is added with symmetrize=None -> rabbit stores logkhalfdiff +# != 0, so this is the only nuisance for which down vs up should disagree at +# sigma=1 beyond minimizer noise. (slope_signal_ch0 also has kfactor=1.2 but +# is an NOI, treated separately by the model.) +ASYMMETRIC_NUIS = "slope_2_signal_ch1" +# The SymAvg/SymDiff partners produced by the symmetrize transforms have +# logkhalfdiff identically zero: structurally symmetric internal nuisances. +SYMMETRIC_NUIS = [ + "slope_lin_signal_ch0SymAvg", + "slope_lin_signal_ch0SymDiff", + "slope_quad_signal_ch0SymAvg", + "slope_quad_signal_ch0SymDiff", + "slope_signal_ch1", # symmetrize="conservative" -> single sym nuisance + "norm", # pure norm systematics + "bkg_norm", + "bkg_2_norm", +] +UNCONSTRAINED_NUIS = ["slope_background", "slope_signal"] + + +def run(cmd: list[str]) -> None: + print("$", " ".join(cmd), flush=True) + subprocess.run(cmd, check=True) + + +def make_tensor() -> None: + OUTDIR.mkdir(parents=True, exist_ok=True) + run([sys.executable, "tests/make_tensor.py", "-o", str(OUTDIR) + "/"]) + assert TENSOR.exists(), TENSOR + + +def fit(label: str, *extra: str) -> Path: + out = OUTDIR / label + if out.exists(): + shutil.rmtree(out) + out.mkdir() + run( + [ + "rabbit_fit.py", + str(TENSOR), + "-o", + str(out), + "--seed", + "42", + *extra, + ] + ) + return out / "fitresults.hdf5" + + +def load_results(path: Path) -> dict: + return io_tools.get_fitresult(str(path)) + + +def get_impact_array(res: dict, name: str) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """Return (values, parms_axis, impacts_axis) for an impacts hist.""" + h = res[name].get() + impacts_axis = np.array(h.axes["impacts"]) + parms_axis = np.array(h.axes["parms"]) + return h.values(), parms_axis, impacts_axis + + +def get_nuis_idx(impacts_axis: np.ndarray, name: str) -> int: + matches = np.where(impacts_axis.astype(str) == name)[0] + if len(matches) == 0: + raise KeyError( + f"{name!r} not in impacts axis: {list(impacts_axis.astype(str))}" + ) + return int(matches[0]) + + +def get_poi_idx(parms_axis: np.ndarray, name: str) -> int: + return int(np.where(parms_axis.astype(str) == name)[0][0]) + + +# --------------------------------------------------------------------------- + + +def t1_gaussian_closure(res_small: dict, res_gauss: dict) -> None: + """At sigma=0.01, globalAsym up/sigma must match gaussianGlobal.""" + print("\n[T1] Gaussian-limit closure (sigma=0.01)") + + sigma = 0.01 + asym, asym_parms, asym_impacts = get_impact_array(res_small, "global_impacts_asym") + # asym shape: (impacts, downUpVar, parms) + # gauss shape: (parms, impacts) + gauss_h = res_gauss["gaussian_global_impacts"].get() + g_vals = gauss_h.values() + g_parms = np.array(gauss_h.axes["parms"]).astype(str) + g_impacts = np.array(gauss_h.axes["impacts"]).astype(str) + + poi = "sig" # the only POI in make_tensor.py default config + poi_idx_asym = get_poi_idx(asym_parms, poi) + poi_idx_gauss = int(np.where(g_parms == poi)[0][0]) + + common = [n for n in asym_impacts.astype(str) if n in g_impacts] + print(f" comparing {len(common)} nuisances on POI {poi!r}") + + max_rel = 0.0 + worst = None + for name in common: + i_asym = get_nuis_idx(asym_impacts, name) + i_gauss = int(np.where(g_impacts == name)[0][0]) + up = asym[i_asym, 1, poi_idx_asym] / sigma + ref = g_vals[poi_idx_gauss, i_gauss] + # gaussian impact is unsigned; compare magnitudes + denom = max(abs(ref), 1e-12) + rel = abs(abs(up) - abs(ref)) / denom + if rel > max_rel: + max_rel, worst = rel, (name, up, ref) + + print(f" max relative diff: {max_rel:.3e} (worst: {worst})") + assert max_rel < 5e-2, ( + f"T1 FAIL: globalAsym at sigma=0.01 deviates from gaussianGlobal by " + f"{max_rel:.3e} on {worst!r}; expected closure <5%" + ) + print(" T1 PASS") + + +def t2_small_sigma_symmetry(res_small: dict) -> None: + """At sigma=0.01, down ~ -up for every scanned nuisance.""" + print("\n[T2] Small-sigma symmetry (sigma=0.01: up = -down)") + + asym, asym_parms, asym_impacts = get_impact_array(res_small, "global_impacts_asym") + poi_idx = get_poi_idx(asym_parms, "sig") + + max_rel = 0.0 + worst = None + for k, name in enumerate(asym_impacts.astype(str)): + down = asym[k, 0, poi_idx] + up = asym[k, 1, poi_idx] + denom = max(abs(up), abs(down), 1e-12) + rel = abs(up + down) / denom + if rel > max_rel: + max_rel, worst = rel, (name, up, down) + + print(f" max |up+down|/max(|up|,|down|): {max_rel:.3e} (worst: {worst})") + assert max_rel < 5e-2, f"T2 FAIL: not symmetric at small sigma; worst={worst}" + print(" T2 PASS") + + +def t3_asymmetry_at_full_sigma(res_full: dict, res_small: dict) -> None: + """At sigma=1.0 the method must produce non-trivial asymmetry, and the + asymmetry must scale with sigma. + + Note: even structurally symmetric nuisances (logkhalfdiff=0) can show + asymmetric impacts at sigma=1.0 with systematic_type="log_normal", because + the kernel exp(+theta * logkavg) is asymmetric in theta away from the + quadratic basin. This is correct behaviour, not a bug. We therefore assert: + a) the structurally asymmetric nuisance shows clearly non-zero + asymmetry at sigma=1.0; + b) every nuisance that was symmetric at sigma=0.01 (validated by T2) + shows MORE asymmetry at sigma=1.0 than at sigma=0.01 -- i.e. the + asymmetry is a real sigma-dependent effect, not noise. + """ + print("\n[T3] Asymmetry detection at sigma=1.0") + + asym_full, asym_parms, asym_impacts = get_impact_array( + res_full, "global_impacts_asym" + ) + asym_small, _, asym_impacts_small = get_impact_array( + res_small, "global_impacts_asym" + ) + poi_idx = get_poi_idx(asym_parms, "sig") + + def asymmetry(arr: np.ndarray, axis: np.ndarray, name: str) -> float: + k = get_nuis_idx(axis, name) + d, u = arr[k, 0, poi_idx], arr[k, 1, poi_idx] + return abs(u + d) / max(abs(u), abs(d), 1e-12) + + # (a) structurally asymmetric nuisance must show measurable asymmetry. + impacts_str = asym_impacts.astype(str) + if ASYMMETRIC_NUIS not in impacts_str: + print( + f" skip (a): {ASYMMETRIC_NUIS} not in impacts axis; " + f"list = {list(impacts_str)}" + ) + else: + a = asymmetry(asym_full, asym_impacts, ASYMMETRIC_NUIS) + print(f" asymmetry({ASYMMETRIC_NUIS}) at sigma=1.0 = {a:.3e}") + assert a > 1e-2, ( + f"T3 FAIL (a): {ASYMMETRIC_NUIS} should be measurably asymmetric " + f"at sigma=1.0, got {a:.3e}" + ) + + # (b) every nuisance: asymmetry grows with sigma (within minimizer noise). + growth_floor = 1e-3 + for name in asym_impacts.astype(str): + if name not in asym_impacts_small.astype(str): + continue + a_small = asymmetry(asym_small, asym_impacts_small, name) + a_full = asymmetry(asym_full, asym_impacts, name) + print( + f" asymmetry({name}): sigma=0.01 -> {a_small:.3e}, " + f"sigma=1.0 -> {a_full:.3e}" + ) + # Allow noise floor: even if the asymmetry is tiny, it should not + # *shrink* with sigma if sigma=0.01 was already at convergence noise. + assert a_full + growth_floor >= a_small, ( + f"T3 FAIL (b): {name} asymmetry shrinks with sigma " + f"({a_small:.3e} -> {a_full:.3e})" + ) + print(" T3 PASS") + + +def t4_state_restoration(res_baseline: dict, res_full: dict) -> None: + """After the asym scan, postfit `parms` must match a vanilla fit's parms.""" + print("\n[T4] State restoration (parms unchanged vs baseline)") + + p0 = res_baseline["parms"].get().values() + p1 = res_full["parms"].get().values() + diff = np.max(np.abs(p0 - p1)) + print(f" max |parms_baseline - parms_after_asym| = {diff:.3e}") + # Same seed, same fit -> should be bit-identical, but allow a tiny float + # margin in case of any nondeterminism in the env. + assert diff < 1e-9, "T4 FAIL: postfit parms changed after the asym loop" + print(" T4 PASS") + + +def t5_output_schema(res_full: dict) -> None: + print("\n[T5] Output schema") + assert "global_impacts_asym" in res_full, list(res_full.keys()) + assert "global_impacts_asym_grouped" in res_full, list(res_full.keys()) + + h = res_full["global_impacts_asym"].get() + names = [a.name for a in h.axes] + print(f" axes = {names}") + assert names == ["impacts", "downUpVar", "parms"], names + + g = res_full["global_impacts_asym_grouped"].get() + gnames = [a.name for a in g.axes] + print(f" grouped axes = {gnames}") + assert gnames == ["impacts", "downUpVar", "parms"], gnames + + groups = list(np.array(g.axes["impacts"]).astype(str)) + print(f" groups = {groups}") + assert groups[-1] == "Total", groups + print(" T5 PASS") + + +def t6_unconstrained_skipped(res_full: dict) -> None: + print("\n[T6] Unconstrained nuisances skipped") + _, _, impacts_axis = get_impact_array(res_full, "global_impacts_asym") + impacts_str = list(impacts_axis.astype(str)) + print(f" scanned nuisances = {impacts_str}") + for n in UNCONSTRAINED_NUIS: + assert n not in impacts_str, f"T6 FAIL: unconstrained {n!r} should be skipped" + print(" T6 PASS") + + +# --------------------------------------------------------------------------- + + +def main() -> None: + os.chdir(os.environ.get("RABBIT_BASE", ".")) + make_tensor() + + common = [ + "--doImpacts", + "--gaussianGlobalImpacts", + ] + + # baseline: just the fit, no asym scan -> reference parms + base_h5 = fit("baseline", *common) + + # small-sigma asym for Gaussian closure & symmetry + small_h5 = fit( + "small_sigma", + *common, + "--globalAsymImpacts", + "--globalAsymImpactsSigma", + "0.01", + ) + + # full-sigma asym for asymmetry detection + full_h5 = fit( + "full_sigma", + *common, + "--globalAsymImpacts", + "--globalAsymImpactsSigma", + "1.0", + ) + + res_base = load_results(base_h5) + res_small = load_results(small_h5) + res_full = load_results(full_h5) + + t5_output_schema(res_full) + t6_unconstrained_skipped(res_full) + t1_gaussian_closure(res_small, res_small) # gaussianGlobal is in same file + t2_small_sigma_symmetry(res_small) + t3_asymmetry_at_full_sigma(res_full, res_small) + t4_state_restoration(res_base, res_full) + + print("\nAll tests passed.") + + +if __name__ == "__main__": + main() From 843fbbd0a9bbf48ab7726e0cc9cb0a9e5629616a Mon Sep 17 00:00:00 2001 From: Luca Lavezzo Date: Fri, 29 May 2026 16:40:58 -0400 Subject: [PATCH 4/7] improve timinig and logging; --noEDM; minimizer settings via CLI; gaussian constraints on params in ParamModels --- bin/rabbit_fit.py | 73 ++++++++++++--- rabbit/fitter.py | 141 ++++++++++++++++++++++++++++- rabbit/param_models/param_model.py | 10 ++ rabbit/parsing.py | 43 +++++++++ 4 files changed, 249 insertions(+), 18 deletions(-) diff --git a/bin/rabbit_fit.py b/bin/rabbit_fit.py index e30c145..b40f670 100755 --- a/bin/rabbit_fit.py +++ b/bin/rabbit_fit.py @@ -566,17 +566,26 @@ def fit(args, fitter, ws, dofit=True): ws.add_1D_integer_hist([max_curvature], "best", "lcurve") if dofit: + _t_min = time.perf_counter() cb = fitter.minimize() + logger.debug(f"[timing] fitter.minimize(): {time.perf_counter() - _t_min:.1f} s") # force profiling of beta with final parameter values # TODO avoid the extra calculation and jitting if possible since the relevant calculation # usually would have been done during the minimization if fitter.binByBinStat and not args.noPostfitProfileBB: + _t_bb = time.perf_counter() fitter._profile_beta() + logger.debug(f"[timing] _profile_beta(): {time.perf_counter() - _t_bb:.1f} s") if cb is not None: ws.add_1D_integer_hist(cb.loss_history, "epoch", "loss") ws.add_1D_integer_hist(cb.time_history, "epoch", "time") + try: + n_iter = len(cb.loss_history) + logger.debug(f"[timing] L-BFGS: {n_iter} epochs recorded") + except Exception: + pass # prefit variances as the default fallback for add_parms_hist below parms_variances = None @@ -628,21 +637,33 @@ def fit(args, fitter, ws, dofit=True): # Hessian-vector products. The CG solves touch O(npar) memory # per call instead of O(npar^2), so this works on problems # where the full covariance would be infeasible. - _, grad = fitter.loss_val_grad() - npoi = int(fitter.poi_model.npoi) - noi_idx_in_x = np.asarray(fitter.indata.noiidxs, dtype=np.int64) + npoi - poi_noi_idx = np.concatenate([np.arange(npoi, dtype=np.int64), noi_idx_in_x]) - edmval, cov_rows = fitter.edmval_cov_rows_hessfree(grad, poi_noi_idx) - logger.info(f"edmval: {edmval}") - - # Build a full-length variance vector with the POI+NOI entries - # populated from the diagonal of the CG-solved rows and the rest - # left as NaN (we did not compute those). add_parms_hist stores - # the vector verbatim into the workspace. n = int(fitter.x.shape[0]) - parms_variances_np = np.full(n, np.nan, dtype=np.float64) - for k, i in enumerate(poi_noi_idx): - parms_variances_np[int(i)] = cov_rows[k, int(i)] + if getattr(args, "noEDM", False): + logger.debug("[timing] --noEDM set: skipping Hessian-free EDM/CG step") + edmval = None + parms_variances_np = np.full(n, np.nan, dtype=np.float64) + else: + _t_lg = time.perf_counter() + _, grad = fitter.loss_val_grad() + logger.debug(f"[timing] loss_val_grad() (postfit): {time.perf_counter() - _t_lg:.1f} s") + + npoi = int(fitter.param_model.npoi) + noi_idx_in_x = np.asarray(fitter.indata.noiidxs, dtype=np.int64) + npoi + poi_noi_idx = np.concatenate([np.arange(npoi, dtype=np.int64), noi_idx_in_x]) + logger.debug(f"[timing] EDM CG: solving for {len(poi_noi_idx)} POI+NOI rows") + + _t_edm = time.perf_counter() + edmval, cov_rows = fitter.edmval_cov_rows_hessfree(grad, poi_noi_idx) + logger.debug(f"[timing] edmval_cov_rows_hessfree(): {time.perf_counter() - _t_edm:.1f} s") + logger.info(f"edmval: {edmval}") + + # Build a full-length variance vector with the POI+NOI entries + # populated from the diagonal of the CG-solved rows and the rest + # left as NaN (we did not compute those). add_parms_hist stores + # the vector verbatim into the workspace. + parms_variances_np = np.full(n, np.nan, dtype=np.float64) + for k, i in enumerate(poi_noi_idx): + parms_variances_np[int(i)] = cov_rows[k, int(i)] parms_variances = tf.constant(parms_variances_np, dtype=fitter.indata.dtype) nllvalreduced = fitter.reduced_nll().numpy() @@ -888,6 +909,17 @@ def main(): "nois": ifitter.parms[ifitter.param_model.nparams :][indata.noiidxs], } + # ParamModel Gaussian priors (if --paramModelPriors was used and the + # model declared sigmas). Stored as a small dict so downstream tooling + # can see what was applied without parsing the rabbit log. + if getattr(ifitter, "param_prior_active", False): + meta["param_priors"] = { + "params": ifitter.param_model.params, # all nparams names + "mask": ifitter.param_prior_mask_np, # bool array + "sigmas": ifitter.param_prior_sigmas_np, # NaN where mask False + "means": ifitter.param_prior_means_np, # NaN where mask False + } + with workspace.Workspace( args.outpath, args.outname, @@ -922,6 +954,19 @@ def main(): for j, (data_values, data_variances) in enumerate(datasets): ifitter.defaultassign() + # If --externalPostfit was supplied AND this is an + # asimov toy (-t -1), load the postfit values BEFORE we + # compute the expected yield. Otherwise expected_yield() + # below is evaluated at the prefit point (because + # defaultassign() just reset x), so we'd get a + # prefit-Asimov regardless of --externalPostfit. The + # in-fit() load (around line 551) then re-runs harmlessly. + if ifit == -1 and args.externalPostfit is not None: + ifitter.load_fitresult( + args.externalPostfit, + args.externalPostfitResult, + profile=not args.noPostfitProfileBB, + ) if ifit == -1: group.append("asimov") ifitter.set_nobs(ifitter.expected_yield()) diff --git a/rabbit/fitter.py b/rabbit/fitter.py index 8c8bb75..40290ce 100644 --- a/rabbit/fitter.py +++ b/rabbit/fitter.py @@ -66,7 +66,14 @@ def __init__(self, xv, early_stopping=-1): def __call__(self, intermediate_result): loss = intermediate_result.fun - logger.debug(f"Iteration {self.iiter}: loss value {loss}") + elapsed = time.time() - self.t0 + prev = self.time_history[-1] if self.time_history else 0.0 + dt = elapsed - prev + + logger.debug( + f"Iteration {self.iiter}: loss {loss} " + f"[dt={dt:.2f}s elapsed={elapsed:.2f}s]" + ) if np.isnan(loss): raise ValueError(f"Loss value is NaN at iteration {self.iiter}") @@ -80,7 +87,7 @@ def __call__(self, intermediate_result): ) self.loss_history.append(loss) - self.time_history.append(time.time() - self.t0) + self.time_history.append(elapsed) self.xval = intermediate_result.x self.iiter += 1 @@ -131,6 +138,12 @@ def __init__( self.diagnostics = options.diagnostics self.minimizer_method = options.minimizerMethod + # Optional scipy.optimize.minimize tolerances. None entries are + # skipped at call site so each method falls back to scipy defaults + # for what's not set. + self.minimizer_maxiter = getattr(options, "minimizerMaxiter", None) + self.minimizer_gtol = getattr(options, "minimizerGtol", None) + self.minimizer_ftol = getattr(options, "minimizerFtol", None) self.hvp_method = getattr(options, "hvpMethod", "revrev") # jitCompile accepts "auto" (the default), "on", or "off". # True / False from programmatic callers are accepted as @@ -190,6 +203,7 @@ def __init__( unblind=options.unblind, blinding_group=getattr(options, "blindingGroup", []), freeze_parameters=options.freezeParameters, + options=options, ) # --- observed number of events per bin @@ -312,6 +326,7 @@ def init_fit_parms( unblind=False, blinding_group=[], freeze_parameters=[], + options=None, ): self.param_model = param_model @@ -391,6 +406,75 @@ def init_fit_parms( self.indata.dtype, ) + # ParamModel Gaussian priors (opt-in via --paramModelPriors). + # When enabled, the Fitter reads two optional attributes from the + # ParamModel: + # - prior_sigmas: np.ndarray shape (nparams,). Entries that are + # finite and > 0 are Gaussian-constrained at that width; NaN / + # non-finite / 0 entries leave the corresponding parameter free. + # - prior_means : np.ndarray shape (nparams,), optional. + # Defaults to param_model.xparamdefault when not provided. + # The penalty 0.5 * (x - μ)^2 / σ^2 is added to _compute_lc. + self.param_prior_active = False + self.param_prior_mask_np = None + self.param_prior_sigmas_np = None + self.param_prior_means_np = None + if getattr(options, "paramModelPriors", False): + pm = self.param_model + sigmas = getattr(pm, "prior_sigmas", None) + if sigmas is not None: + sigmas = np.asarray(sigmas, dtype=np.float64) + if sigmas.shape != (pm.nparams,): + raise ValueError( + f"param_model.prior_sigmas must have shape ({pm.nparams},); " + f"got {sigmas.shape}" + ) + means = getattr(pm, "prior_means", None) + if means is None: + means = pm.xparamdefault.numpy().astype(np.float64) + else: + means = np.asarray(means, dtype=np.float64) + if means.shape != (pm.nparams,): + raise ValueError( + f"param_model.prior_means must have shape ({pm.nparams},); " + f"got {means.shape}" + ) + mask = np.isfinite(sigmas) & (sigmas > 0) + n_priored = int(mask.sum()) + if n_priored > 0: + self.param_prior_active = True + self.param_prior_mask_np = mask.copy() + self.param_prior_sigmas_np = np.where(mask, sigmas, np.nan) + self.param_prior_means_np = np.where(mask, means, np.nan) + inv2_np = np.where(mask, 1.0 / sigmas**2, 0.0) + self.param_prior_mask = tf.constant(mask, dtype=tf.bool) + self.param_prior_inv2 = tf.constant(inv2_np, dtype=self.indata.dtype) + self.param_prior_means = tf.constant( + np.where(mask, means, 0.0), dtype=self.indata.dtype + ) + self.param_prior_log_sigma2 = tf.constant( + np.where(mask, np.log(np.where(mask, sigmas, 1.0) ** 2), 0.0), + dtype=self.indata.dtype, + ) + logger.info( + f"[paramPriors] applying Gaussian priors to " + f"{n_priored}/{pm.nparams} ParamModel params:" + ) + for i, p in enumerate(pm.params): + if mask[i]: + name = p.decode() if isinstance(p, bytes) else str(p) + logger.info(f" {name}: μ={means[i]:.4g} σ={sigmas[i]:.4g}") + else: + logger.info( + "[paramPriors] --paramModelPriors set but prior_sigmas " + "has no finite entries; no priors applied." + ) + else: + logger.info( + "[paramPriors] --paramModelPriors set but ParamModel does " + "not declare prior_sigmas; no priors applied." + ) + # constraint minima for nuisance parameters self.theta0 = tf.Variable( self.theta0default, @@ -655,7 +739,17 @@ def get_theta(self): def get_model_nui(self): npoi = self.param_model.npoi npou = self.param_model.npou - return self.x[npoi : npoi + npou] + nui = self.x[npoi : npoi + npou] + # Apply frozen_params_mask the same way get_poi() and get_theta() do. + # Without this, --freezeParameters silently fails to freeze POUs: + # the param is registered as frozen but no stop_gradient is applied, + # so the optimizer updates it anyway. + nui = tf.where( + self.frozen_params_mask[npoi : npoi + npou], + tf.stop_gradient(nui), + nui, + ) + return nui def get_poi(self): xpoi = self.x[: self.param_model.npoi] @@ -2382,7 +2476,28 @@ def _compute_lc(self, full_nll=False): # normalization factor for normal distribution: log(1/sqrt(2*pi)) = -0.9189385332046727 lc = lc + 0.9189385332046727 * self.indata.constraintweights - return tf.reduce_sum(lc) + total = tf.reduce_sum(lc) + + # ParamModel Gaussian priors (opt-in via --paramModelPriors). + # Per-entry mask is baked into param_prior_inv2 (zero where mask is + # False), so this is a single elementwise op with no branching. + if self.param_prior_active: + pm_params = tf.concat( + [self.get_poi(), self.get_model_nui()], axis=0 + ) + lc_pm = 0.5 * self.param_prior_inv2 * tf.square( + pm_params - self.param_prior_means + ) + if full_nll: + # 0.5*log(2π σ²) on entries with a prior. The log(σ²) term is + # precomputed (zero where masked off). + mask_dtype = tf.cast(self.param_prior_mask, lc_pm.dtype) + lc_pm = lc_pm + mask_dtype * ( + 0.9189385332046727 + 0.5 * self.param_prior_log_sigma2 + ) + total = total + tf.reduce_sum(lc_pm) + + return total def _compute_lbeta(self, beta, full_nll=False): if self.binByBinStat: @@ -2671,6 +2786,23 @@ def scipy_hess(xval): else: info_minimize = dict() + # Build scipy.optimize.minimize options from --minimizerMaxiter/ + # --minimizerGtol/--minimizerFtol. Anything not set explicitly falls + # back to the historical tol=0.0 baseline below (run to the tightest + # internal criteria), since `tol` only fills options keys that are not + # already present. Methods that don't recognize an option ignore it + # with an OptimizeWarning (no crash). + sci_opts = {} + if self.minimizer_maxiter is not None: + sci_opts["maxiter"] = int(self.minimizer_maxiter) + if self.minimizer_gtol is not None: + sci_opts["gtol"] = float(self.minimizer_gtol) + if self.minimizer_ftol is not None: + sci_opts["ftol"] = float(self.minimizer_ftol) + logger.info( + f"[minimize] method={self.minimizer_method} options={sci_opts}" + ) + try: res = scipy.optimize.minimize( scipy_loss, @@ -2679,6 +2811,7 @@ def scipy_hess(xval): jac=True, tol=0.0, callback=callback, + options=sci_opts, **info_minimize, ) except Exception as ex: diff --git a/rabbit/param_models/param_model.py b/rabbit/param_models/param_model.py index abd0c70..3752daf 100644 --- a/rabbit/param_models/param_model.py +++ b/rabbit/param_models/param_model.py @@ -17,6 +17,16 @@ def __init__(self, indata, *args, **kwargs): # self.xparamdefault = # default values for all parameters (length nparams) # self.is_linear = # define if the model is linear in the parameters # self.allowNegativeParam = # define if the POI parameters can be negative or not + # + # # optional: Gaussian priors on the model's parameters + # # consumed by Fitter when --paramModelPriors is set (otherwise ignored). + # self.prior_sigmas = # np.ndarray, shape (nparams,). Entries that are + # # finite and > 0 are Gaussian-constrained at that + # # width; NaN / non-finite / 0 entries leave the + # # corresponding parameter free. Convention is up + # # to the model (e.g. POIs free, POUs constrained). + # self.prior_means = # np.ndarray, shape (nparams,). Optional; defaults + # # to self.xparamdefault when not provided. @property def nparams(self): diff --git a/rabbit/parsing.py b/rabbit/parsing.py index 6e8a0b8..963c29a 100644 --- a/rabbit/parsing.py +++ b/rabbit/parsing.py @@ -240,6 +240,49 @@ def common_parser(): action="store_true", help="Don't compute the hessian of parameters", ) + parser.add_argument( + "--noEDM", + default=False, + action="store_true", + help="Skip the Hessian-free EDM/CG postfit step (only meaningful with " + "--noHessian). Skips both the edmval estimate and the POI+NOI " + "uncertainty rows; their entries in the output are NaN.", + ) + parser.add_argument( + "--paramModelPriors", + default=False, + action="store_true", + help="Opt in to Gaussian priors on ParamModel parameters. When set, " + "the Fitter reads param_model.prior_sigmas (np.ndarray of shape " + "(nparams,), entries finite > 0 → prior of that width, NaN → free) " + "and param_model.prior_means (defaults to param_model.xparamdefault) " + "and adds the corresponding penalty to the NLL constraint term. By " + "default the feature is off and ParamModel params float free.", + ) + parser.add_argument( + "--minimizerMaxiter", + type=int, + default=None, + help="Cap the number of scipy.optimize.minimize iterations. " + "Passed as options={'maxiter': N} to the minimizer. None (default) " + "uses scipy's method-specific default (typically 1000+).", + ) + parser.add_argument( + "--minimizerGtol", + type=float, + default=None, + help="Gradient-norm tolerance for the minimizer. Passed as " + "options={'gtol': X} (recognized by trust-krylov, L-BFGS-B, BFGS, CG). " + "None (default) uses scipy's per-method default.", + ) + parser.add_argument( + "--minimizerFtol", + type=float, + default=None, + help="Function-value tolerance for the minimizer (L-BFGS-B only). " + "Passed as options={'ftol': X}. None (default) uses scipy's default " + "(2.22e-9 for L-BFGS-B). Ignored by methods that don't honor 'ftol'.", + ) parser.add_argument( "--prefitUnconstrainedNuisanceUncertainty", default=0.0, From b2dd17f4e31b4d366c7d39fc18374205c24bfea9 Mon Sep 17 00:00:00 2001 From: Luca Lavezzo Date: Wed, 10 Jun 2026 16:26:15 -0400 Subject: [PATCH 5/7] changes for NP scetlib model --- bin/rabbit_print_pulls_and_constraints.py | 37 ++++++++++----- rabbit/fitter.py | 50 +++++++++++++++++++-- rabbit/impacts/traditional_impacts.py | 55 ++++++++++++++++++----- rabbit/workspace.py | 18 +++++++- 4 files changed, 132 insertions(+), 28 deletions(-) diff --git a/bin/rabbit_print_pulls_and_constraints.py b/bin/rabbit_print_pulls_and_constraints.py index 2f435b0..be4430f 100755 --- a/bin/rabbit_print_pulls_and_constraints.py +++ b/bin/rabbit_print_pulls_and_constraints.py @@ -50,6 +50,12 @@ def make_parser(): default=None, help="Exclude nuisances by regular expression", ) + parser.add_argument( + "--noPrefit", + default=False, + action="store_true", + help="Suppress the prefit pull and constraint columns from the printout", + ) return parser @@ -109,15 +115,22 @@ def main(): constraints_asym = constraints_asym[order] nround = 5 + prefit_header = ( + "" if args.noPrefit else f" ({'pull prefit':>11} +/- {'constraint prefit':>17})" + ) if args.asym: - print( - f" {'Parameter':<30} {'pull':>6} +/- {'constraint':>10} + {'up':>10} - {'down':>10} ({'pull prefit':>11} +/- {'constraint prefit':>17})" - ) - print(" " + "-" * 100) + header = f" {'Parameter':<30} {'pull':>6} +/- {'constraint':>10} + {'up':>10} - {'down':>10}{prefit_header}" + print(header) + print(" " + "-" * (len(header) - 3)) print( "\n".join( [ - f" {l:<30} {round(p, nround):>6} +/- {round(c, nround):>10} + {round(c_asym[0], nround):>10} - {round(c_asym[1], nround):>10} ({round(pp, nround):>11} +/- {round(pc, nround):>17})" + f" {l:<30} {round(p, nround):>6} +/- {round(c, nround):>10} + {round(c_asym[0], nround):>10} - {round(c_asym[1], nround):>10}" + + ( + "" + if args.noPrefit + else f" ({round(pp, nround):>11} +/- {round(pc, nround):>17})" + ) for l, p, c, c_asym, pp, pc in zip( labels, pulls, @@ -130,14 +143,18 @@ def main(): ) ) else: - print( - f" {'Parameter':<30} {'pull':>6} +/- {'constraint':>10} ({'pull prefit':>11} +/- {'constraint prefit':>17})" - ) - print(" " + "-" * 100) + header = f" {'Parameter':<30} {'pull':>6} +/- {'constraint':>10}{prefit_header}" + print(header) + print(" " + "-" * (len(header) - 3)) print( "\n".join( [ - f" {l:<30} {round(p, nround):>6} +/- {round(c, nround):>10} ({round(pp, nround):>11} +/- {round(pc, nround):>17})" + f" {l:<30} {round(p, nround):>6} +/- {round(c, nround):>10}" + + ( + "" + if args.noPrefit + else f" ({round(pp, nround):>11} +/- {round(pc, nround):>17})" + ) for l, p, c, pp, pc in zip( labels, pulls, constraints, pulls_prefit, constraints_prefit ) diff --git a/rabbit/fitter.py b/rabbit/fitter.py index 40290ce..c794c59 100644 --- a/rabbit/fitter.py +++ b/rabbit/fitter.py @@ -1158,6 +1158,47 @@ def _hvp_np(p_np): return edmval, cov_rows + def _resolved_param_impact_groups(self): + """ + ParamModel impact groups resolved to floating full-x parameter indices. + """ + groups = getattr(self.param_model, "param_impact_groups", None) + if not groups: + return [] + parms = self.parms.astype(str) + name_to_idx = {p: i for i, p in enumerate(parms)} + frozen = set(int(i) for i in np.atleast_1d(self.frozen_indices)) + resolved = [] + for label, pnames in groups.items(): + idxs = [ + name_to_idx[p] + for p in pnames + if p in name_to_idx and name_to_idx[p] not in frozen + ] + if idxs: + resolved.append((label, np.array(idxs, dtype=np.int32))) + return resolved + + def _cov_stat_floating(self, hess, nstat): + """ + Invert the stat sub-Hessian hess[:nstat, :nstat], excluding frozen params. + """ + hess_stat = hess[:nstat, :nstat] + if len(self.frozen_params) == 0: + return tf.linalg.inv(hess_stat) + stat_float = self.floating_indices[self.floating_indices < nstat] + sub = tf.gather(tf.gather(hess_stat, stat_float, axis=0), stat_float, axis=1) + sub_inv = tf.linalg.inv(sub) + coords = tf.reshape( + tf.stack(tf.meshgrid(stat_float, stat_float, indexing="ij"), axis=-1), + [-1, 2], + ) + return tf.scatter_nd( + tf.cast(coords, tf.int64), + tf.reshape(sub_inv, [-1]), + tf.constant([nstat, nstat], dtype=tf.int64), + ) + @tf.function def impacts_parms(self, hess): @@ -1166,18 +1207,17 @@ def impacts_parms(self, hess): + self.param_model.npou + self.indata.nsystnoconstraint ) - hess_stat = hess[:nstat, :nstat] - cov_stat = tf.linalg.inv(hess_stat) + cov_stat = self._cov_stat_floating(hess, nstat) if self.binByBinStat: val_no_bbb, grad_no_bbb, hess_no_bbb = self.loss_val_grad_hess( profile=False ) - hess_stat_no_bbb = hess_no_bbb[:nstat, :nstat] - cov_stat_no_bbb = tf.linalg.inv(hess_stat_no_bbb) + cov_stat_no_bbb = self._cov_stat_floating(hess_no_bbb, nstat) else: cov_stat_no_bbb = None + param_groups = self._resolved_param_impact_groups() impacts, impacts_grouped = traditional_impacts.impacts_parms( self.cov, cov_stat, @@ -1185,6 +1225,8 @@ def impacts_parms(self, hess): self.param_model.npoi, self.indata.noiidxs, self.indata.systgroupidxs, + nmodel_params=self.param_model.npoi + self.param_model.npou, + param_groupidxs=[idxs for _, idxs in param_groups], ) return impacts, impacts_grouped diff --git a/rabbit/impacts/traditional_impacts.py b/rabbit/impacts/traditional_impacts.py index f02cdd1..8236cd6 100644 --- a/rabbit/impacts/traditional_impacts.py +++ b/rabbit/impacts/traditional_impacts.py @@ -16,58 +16,76 @@ def _compute_impact_group(cov, v, idxs, nsignal_params=0): return tf.sqrt(v_invC_v) -def _gather_poi_noi_vector(v, noiidxs, nsignal_params=0): +def _gather_poi_noi_vector(v, noiidxs, nsignal_params=0, nmodel_params=None): + # nmodel_params = npoi + npou + if nmodel_params is None: + nmodel_params = nsignal_params v_poi = v[:nsignal_params] # protection for constained NOIs, set them to 0 - mask = (noiidxs >= 0) & (noiidxs < tf.shape(v[nsignal_params:])[0]) + mask = (noiidxs >= 0) & (noiidxs < tf.shape(v[nmodel_params:])[0]) safe_idxs = tf.where(mask, noiidxs, 0) mask = tf.cast(mask, v.dtype) mask = tf.reshape( mask, tf.concat([tf.shape(mask), tf.ones(tf.rank(v) - 1, dtype=tf.int32)], axis=0), ) - v_noi = tf.gather(v[nsignal_params:], safe_idxs) * mask + v_noi = tf.gather(v[nmodel_params:], safe_idxs) * mask v_gathered = tf.concat([v_poi, v_noi], axis=0) return v_gathered def impacts_parms( - cov, cov_stat, cov_stat_no_bbb, nsignal_params=0, noiidxs=[], systgroupidxs=[] + cov, + cov_stat, + cov_stat_no_bbb, + nsignal_params=0, + noiidxs=[], + systgroupidxs=[], + nmodel_params=None, + param_groupidxs=None, ): """ Gaussian approximation """ + if nmodel_params is None: + nmodel_params = nsignal_params # impact for poi at index i in covariance matrix from nuisance with index j is C_ij/sqrt(C_jj) = /sqrt() - v = _gather_poi_noi_vector(cov, noiidxs, nsignal_params) - impacts = v / tf.reshape(tf.sqrt(tf.linalg.diag_part(cov)), [1, -1]) + v = _gather_poi_noi_vector(cov, noiidxs, nsignal_params, nmodel_params) + # Frozen / fixed parameters have zero variance (cov diag == 0); their impact + # is zero, not 0/0 = nan. Guard the denominator so those columns read 0. + sigma = tf.sqrt(tf.linalg.diag_part(cov)) + safe_sigma = tf.where(sigma > 0, sigma, tf.ones_like(sigma)) + impacts = v / tf.reshape(safe_sigma, [1, -1]) if cov_stat_no_bbb is not None: # impact bin-by-bin stat impacts_data_stat = tf.sqrt(tf.linalg.diag_part(cov_stat_no_bbb)) impacts_data_stat = _gather_poi_noi_vector( - impacts_data_stat, noiidxs, nsignal_params + impacts_data_stat, noiidxs, nsignal_params, nmodel_params ) impacts_data_stat = tf.reshape(impacts_data_stat, (-1, 1)) impacts_bbb_sq = tf.linalg.diag_part(cov_stat - cov_stat_no_bbb) - impacts_bbb_sq = _gather_poi_noi_vector(impacts_bbb_sq, noiidxs, nsignal_params) + impacts_bbb_sq = _gather_poi_noi_vector( + impacts_bbb_sq, noiidxs, nsignal_params, nmodel_params + ) impacts_bbb = tf.sqrt(tf.nn.relu(impacts_bbb_sq)) # max(0,x) impacts_bbb = tf.reshape(impacts_bbb, (-1, 1)) impacts_grouped = tf.concat([impacts_data_stat, impacts_bbb], axis=1) else: impacts_data_stat = tf.sqrt(tf.linalg.diag_part(cov_stat)) impacts_data_stat = _gather_poi_noi_vector( - impacts_data_stat, noiidxs, nsignal_params + impacts_data_stat, noiidxs, nsignal_params, nmodel_params ) impacts_data_stat = tf.reshape(impacts_data_stat, (-1, 1)) impacts_grouped = impacts_data_stat if len(systgroupidxs): + # systgroupidxs are syst-relative -> shift into full-x by nmodel_params + # and gather from the full covariance (nsignal_params=0 path). impacts_grouped_syst = tf.map_fn( - lambda idxs: _compute_impact_group( - cov, v[:, nsignal_params:], idxs, nsignal_params - ), + lambda idxs: _compute_impact_group(cov, v, nmodel_params + idxs, 0), tf.ragged.constant(systgroupidxs, dtype=tf.int32), fn_output_signature=tf.TensorSpec( shape=(impacts.shape[0],), dtype=tf.float64 @@ -76,4 +94,17 @@ def impacts_parms( impacts_grouped_syst = tf.transpose(impacts_grouped_syst) impacts_grouped = tf.concat([impacts_grouped_syst, impacts_grouped], axis=1) + if param_groupidxs: + # ParamModel parameter groups, already in full-x space (floating-filtered). + # Appended at the END so the grouped-impact axis is + # [syst groups, stat, (binByBinStat), param groups]. + param_cols = [ + tf.reshape( + _compute_impact_group(cov, v, tf.constant(g, dtype=tf.int32), 0), + (-1, 1), + ) + for g in param_groupidxs + ] + impacts_grouped = tf.concat([impacts_grouped, *param_cols], axis=1) + return impacts, impacts_grouped diff --git a/rabbit/workspace.py b/rabbit/workspace.py index a0db99c..08fa34b 100644 --- a/rabbit/workspace.py +++ b/rabbit/workspace.py @@ -13,13 +13,19 @@ ) -def getGroupedImpactsAxes(indata, bin_by_bin_stat=False, per_process=False): +def getGroupedImpactsAxes( + indata, bin_by_bin_stat=False, per_process=False, extra_groups=None +): impact_names = list(indata.systgroups.astype(str)) impact_names.append("stat") if bin_by_bin_stat: impact_names.append("binByBinStat") if per_process: impact_names.extend([f"binByBinStat{p}" for p in indata.procs.astype(str)]) + # ParamModel parameter groups are appended at the end, matching the column + # order produced by traditional_impacts.impacts_parms. + if extra_groups: + impact_names.extend(extra_groups) return hist.axis.StrCategory(impact_names, name="impacts") @@ -59,8 +65,16 @@ def __init__(self, outdir, outname, fitter, postfix=None): parms = list(fitter.parms.astype(str)) self.impact_axis = hist.axis.StrCategory(parms, name="impacts") self.global_impact_axis = hist.axis.StrCategory(systs, name="impacts") + # ParamModel impact groups (e.g. SCETlib NP gamma_nu / F_eff) are only + # computed by the traditional impacts, so extend that axis only. + param_impact_group_names = [ + name for name, _ in fitter._resolved_param_impact_groups() + ] self.grouped_impact_axis = getGroupedImpactsAxes( - fitter.indata, bin_by_bin_stat=fitter.binByBinStat, per_process=False + fitter.indata, + bin_by_bin_stat=fitter.binByBinStat, + per_process=False, + extra_groups=param_impact_group_names, ) self.grouped_global_impact_axis = getGroupedImpactsAxes( fitter.indata, From 2fa3ae52a4f67bad6d8e8bb9a585df7e2beb4244 Mon Sep 17 00:00:00 2001 From: Luca Lavezzo Date: Thu, 11 Jun 2026 09:36:00 -0400 Subject: [PATCH 6/7] format --- bin/rabbit_fit.py | 16 +++++++----- bin/rabbit_print_pulls_and_constraints.py | 4 ++- rabbit/fitter.py | 32 +++++++++++++---------- rabbit/io_tools.py | 6 +---- 4 files changed, 32 insertions(+), 26 deletions(-) diff --git a/bin/rabbit_fit.py b/bin/rabbit_fit.py index fdbafda..8a6f33e 100755 --- a/bin/rabbit_fit.py +++ b/bin/rabbit_fit.py @@ -568,7 +568,9 @@ def fit(args, fitter, ws, dofit=True): if dofit: _t_min = time.perf_counter() cb = fitter.minimize() - logger.debug(f"[timing] fitter.minimize(): {time.perf_counter() - _t_min:.1f} s") + logger.debug( + f"[timing] fitter.minimize(): {time.perf_counter() - _t_min:.1f} s" + ) # force profiling of beta with final parameter values # TODO avoid the extra calculation and jitting if possible since the relevant calculation @@ -577,7 +579,9 @@ def fit(args, fitter, ws, dofit=True): _t_bb = time.perf_counter() logger.info(f"profile beta") fitter._profile_beta() - logger.debug(f"[timing] _profile_beta(): {time.perf_counter() - _t_bb:.1f} s") + logger.debug( + f"[timing] _profile_beta(): {time.perf_counter() - _t_bb:.1f} s" + ) if cb is not None: ws.add_1D_integer_hist(cb.loss_history, "epoch", "loss") @@ -922,10 +926,10 @@ def main(): # can see what was applied without parsing the rabbit log. if getattr(ifitter, "param_prior_active", False): meta["param_priors"] = { - "params": ifitter.param_model.params, # all nparams names - "mask": ifitter.param_prior_mask_np, # bool array - "sigmas": ifitter.param_prior_sigmas_np, # NaN where mask False - "means": ifitter.param_prior_means_np, # NaN where mask False + "params": ifitter.param_model.params, # all nparams names + "mask": ifitter.param_prior_mask_np, # bool array + "sigmas": ifitter.param_prior_sigmas_np, # NaN where mask False + "means": ifitter.param_prior_means_np, # NaN where mask False } with workspace.Workspace( diff --git a/bin/rabbit_print_pulls_and_constraints.py b/bin/rabbit_print_pulls_and_constraints.py index be4430f..e29daee 100755 --- a/bin/rabbit_print_pulls_and_constraints.py +++ b/bin/rabbit_print_pulls_and_constraints.py @@ -143,7 +143,9 @@ def main(): ) ) else: - header = f" {'Parameter':<30} {'pull':>6} +/- {'constraint':>10}{prefit_header}" + header = ( + f" {'Parameter':<30} {'pull':>6} +/- {'constraint':>10}{prefit_header}" + ) print(header) print(" " + "-" * (len(header) - 3)) print( diff --git a/rabbit/fitter.py b/rabbit/fitter.py index 9367954..b16c138 100644 --- a/rabbit/fitter.py +++ b/rabbit/fitter.py @@ -41,7 +41,9 @@ def match_regexp_params(regular_expressions, parameter_names): seen = set() for s in parameter_names: decoded = s.decode() if hasattr(s, "decode") else s - if decoded in exact_lookup or any(r.match(decoded) for r in compiled_expressions): + if decoded in exact_lookup or any( + r.match(decoded) for r in compiled_expressions + ): if decoded not in seen: seen.add(decoded) matched.append(s) @@ -112,8 +114,8 @@ def __init__( # skipped at call site so each method falls back to scipy defaults # for what's not set. self.minimizer_maxiter = getattr(options, "minimizerMaxiter", None) - self.minimizer_gtol = getattr(options, "minimizerGtol", None) - self.minimizer_ftol = getattr(options, "minimizerFtol", None) + self.minimizer_gtol = getattr(options, "minimizerGtol", None) + self.minimizer_ftol = getattr(options, "minimizerFtol", None) self.hvp_method = getattr(options, "hvpMethod", "revrev") # jitCompile accepts "auto" (the default), "on", or "off". # True / False from programmatic callers are accepted as @@ -375,8 +377,10 @@ def init_fit_parms( self.param_prior_sigmas_np = np.where(mask, sigmas, np.nan) self.param_prior_means_np = np.where(mask, means, np.nan) inv2_np = np.where(mask, 1.0 / sigmas**2, 0.0) - self.param_prior_mask = tf.constant(mask, dtype=tf.bool) - self.param_prior_inv2 = tf.constant(inv2_np, dtype=self.indata.dtype) + self.param_prior_mask = tf.constant(mask, dtype=tf.bool) + self.param_prior_inv2 = tf.constant( + inv2_np, dtype=self.indata.dtype + ) self.param_prior_means = tf.constant( np.where(mask, means, 0.0), dtype=self.indata.dtype ) @@ -391,7 +395,9 @@ def init_fit_parms( for i, p in enumerate(pm.params): if mask[i]: name = p.decode() if isinstance(p, bytes) else str(p) - logger.info(f" {name}: μ={means[i]:.4g} σ={sigmas[i]:.4g}") + logger.info( + f" {name}: μ={means[i]:.4g} σ={sigmas[i]:.4g}" + ) else: logger.info( "[paramPriors] --paramModelPriors set but prior_sigmas " @@ -1950,11 +1956,11 @@ def _compute_lc(self, full_nll=False): # Per-entry mask is baked into param_prior_inv2 (zero where mask is # False), so this is a single elementwise op with no branching. if self.param_prior_active: - pm_params = tf.concat( - [self.get_poi(), self.get_model_nui()], axis=0 - ) - lc_pm = 0.5 * self.param_prior_inv2 * tf.square( - pm_params - self.param_prior_means + pm_params = tf.concat([self.get_poi(), self.get_model_nui()], axis=0) + lc_pm = ( + 0.5 + * self.param_prior_inv2 + * tf.square(pm_params - self.param_prior_means) ) if full_nll: # 0.5*log(2π σ²) on entries with a prior. The log(σ²) term is @@ -2214,9 +2220,7 @@ def scipy_hess(xval): sci_opts["gtol"] = float(self.minimizer_gtol) if self.minimizer_ftol is not None: sci_opts["ftol"] = float(self.minimizer_ftol) - logger.info( - f"[minimize] method={self.minimizer_method} options={sci_opts}" - ) + logger.info(f"[minimize] method={self.minimizer_method} options={sci_opts}") try: res = scipy.optimize.minimize( diff --git a/rabbit/io_tools.py b/rabbit/io_tools.py index fabbd2c..b881ed9 100644 --- a/rabbit/io_tools.py +++ b/rabbit/io_tools.py @@ -47,11 +47,7 @@ def read_impacts_poi( ): # read impacts of a single POI - if ( - asym - and impact_type == "traditional" - and "impacts_asym" not in fitresult.keys() - ): + if asym and impact_type == "traditional" and "impacts_asym" not in fitresult.keys(): # Fallback: read asymmetric traditional impacts from a generic # contour scan output if --asymImpacts wasn't run. h_impacts = fitresult["contour_scans"].get()[{"confidence_level": "1.0"}] From 6bfa4ffcbad0e80d58423e1ecbca41654131b2c8 Mon Sep 17 00:00:00 2001 From: Luca Lavezzo Date: Thu, 11 Jun 2026 09:48:48 -0400 Subject: [PATCH 7/7] merge main --- bin/rabbit_fit.py | 2 +- rabbit/impacts/global_asym_impacts.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/bin/rabbit_fit.py b/bin/rabbit_fit.py index 8a6f33e..607990a 100755 --- a/bin/rabbit_fit.py +++ b/bin/rabbit_fit.py @@ -575,7 +575,7 @@ def fit(args, fitter, ws, dofit=True): # force profiling of beta with final parameter values # TODO avoid the extra calculation and jitting if possible since the relevant calculation # usually would have been done during the minimization - if fitter.binByBinStat and not args.noPostfitProfileBB: + if fitter.bbstat.enabled and not args.noPostfitProfileBB: _t_bb = time.perf_counter() logger.info(f"profile beta") fitter._profile_beta() diff --git a/rabbit/impacts/global_asym_impacts.py b/rabbit/impacts/global_asym_impacts.py index f928d3a..a169cc0 100644 --- a/rabbit/impacts/global_asym_impacts.py +++ b/rabbit/impacts/global_asym_impacts.py @@ -152,7 +152,7 @@ def global_asym_impacts_parms( try: fitter.minimize() - if fitter.binByBinStat: + if fitter.bbstat.enabled: fitter._profile_beta() impacts[k, j] = (fitter.x.value() - x_nom).numpy() except Exception as e: @@ -167,7 +167,7 @@ def global_asym_impacts_parms( # Restore the fit state so downstream postfit computations see the nominal. fitter.theta0.assign(theta0_nom) fitter.x.assign(x_nom) - if fitter.binByBinStat: + if fitter.bbstat.enabled: fitter._profile_beta() if n_scanned > 0: