Skip to content

Update external postfit functionality #102

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
davidwalter2 wants to merge 5 commits into
base: main
Choose a base branch
from
+66 −58
Open

Update external postfit functionality #102

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
d596b4a
Allow to use externalPostfit as starting point for fit; profile beta …
davidwalter2 Jan 6, 2026
54a575c
Add option to skip postfit profiling (e.g. only using externalPostfit…
davidwalter2 Jan 25, 2026
c061cbd
Minor fix on pulls labeling
davidwalter2 Jan 25, 2026
7036274
Make it possible to really disable postfit profiling
davidwalter2 Jan 26, 2026
79ff7ae
Address Luca and Marco's comments
davidwalter2 Jan 26, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
65 changes: 38 additions & 27 deletions bin/rabbit_fit.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 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).",
)
Copy link
Collaborator

@lucalavezzo lucalavezzo Jan 26, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It may be helpful to be more descriptive in this help section (or perhaps start a documentation?), because honestly I don't always follow the intricacies of these settings; in particular, here, I'm not sure what the use case for this is, and even a short sentence to give more context may help the layman user :)

Copy link
Collaborator Author

@davidwalter2 davidwalter2 Jan 26, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've modified the help message. In general I agree that we should come up with a better documentation.

parser.add_argument(
"--doImpacts",
default=False,
Expand Down Expand Up @@ -276,41 +282,46 @@ 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 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 not args.noHessian:
# compute the covariance matrix and estimated distance to minimum
if cb is not None:
ws.add_loss_time_hist(cb.loss_history, cb.time_history)

val, grad, hess = fitter.loss_val_grad_hess()
edmval, cov = edmval_cov(grad, hess)
if not args.noHessian:
# compute the covariance matrix and estimated distance to minimum

logger.info(f"edmval: {edmval}")
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()

Expand All @@ -331,7 +342,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,
}
)

Expand Down Expand Up @@ -559,7 +570,7 @@ def main():
ifitter,
ws,
prefit=False,
profile=args.externalPostfit is None,
profile=not args.noPostfitProfile,
)
else:
fit_time.append(time.time())
Expand Down
19 changes: 6 additions & 13 deletions bin/rabbit_plot_hists.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -817,16 +817,16 @@ def make_plot(
h_den = h_inclusive

if diff:
h1 = hh.addHists(h_inclusive, h_den, scale2=-1)
h0 = hh.addHists(h_num, h_num, scale2=-1)
Copy link
Collaborator

@lucalavezzo lucalavezzo Jan 26, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I guess you're using this to get an empty histogram with the same size as h_num, does this treat the variances correctly (if that matters)? Depending on the storage type, it may be summing the variances to find the final ones on h0.

Copy link
Collaborator Author

@davidwalter2 davidwalter2 Jan 26, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is actually only used to draw the grey line around 1 (or 0 in case the difference is plot), without errors.

h2 = hh.addHists(h_data, h_den, scale2=-1)
if h_data_stat is not None:
h2_stat = hh.divideHists(
h_data_stat, h_den, cutoff=cutoff, rel_unc=True
)
else:
h1 = hh.divideHists(
h_inclusive,
h_den,
h0 = hh.divideHists(
h_num,
h_num,
cutoff=1e-8,
rel_unc=True,
flow=False,
Expand All @@ -839,7 +839,7 @@ def make_plot(
)

hep.histplot(
h1,
h0,
histtype="step",
color="grey",
alpha=0.5,
Expand Down Expand Up @@ -1311,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

Expand Down
19 changes: 12 additions & 7 deletions bin/rabbit_plot_pulls_and_impacts.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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(
Expand Down Expand Up @@ -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",
Expand All @@ -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,
Expand All @@ -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,
Expand All @@ -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,
),
Expand Down Expand Up @@ -325,7 +328,7 @@ def make_bar(
size=8,
),
error_x=error_x,
name="Pulls ± Constraints",
name=f"Pulls ± Constraints{name_suffix}",
showlegend=include_ref,
),
row=1,
Expand All @@ -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,
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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(
Expand Down
21 changes: 10 additions & 11 deletions rabbit/fitter.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -288,34 +288,39 @@ 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)
h_parms_ext = h5results_ext["parms"].get()

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
common_elements, idxs, idxs_ext = np.intersect1d(
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)
Expand Down Expand Up @@ -2152,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):
Expand Down