diff --git a/bin/rabbit_fit.py b/bin/rabbit_fit.py index b5f973d..4ea79cc 100755 --- a/bin/rabbit_fit.py +++ b/bin/rabbit_fit.py @@ -434,24 +434,47 @@ 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") + 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 - - # take covariance from externalPostfit in case the fit was skipped - if dofit or args.externalPostfit is 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 + # 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() @@ -499,13 +522,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 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..a7da4fe 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 @@ -386,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(): @@ -426,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() @@ -1862,6 +1880,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 +1903,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",