fit ops qol

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

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
                     )

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:

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",