From 1f3ec520893d3dcb9d394db9b6d93a610d4a44af Mon Sep 17 00:00:00 2001 From: Luca Lavezzo Date: Thu, 11 Jun 2026 17:50:23 -0400 Subject: [PATCH 1/6] Fit ops QoL: minimizer CLI controls, timing logs, externalPostfit fixes New CLI controls --minimizerMaxiter / --minimizerGtol / --minimizerFtol forwarded to scipy.optimize.minimize as options={...}; unset values keep the historical tol=0.0 behavior. Timing instrumentation at debug verbosity: per-iteration dt/elapsed in FitterCallback and [timing] lines around minimize(), _profile_beta(), and the Hessian-free EDM/CG postfit step. Two --externalPostfit fixes: (1) an Asimov dataset (-t -1) now loads the external postfit values before computing the expected yield, so the Asimov is generated at the postfit point rather than silently at prefit; (2) the postfit covariance is recomputed at the loaded postfit point when a Hessian is wanted (no --noHessian), enabling the two-pass covariance recipe with --externalPostfit --noFit. The --noEDM help text is updated to describe what is actually skipped. rabbit_print_pulls_and_constraints.py gains --noPrefit to suppress the prefit pull/constraint columns; the header separator now matches the actual header width. Co-Authored-By: Claude Fable 5 --- bin/rabbit_fit.py | 47 ++++++++++++++++++++++- bin/rabbit_print_pulls_and_constraints.py | 37 +++++++++++++----- rabbit/fitter.py | 33 +++++++++++++++- rabbit/parsing.py | 28 +++++++++++++- 4 files changed, 131 insertions(+), 14 deletions(-) diff --git a/bin/rabbit_fit.py b/bin/rabbit_fit.py index b5f973d..a7d266f 100755 --- a/bin/rabbit_fit.py +++ b/bin/rabbit_fit.py @@ -434,24 +434,41 @@ 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.bbstat.enabled and not args.noPostfitProfileBB: + _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" + ) 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 - # take covariance from externalPostfit in case the fit was skipped - if dofit or args.externalPostfit is None: + # take covariance from externalPostfit in case the fit was skipped — + # unless a Hessian is explicitly wanted (no --noHessian): the two-pass + # straight-through covariance recipe (--externalPostfit ... --noFit, + # WITHOUT --noHessian) relies on recomputing the Hessian at the loaded + # postfit point. + if dofit or args.externalPostfit is None or not args.noHessian: if not args.noEDM and not args.noHessian: # compute the covariance matrix and estimated distance to minimum _, grad, hess = fitter.loss_val_grad_hess() @@ -499,13 +516,26 @@ 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. + _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 @@ -748,6 +778,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/bin/rabbit_print_pulls_and_constraints.py b/bin/rabbit_print_pulls_and_constraints.py index 2f435b0..e29daee 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,20 @@ def main(): ) ) else: - print( - f" {'Parameter':<30} {'pull':>6} +/- {'constraint':>10} ({'pull prefit':>11} +/- {'constraint prefit':>17})" + header = ( + f" {'Parameter':<30} {'pull':>6} +/- {'constraint':>10}{prefit_header}" ) - print(" " + "-" * 100) + 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 1906542..eb47528 100644 --- a/rabbit/fitter.py +++ b/rabbit/fitter.py @@ -54,7 +54,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}") @@ -68,7 +75,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 @@ -92,6 +99,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 @@ -1862,6 +1875,21 @@ 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, @@ -1870,6 +1898,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/parsing.py b/rabbit/parsing.py index 085ba1e..f8821e7 100644 --- a/rabbit/parsing.py +++ b/rabbit/parsing.py @@ -244,7 +244,33 @@ def common_parser(): "--noEDM", default=False, action="store_true", - help="Don't compute the estimated distance to minimum as fit quality evaluation", + 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( + "--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( "--forceLinear", From 96498b11b14b54bf11920cfb151d7e0e8f6e0a93 Mon Sep 17 00:00:00 2001 From: Luca Lavezzo Date: Thu, 11 Jun 2026 17:57:19 -0400 Subject: [PATCH 2/6] Address review: simplify minimizer-iterations debug log Per review on WMass/rabbit#133: drop the unnecessary try/except and the incorrect L-BFGS label (the callback counts iterations of whichever minimizer method is in use). Co-Authored-By: Claude Fable 5 --- bin/rabbit_fit.py | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/bin/rabbit_fit.py b/bin/rabbit_fit.py index a7d266f..0f537ba 100755 --- a/bin/rabbit_fit.py +++ b/bin/rabbit_fit.py @@ -454,11 +454,9 @@ def fit(args, fitter, ws, dofit=True): 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 + logger.debug( + f"[timing] minimizer: {len(cb.loss_history)} iterations recorded" + ) # prefit variances as the default fallback for add_parms_hist below parms_variances = None From b5c3713fb482e36534e5178cdf731b805ad2240e Mon Sep 17 00:00:00 2001 From: Luca Lavezzo Date: Thu, 11 Jun 2026 18:01:44 -0400 Subject: [PATCH 3/6] Address review: drop forced Hessian recompute under --externalPostfit Per review on WMass/rabbit#133 ('or not args.noHessian' should be removed): with --externalPostfit --noFit the covariance is taken from the loaded fitresult again, as before. The recompute also failed with a non-positive-definite Hessian on the test tensor's Asimov workflow. Co-Authored-By: Claude Fable 5 --- bin/rabbit_fit.py | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/bin/rabbit_fit.py b/bin/rabbit_fit.py index 0f537ba..2894a8b 100755 --- a/bin/rabbit_fit.py +++ b/bin/rabbit_fit.py @@ -461,12 +461,8 @@ def fit(args, fitter, ws, dofit=True): # prefit variances as the default fallback for add_parms_hist below parms_variances = None - # take covariance from externalPostfit in case the fit was skipped — - # unless a Hessian is explicitly wanted (no --noHessian): the two-pass - # straight-through covariance recipe (--externalPostfit ... --noFit, - # WITHOUT --noHessian) relies on recomputing the Hessian at the loaded - # postfit point. - if dofit or args.externalPostfit is None or not args.noHessian: + # take covariance from externalPostfit in case the fit was skipped + if dofit or args.externalPostfit is None: if not args.noEDM and not args.noHessian: # compute the covariance matrix and estimated distance to minimum _, grad, hess = fitter.loss_val_grad_hess() From 9a5767fdef8fdf5ee3c56383be63a4a5aa55e070 Mon Sep 17 00:00:00 2001 From: Luca Lavezzo Date: Fri, 12 Jun 2026 07:49:20 -0400 Subject: [PATCH 4/6] externalPostfit: recompute Hessian when the fitresult has no covariance Reconciles the review discussion on WMass/rabbit#133 (rabbit_fit.py:603): when --externalPostfit provides a covariance it is used as before; when it does not (the fitresult was produced with --noHessian), the Hessian is recomputed at the loaded postfit point instead of silently writing no covariance. This enables the two-pass recipe for models whose full Hessian is infeasible during the fit: pass 1 fits with --noHessian, pass 2 reruns with --externalPostfit ... --noFit (without --noHessian) to obtain the covariance. Verified on the test tensor that the two-pass covariance is identical to the single-pass one, and that a fitresult with a covariance is still consumed without recomputation. Co-Authored-By: Claude Fable 5 --- bin/rabbit_fit.py | 13 +++++++++++-- rabbit/fitter.py | 5 +++++ 2 files changed, 16 insertions(+), 2 deletions(-) diff --git a/bin/rabbit_fit.py b/bin/rabbit_fit.py index 2894a8b..a280bad 100755 --- a/bin/rabbit_fit.py +++ b/bin/rabbit_fit.py @@ -461,8 +461,17 @@ def fit(args, fitter, ws, dofit=True): # prefit variances as the default fallback for add_parms_hist below parms_variances = None - # take covariance from externalPostfit in case the fit was skipped - if dofit or args.externalPostfit is None: + # take covariance from externalPostfit in case the fit was skipped. + # If the external fitresult does not contain a covariance (e.g. it was + # produced with --noHessian), recompute the Hessian at the loaded postfit + # point instead — the two-pass recipe: fit once with --noHessian, then + # rerun with --externalPostfit ... --noFit (without --noHessian) to get + # the covariance. + if ( + dofit + or args.externalPostfit is None + or not getattr(fitter, "external_cov_loaded", False) + ): if not args.noEDM and not args.noHessian: # compute the covariance matrix and estimated distance to minimum _, grad, hess = fitter.loss_val_grad_hess() diff --git a/rabbit/fitter.py b/rabbit/fitter.py index eb47528..a7da4fe 100644 --- a/rabbit/fitter.py +++ b/rabbit/fitter.py @@ -399,6 +399,10 @@ def __deepcopy__(self, memo): def load_fitresult(self, fitresult_file, fitresult_key, profile=True): # load results from external fit and set postfit value and covariance elements for common parameters + # external_cov_loaded records whether the loaded fitresult carried a + # covariance (it does not when it was produced with --noHessian), so + # callers can decide to recompute the Hessian at the loaded point. + self.external_cov_loaded = False cov_ext = None with h5py.File(fitresult_file, "r") as fext: if "x" in fext.keys(): @@ -439,6 +443,7 @@ def load_fitresult(self, fitresult_file, fitresult_key, profile=True): covval = self.cov.numpy() covval[np.ix_(idxs, idxs)] = cov_ext[np.ix_(idxs_ext, idxs_ext)] self.cov.assign(tf.constant(covval)) + self.external_cov_loaded = True if profile: self._profile_beta() From b20990064e8a6395a7d8979bf358b467cf73a4f2 Mon Sep 17 00:00:00 2001 From: Luca Lavezzo Date: Fri, 12 Jun 2026 11:16:56 -0400 Subject: [PATCH 5/6] parms hist: always Weight storage, NaN variances when not computed parms_variances now defaults to a NaN-filled vector instead of None, so the parms hist is always written with Weight storage and entries whose uncertainty was not computed (e.g. --noHessian with --noEDM) read NaN from .variances() downstream, instead of a silently absent or plausible-looking value. Co-Authored-By: Claude Fable 5 --- bin/rabbit_fit.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/bin/rabbit_fit.py b/bin/rabbit_fit.py index a280bad..64e6f99 100755 --- a/bin/rabbit_fit.py +++ b/bin/rabbit_fit.py @@ -458,8 +458,11 @@ def fit(args, fitter, ws, dofit=True): f"[timing] minimizer: {len(cb.loss_history)} iterations recorded" ) - # prefit variances as the default fallback for add_parms_hist below - parms_variances = None + # default for add_parms_hist below: NaN variances so the parms hist is + # always written with Weight storage, and entries whose uncertainty was + # not computed (e.g. --noHessian with --noEDM) read NaN downstream + # instead of a silently absent or plausible-looking value. + parms_variances = np.full(len(fitter.parms), np.nan) # take covariance from externalPostfit in case the fit was skipped. # If the external fitresult does not contain a covariance (e.g. it was From 19b5b5f80a6464abc6b565a121b1501bc3860740 Mon Sep 17 00:00:00 2001 From: Luca Lavezzo Date: Mon, 22 Jun 2026 17:35:57 -0400 Subject: [PATCH 6/6] Revert externalPostfit Asimov-at-postfit change Build the Asimov dataset at the prefit point regardless of --externalPostfit, as before. The change had made `-t -1 --externalPostfit X` generate the Asimov at the loaded postfit point, but that breaks the intended use cases: - Two-pass covariance recipe (minimize, then Hessian via --externalPostfit ... --noFit): step 2 must rebuild the same likelihood step 1 minimized. A one-step Asimov fit holds nobs fixed at the prefit expectation while the parameters minimize, so a faithful split must use the prefit Asimov in step 2. Regenerating it at the postfit point yields a different likelihood whose minimum is no longer the loaded point (EDM != 0 once regularization/priors move the postfit off prefit), giving the wrong covariance. - Blinding operates on the real-data fit (-t 0); Asimov fits are never blinded, so this path is irrelevant there. --externalPostfit still loads the postfit for the postfit step via the in-fit() load. A genuine postfit-Asimov (expected sensitivity at the data-driven point) is a separate use case that should be an explicit opt-in flag, not silent default behavior. Co-Authored-By: Claude Opus 4.8 (1M context) --- bin/rabbit_fit.py | 13 ------------- 1 file changed, 13 deletions(-) diff --git a/bin/rabbit_fit.py b/bin/rabbit_fit.py index 64e6f99..4ea79cc 100755 --- a/bin/rabbit_fit.py +++ b/bin/rabbit_fit.py @@ -784,19 +784,6 @@ 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())