From d596b4a270ad3c19087bfd5c1e34724e3235b920 Mon Sep 17 00:00:00 2001 From: davidwalter2 Date: Tue, 6 Jan 2026 15:16:20 -0500 Subject: [PATCH 1/6] Allow to use externalPostfit as starting point for fit; profile beta parameters after loading externalPostfit --- bin/rabbit_fit.py | 52 ++++++++++++++++++++++++----------------------- 1 file changed, 27 insertions(+), 25 deletions(-) diff --git a/bin/rabbit_fit.py b/bin/rabbit_fit.py index 58d7170..968c396 100755 --- a/bin/rabbit_fit.py +++ b/bin/rabbit_fit.py @@ -276,41 +276,43 @@ def fit(args, fitter, ws, dofit=True): if args.externalPostfit is not None: fitter.load_fitresult(args.externalPostfit, args.externalPostfitResult) - else: - if dofit: - cb = fitter.minimize() - if cb is not None: - ws.add_loss_time_hist(cb.loss_history, cb.time_history) + if fitter.binByBinStat: + fitter._profile_beta() - if not args.noHessian: - # compute the covariance matrix and estimated distance to minimum + if dofit: + cb = fitter.minimize() - val, grad, hess = fitter.loss_val_grad_hess() - edmval, cov = edmval_cov(grad, hess) + if cb is not None: + ws.add_loss_time_hist(cb.loss_history, cb.time_history) - logger.info(f"edmval: {edmval}") + if not args.noHessian: + # compute the covariance matrix and estimated distance to minimum + + val, grad, hess = fitter.loss_val_grad_hess() + edmval, cov = edmval_cov(grad, hess) + logger.info(f"edmval: {edmval}") - fitter.cov.assign(cov) + fitter.cov.assign(cov) - del cov + del cov - if fitter.binByBinStat and fitter.diagnostics: - # This is the estimated distance to minimum with respect to variations of - # the implicit binByBinStat nuisances beta at fixed parameter values. - # It should be near-zero by construction as long as the analytic profiling is - # correct - _, gradbeta, hessbeta = fitter.loss_val_grad_hess_beta() - edmvalbeta, covbeta = edmval_cov(gradbeta, hessbeta) - logger.info(f"edmvalbeta: {edmvalbeta}") + if fitter.binByBinStat and fitter.diagnostics: + # This is the estimated distance to minimum with respect to variations of + # the implicit binByBinStat nuisances beta at fixed parameter values. + # It should be near-zero by construction as long as the analytic profiling is + # correct + _, gradbeta, hessbeta = fitter.loss_val_grad_hess_beta() + edmvalbeta, covbeta = edmval_cov(gradbeta, hessbeta) + logger.info(f"edmvalbeta: {edmvalbeta}") - if args.doImpacts: - ws.add_impacts_hists(*fitter.impacts_parms(hess)) + if args.doImpacts: + ws.add_impacts_hists(*fitter.impacts_parms(hess)) - del hess + del hess - if args.globalImpacts: - ws.add_global_impacts_hists(*fitter.global_impacts_parms()) + if args.globalImpacts: + ws.add_global_impacts_hists(*fitter.global_impacts_parms()) nllvalreduced = fitter.reduced_nll().numpy() From 54a575ce90a4e8909347527910ed25b77f0ca7a6 Mon Sep 17 00:00:00 2001 From: davidwalter2 Date: Sun, 25 Jan 2026 08:59:17 -0500 Subject: [PATCH 2/6] Add option to skip postfit profiling (e.g. only using externalPostfit); Add possibility to only load parameters, not covariance from external postfit; always draw zero line in hist plot --- bin/rabbit_fit.py | 10 ++++++++-- bin/rabbit_plot_hists.py | 11 ++++++++++- rabbit/fitter.py | 15 ++++++++++----- 3 files changed, 28 insertions(+), 8 deletions(-) diff --git a/bin/rabbit_fit.py b/bin/rabbit_fit.py index 968c396..2bb289d 100755 --- a/bin/rabbit_fit.py +++ b/bin/rabbit_fit.py @@ -149,6 +149,12 @@ def make_parser(): type=str, help="Specify result from external postfit file", ) + parser.add_argument( + "--noPostfitProfile", + default=False, + action="store_true", + help="Do not profile beta parameters in the postfit (e.g. when using --externalPostfit).", + ) parser.add_argument( "--doImpacts", default=False, @@ -333,7 +339,7 @@ def fit(args, fitter, ws, dofit=True): "nllvalreduced": nllvalreduced, "ndfsat": ndfsat, "edmval": edmval, - "postfit_profile": args.externalPostfit is None, + "postfit_profile": not args.noPostfitProfile, } ) @@ -561,7 +567,7 @@ def main(): ifitter, ws, prefit=False, - profile=args.externalPostfit is None, + profile=not args.noPostfitProfile, ) else: fit_time.append(time.time()) diff --git a/bin/rabbit_plot_hists.py b/bin/rabbit_plot_hists.py index bb66e0a..aebf7af 100755 --- a/bin/rabbit_plot_hists.py +++ b/bin/rabbit_plot_hists.py @@ -817,6 +817,7 @@ def make_plot( h_den = h_inclusive if diff: + h0 = hh.addHists(h_num, h_num, scale2=-1) h1 = hh.addHists(h_inclusive, h_den, scale2=-1) h2 = hh.addHists(h_data, h_den, scale2=-1) if h_data_stat is not None: @@ -824,6 +825,14 @@ def make_plot( h_data_stat, h_den, cutoff=cutoff, rel_unc=True ) else: + h0 = hh.divideHists( + h_num, + h_num, + cutoff=1e-8, + rel_unc=True, + flow=False, + by_ax_name=False, + ) h1 = hh.divideHists( h_inclusive, h_den, @@ -839,7 +848,7 @@ def make_plot( ) hep.histplot( - h1, + h0, histtype="step", color="grey", alpha=0.5, diff --git a/rabbit/fitter.py b/rabbit/fitter.py index faba62a..e754ca2 100644 --- a/rabbit/fitter.py +++ b/rabbit/fitter.py @@ -288,12 +288,14 @@ def __init__( def load_fitresult(self, fitresult_file, fitresult_key): # load results from external fit and set postfit value and covariance elements for common parameters + cov_ext = None with h5py.File(fitresult_file, "r") as fext: if "x" in fext.keys(): # fitresult from combinetf x_ext = fext["x"][...] parms_ext = fext["parms"][...].astype(str) - cov_ext = fext["cov"][...] + if "cov" in fext.keys(): + cov_ext = fext["cov"][...] else: # fitresult from rabbit h5results_ext = io_tools.get_fitresult(fext, fitresult_key) @@ -301,10 +303,10 @@ def load_fitresult(self, fitresult_file, fitresult_key): x_ext = h_parms_ext.values() parms_ext = np.array(h_parms_ext.axes["parms"]) - cov_ext = h5results_ext["cov"].get().values() + if "cov" in h5results_ext.keys(): + cov_ext = h5results_ext["cov"].get().values() xvals = self.x.numpy() - covval = self.cov.numpy() parms = self.parms.astype(str) # Find common elements with their matching indices @@ -312,10 +314,13 @@ def load_fitresult(self, fitresult_file, fitresult_key): parms, parms_ext, assume_unique=True, return_indices=True ) xvals[idxs] = x_ext[idxs_ext] - covval[np.ix_(idxs, idxs)] = cov_ext[np.ix_(idxs_ext, idxs_ext)] self.x.assign(xvals) - self.cov.assign(tf.constant(covval)) + + if cov_ext is not None: + covval = self.cov.numpy() + covval[np.ix_(idxs, idxs)] = cov_ext[np.ix_(idxs_ext, idxs_ext)] + self.cov.assign(tf.constant(covval)) def update_frozen_params(self): new_mask_np = np.isin(self.parms, self.frozen_params) From c061cbd42acf967ee982283cb0d7dde9b414aaa3 Mon Sep 17 00:00:00 2001 From: davidwalter2 Date: Sun, 25 Jan 2026 11:48:37 -0500 Subject: [PATCH 3/6] Minor fix on pulls labeling --- bin/rabbit_plot_pulls_and_impacts.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/bin/rabbit_plot_pulls_and_impacts.py b/bin/rabbit_plot_pulls_and_impacts.py index d70b50b..e56370d 100755 --- a/bin/rabbit_plot_pulls_and_impacts.py +++ b/bin/rabbit_plot_pulls_and_impacts.py @@ -325,7 +325,7 @@ def make_bar( size=8, ), error_x=error_x, - name="Pulls ± Constraints", + name=f"Pulls ± Constraints ({name})", showlegend=include_ref, ), row=1, From 703627493739dfeb14598ba7b08950bf6c5ae70b Mon Sep 17 00:00:00 2001 From: davidwalter2 Date: Mon, 26 Jan 2026 13:30:32 -0500 Subject: [PATCH 4/6] Make it possible to really disable postfit profiling --- bin/rabbit_fit.py | 9 ++++++--- rabbit/fitter.py | 6 ------ 2 files changed, 6 insertions(+), 9 deletions(-) diff --git a/bin/rabbit_fit.py b/bin/rabbit_fit.py index 2bb289d..0831a5f 100755 --- a/bin/rabbit_fit.py +++ b/bin/rabbit_fit.py @@ -283,12 +283,15 @@ def fit(args, fitter, ws, dofit=True): if args.externalPostfit is not None: fitter.load_fitresult(args.externalPostfit, args.externalPostfitResult) - if fitter.binByBinStat: - fitter._profile_beta() - if dofit: cb = fitter.minimize() + # 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.noPostfitProfile: + fitter._profile_beta() + if cb is not None: ws.add_loss_time_hist(cb.loss_history, cb.time_history) diff --git a/rabbit/fitter.py b/rabbit/fitter.py index e754ca2..bb5e46a 100644 --- a/rabbit/fitter.py +++ b/rabbit/fitter.py @@ -2157,12 +2157,6 @@ def scipy_hess(xval): self.x.assign(xval) - # 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 self.binByBinStat: - self._profile_beta() - return callback def nll_scan(self, param, scan_range, scan_points, use_prefit=False): From 79ff7ae4736ea6286636956c6a07c7f66165366f Mon Sep 17 00:00:00 2001 From: davidwalter2 Date: Mon, 26 Jan 2026 16:37:58 -0500 Subject: [PATCH 5/6] Address Luca and Marco's comments --- bin/rabbit_fit.py | 2 +- bin/rabbit_plot_hists.py | 18 +----------------- bin/rabbit_plot_pulls_and_impacts.py | 19 ++++++++++++------- 3 files changed, 14 insertions(+), 25 deletions(-) diff --git a/bin/rabbit_fit.py b/bin/rabbit_fit.py index 0831a5f..3cfe9bf 100755 --- a/bin/rabbit_fit.py +++ b/bin/rabbit_fit.py @@ -153,7 +153,7 @@ def make_parser(): "--noPostfitProfile", default=False, action="store_true", - help="Do not profile beta parameters in the postfit (e.g. when using --externalPostfit).", + help="Do not profile the bin-by-bin parameters in the postfit (e.g. if parameters are loaded from another fit using --externalPostfit and/or no data is available to be profiled).", ) parser.add_argument( "--doImpacts", diff --git a/bin/rabbit_plot_hists.py b/bin/rabbit_plot_hists.py index aebf7af..9d92e8d 100755 --- a/bin/rabbit_plot_hists.py +++ b/bin/rabbit_plot_hists.py @@ -818,7 +818,6 @@ def make_plot( if diff: h0 = hh.addHists(h_num, h_num, scale2=-1) - h1 = hh.addHists(h_inclusive, h_den, scale2=-1) h2 = hh.addHists(h_data, h_den, scale2=-1) if h_data_stat is not None: h2_stat = hh.divideHists( @@ -833,14 +832,6 @@ def make_plot( flow=False, by_ax_name=False, ) - h1 = hh.divideHists( - h_inclusive, - h_den, - cutoff=1e-8, - rel_unc=True, - flow=False, - by_ax_name=False, - ) h2 = hh.divideHists(h_data, h_den, cutoff=cutoff, rel_unc=True) if h_data_stat is not None: h2_stat = hh.divideHists( @@ -1320,14 +1311,7 @@ def make_plots( h_data_stat = None if hists_up is not None: - hup = [ - ( - h[{k.replace("Sig", ""): v for k, v in idxs.items()}] - if h is not None - else None - ) - for h in hists_up - ] + hup = [h[idxs] if h is not None else None for h in hists_up] else: hup = None diff --git a/bin/rabbit_plot_pulls_and_impacts.py b/bin/rabbit_plot_pulls_and_impacts.py index e56370d..d109aba 100755 --- a/bin/rabbit_plot_pulls_and_impacts.py +++ b/bin/rabbit_plot_pulls_and_impacts.py @@ -206,6 +206,9 @@ def plotImpacts( else: text_on_bars = True + name_suffix = f" ({name})" if name else "" + ref_name_suffix = f" ({ref_name})" if ref_name else "" + if impacts: def make_bar( @@ -234,7 +237,7 @@ def make_bar( ) sign = "+/-" if (oneSidedImpacts and not any(df["impact_down"] > 0)) else "+" - label = f"{sign}1σ impact ({name})" if name else f"{sign}1σ impact" + label = f"{sign}1σ impact{name_suffix}" fig.add_trace( make_bar( key="impact_up", @@ -249,7 +252,7 @@ def make_bar( fig.add_trace( make_bar( key="impact_up_ref", - name=f"{sign}1σ impact ({ref_name})", + name=f"{sign}1σ impact{ref_name_suffix}", filled=False, ), row=1, @@ -260,7 +263,7 @@ def make_bar( fig.add_trace( make_bar( key="impact_down", - name=f"-1σ impact ({name})" if name else "-1σ impact", + name=f"-1σ impact{name_suffix}", color="#e41a1c", text_on_bars=text_on_bars, opacity=0.5 if include_ref else 1, @@ -272,7 +275,7 @@ def make_bar( fig.add_trace( make_bar( key="impact_down_ref", - name=f"-1σ impact ({ref_name})", + name=f"-1σ impact{ref_name_suffix}", color="#e41a1c", filled=False, ), @@ -325,7 +328,7 @@ def make_bar( size=8, ), error_x=error_x, - name=f"Pulls ± Constraints ({name})", + name=f"Pulls ± Constraints{name_suffix}", showlegend=include_ref, ), row=1, @@ -346,7 +349,7 @@ def make_bar( y=labels, orientation="h", **get_marker(filled=False, color="black"), - name=f"Pulls ± Constraints ({ref_name})", + name=f"Pulls ± Constraints{ref_name_suffix}", showlegend=True, ), row=1, @@ -386,7 +389,7 @@ def make_bar( width=1 ), # Adjust the thickness of the marker lines ), - name=f"Asym. pulls ({ref_name})", + name=f"Asym. pulls{ref_name_suffix}", showlegend=include_ref, ), row=1, @@ -760,11 +763,13 @@ def parseArgs(): parser.add_argument( "--refName", type=str, + default="ref.", help="Name of reference input for legend", ) parser.add_argument( "--name", type=str, + default=None, help="Name of input for legend", ) parser.add_argument( From 108a51d9bf74b6b632833ad73cf345dc9dfa2a6d Mon Sep 17 00:00:00 2001 From: davidwalter2 Date: Tue, 27 Jan 2026 10:26:27 -0500 Subject: [PATCH 6/6] Rename option --- bin/rabbit_fit.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/bin/rabbit_fit.py b/bin/rabbit_fit.py index 3cfe9bf..44893f1 100755 --- a/bin/rabbit_fit.py +++ b/bin/rabbit_fit.py @@ -150,7 +150,7 @@ def make_parser(): help="Specify result from external postfit file", ) parser.add_argument( - "--noPostfitProfile", + "--noPostfitProfileBB", default=False, action="store_true", help="Do not profile the bin-by-bin parameters in the postfit (e.g. if parameters are loaded from another fit using --externalPostfit and/or no data is available to be profiled).", @@ -289,7 +289,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.noPostfitProfile: + if fitter.binByBinStat and not args.noPostfitProfileBB: fitter._profile_beta() if cb is not None: @@ -342,7 +342,7 @@ def fit(args, fitter, ws, dofit=True): "nllvalreduced": nllvalreduced, "ndfsat": ndfsat, "edmval": edmval, - "postfit_profile": not args.noPostfitProfile, + "postfit_profile": not args.noPostfitProfileBB, } ) @@ -570,7 +570,7 @@ def main(): ifitter, ws, prefit=False, - profile=not args.noPostfitProfile, + profile=not args.noPostfitProfileBB, ) else: fit_time.append(time.time())