From a614824386508f7e1f836e960fc074322f2f0caa Mon Sep 17 00:00:00 2001 From: davidwalter2 Date: Fri, 9 Jan 2026 13:04:17 -0500 Subject: [PATCH 1/5] Some fixes to limits stuff --- bin/rabbit_limit.py | 116 ++++++++++++++++++++++++------------------- tests/plot_limits.py | 51 ++++++++++++++----- 2 files changed, 104 insertions(+), 63 deletions(-) diff --git a/bin/rabbit_limit.py b/bin/rabbit_limit.py index 5e77a52..5ae71e0 100755 --- a/bin/rabbit_limit.py +++ b/bin/rabbit_limit.py @@ -24,6 +24,14 @@ def make_parser(): parser = parsing.common_parser() parser.add_argument("--outname", default="limits.hdf5", help="output file name") + parser.add_argument( + "--modes", + nargs="+", + default=["gaussian", "likelihood"], + type=str, + choices=["gaussian", "likelihood"], + help="Different modes of approximations to be used", + ) parser.add_argument( "--asymptoticLimits", nargs="+", @@ -42,10 +50,8 @@ def make_parser(): return parser.parse_args() -def do_asymptotic_limits( - args, fitter, ws, data_values, mapping=None, bkg_only_fit=False -): - if bkg_only_fit: +def do_asymptotic_limits(args, fitter, ws, data_values, mapping=None, fit_data=False): + if fit_data: logger.info("Perform background only fit") # set process to zero and freeze fitter.freeze_params(args.asymptoticLimits) @@ -73,9 +79,10 @@ def do_asymptotic_limits( # Clone fitter to simultaneously minimize on asimov dataset fitter_asimov = copy.deepcopy(fitter) - # unconditional fit to real data - fitter.set_nobs(data_values) - fitter.minimize() + if fit_data: + logger.info("Unconditional fit to real data") + fitter.set_nobs(data_values) + fitter.minimize() # asymptotic limits (CLs) # see: @@ -130,13 +137,15 @@ def do_asymptotic_limits( idx = np.where(fitter.parms.astype(str) == key)[0][0] if key in fitter.poi_model.pois.astype(str): + is_poi = True xbest = fitter_asimov.get_blinded_poi()[idx] xobs = fitter.get_blinded_poi()[idx] elif key in fitter.parms.astype(str): + is_poi = False xbest = fitter_asimov.get_blinded_theta()[ - fitter_asimov.poi_model.npoi + idx + idx - fitter_asimov.poi_model.npoi ] - xobs = fitter.get_blinded_theta()[fitter_asimov.poi_model.npoi + idx] + xobs = fitter.get_blinded_theta()[idx - fitter.poi_model.npoi] xerr = fitter_asimov.cov[idx, idx] ** 0.5 @@ -149,7 +158,7 @@ def do_asymptotic_limits( fitter.cov.assign(cov) xobs_err = fitter.cov[idx, idx] ** 0.5 - if not args.allowNegativePOI: + if is_poi and not args.allowNegativePOI: xerr = 2 * xerr * xbest**0.5 xobs_err = 2 * xobs_err * xobs**0.5 @@ -167,53 +176,58 @@ def do_asymptotic_limits( qmuA = (norm.ppf(cl_b) - norm.ppf(cl_sb)) ** 2 logger.debug(f"Find r with q_(r,A)=-2ln(L)/ln(L0) = {qmuA}") + if "gaussian" in args.modes: + # Gaussian approximation + r = xbest + xerr * qmuA**0.5 + logger.info( + f"Expected (Gaussian) {round((cl_b)*100,1):4.1f}%: {key} < {r}" + ) + limits[i, j, k] = r + + if "likelihood" in args.modes: + # Likelihood approximation + r, _ = fitter_asimov.contour_scan( + key, nllvalreduced_asimov, qmuA, signs=[1], fun=fun + ) + logger.info( + f"Expected (Likelihood) {round((cl_b)*100,1):4.1f}%: {key} < {r}" + ) + limits_nll[i, j, k] = r + + if "gaussian" in args.modes: # Gaussian approximation - r = xbest + xerr * qmuA**0.5 - logger.info( - f"Expected (Gaussian) {round((cl_b)*100,1):4.1f}%: {key} < {r}" + limits_obs[i, j] = asymptotic_limits.compute_gaussian_limit( + key, xobs, xobs_err, xerr, cl_s ) - limits[i, j, k] = r + if "likelihood" in args.modes: # Likelihood approximation - r, _ = fitter_asimov.contour_scan( - key, nllvalreduced_asimov, qmuA, signs=[1], fun=fun - ) - logger.info( - f"Expected (Likelihood) {round((cl_b)*100,1):4.1f}%: {key} < {r}" - ) - limits_nll[i, j, k] = r - - # Gaussian approximation - limits_obs[i, j] = asymptotic_limits.compute_gaussian_limit( - key, xobs, xobs_err, xerr, cl_s - ) + if mapping is not None: + # TODO: implement for channels + pass + else: + # TODO: make it work + # nllvalreduced = fitter.reduced_nll().numpy() + # limits_nll_obs[i, j] = asymptotic_limits.compute_likelihood_limit(fitter, fitter_asimov, nllvalreduced, nllvalreduced_asimov, key, cl_s) + pass - # Likelihood approximation - if mapping is not None: - # TODO: implement for channels - pass - else: - # TODO: make it work - # nllvalreduced = fitter.reduced_nll().numpy() - # limits_nll_obs[i, j] = asymptotic_limits.compute_likelihood_limit(fitter, fitter_asimov, nllvalreduced, nllvalreduced_asimov, key, cl_s) - pass - - ws.add_limits_hist( - limits, - args.asymptoticLimits, - args.levels, - clb_list, - base_name="gaussian_asymptoticLimits_expected", - ) + if "gaussian" in args.modes: + ws.add_limits_hist( + limits, + args.asymptoticLimits, + args.levels, + clb_list, + base_name="gaussian_asymptoticLimits_expected", + ) - ws.add_limits_hist( - limits_obs, - args.asymptoticLimits, - args.levels, - base_name="gaussian_asymptoticLimits_observed", - ) + ws.add_limits_hist( + limits_obs, + args.asymptoticLimits, + args.levels, + base_name="gaussian_asymptoticLimits_observed", + ) - if mapping is not None: + if "likelihood" in args.modes: # TODO: implement for channels ws.add_limits_hist( limits_nll, @@ -348,7 +362,7 @@ def main(): ws, data_values=observations, mapping=mapping, - bkg_only_fit=ifit >= 0, + fit_data=ifit >= 0, ) fit_time.append(time.time()) diff --git a/tests/plot_limits.py b/tests/plot_limits.py index cffff1d..2ea9358 100644 --- a/tests/plot_limits.py +++ b/tests/plot_limits.py @@ -5,6 +5,7 @@ import mplhep as hep import numpy as np +from wums import boostHistHelpers as hh from wums import logging from rabbit import io_tools @@ -59,6 +60,13 @@ def parseArgs(): parser.add_argument( "--scale", type=float, default=1.0, help="Scale limits by this value" ) + parser.add_argument( + "--mode", + type=str, + default="gaussian", + choices=["gaussian", "likelihood"], + help="Different modes of approximations to be used", + ) parser.add_argument( "--config", type=str, @@ -127,30 +135,49 @@ def main(): logger = logging.setup_logger(__file__, args.verbose, args.noColorLogger) + key_expected = f"{args.mode}_asymptoticLimits_expected" + key_observed = f"{args.mode}_asymptoticLimits_observed" + cls_list = set() clb_list = set() limits = [] limits_asimov = [] lumi = None + xvals = [] for infile in args.infiles: fitresult_asimov, meta = io_tools.get_fitresult( infile, result="asimov", meta=True ) - h_limits_asimov = fitresult_asimov["asymptoticLimits"].get() + bsm_sample, mixing = meta["meta_info_input"]["meta_info"]["args"][ + "addBSMMixing" + ] + bsm_name, mass = bsm_sample.split("_") + xvals.append(int(mass)) + + scale = float(mixing) * args.scale + + if key_expected in fitresult_asimov.keys(): + h_limits_asimov = fitresult_asimov[key_expected].get() + h_limits_asimov = hh.scaleHist(h_limits_asimov, scale) + else: + raise RuntimeError( + f"Key {key_expected} not found, available keyes are {[key for key in fitresult_asimov.keys()]}" + ) cls_list = cls_list.union(set([cls for cls in h_limits_asimov.axes["cls"]])) clb_list = clb_list.union(set([clb for clb in h_limits_asimov.axes["clb"]])) limits_asimov.append(h_limits_asimov) - fitresult, meta = io_tools.get_fitresult(infile, meta=True) + fitresult = io_tools.get_fitresult(infile) if fitresult != fitresult_asimov and not args.noObs: if ( - "asymptoticLimits" in fitresult.keys() - and "clb" not in fitresult["asymptoticLimits"].get().axes.name + key_observed in fitresult.keys() + and "clb" not in fitresult[key_observed].get().axes.name ): - h_limits = fitresult["asymptoticLimits"].get() + h_limits = fitresult[key_observed].get() + h_limits = hh.scaleHist(h_limits, scale) cls_list = cls_list.union(set([cls for cls in h_limits.axes["cls"]])) limits.append(h_limits) else: @@ -165,33 +192,33 @@ def main(): if args.xvals is not None: x = np.array(args.xvals) else: - x = np.arange(len(limits_asimov)) + x = np.array(xvals) clb_list = list(clb_list) for cls in cls_list: yexp = np.array( - [l[{"cls": cls, "clb": "0.5"}] * args.scale for l in limits_asimov] + [l[{"cls": cls, "clb": "0.5"}] for l in limits_asimov] ).flatten() yexp_m2 = np.array( - [l[{"cls": cls, "clb": "0.025"}] * args.scale for l in limits_asimov] + [l[{"cls": cls, "clb": "0.025"}] for l in limits_asimov] ).flatten() yexp_m1 = np.array( - [l[{"cls": cls, "clb": "0.16"}] * args.scale for l in limits_asimov] + [l[{"cls": cls, "clb": "0.16"}] for l in limits_asimov] ).flatten() yexp_p1 = np.array( - [l[{"cls": cls, "clb": "0.84"}] * args.scale for l in limits_asimov] + [l[{"cls": cls, "clb": "0.84"}] for l in limits_asimov] ).flatten() yexp_p2 = np.array( - [l[{"cls": cls, "clb": "0.975"}] * args.scale for l in limits_asimov] + [l[{"cls": cls, "clb": "0.975"}] for l in limits_asimov] ).flatten() ylist = [yexp, yexp_m1, yexp_m2, yexp_p1, yexp_p2] if len(limits) > 0: - yobs = np.array([l[{"cls": cls}] * args.scale for l in limits]).flatten() + yobs = np.array([l[{"cls": cls}] for l in limits]).flatten() ylist.append(yobs) ylist = np.array(ylist) From 512940f5cf53b0807cae2619da4f1c1bc445158e Mon Sep 17 00:00:00 2001 From: davidwalter2 Date: Fri, 9 Jan 2026 13:05:22 -0500 Subject: [PATCH 2/5] Improve impacts plots to be able to specify channel for observable --- bin/rabbit_plot_pulls_and_impacts.py | 61 ++++++++++++++++++++++------ 1 file changed, 48 insertions(+), 13 deletions(-) diff --git a/bin/rabbit_plot_pulls_and_impacts.py b/bin/rabbit_plot_pulls_and_impacts.py index b847c30..c67660d 100755 --- a/bin/rabbit_plot_pulls_and_impacts.py +++ b/bin/rabbit_plot_pulls_and_impacts.py @@ -14,12 +14,14 @@ from rabbit import io_tools -from wums import output_tools, plot_tools # isort: skip +from wums import logging, output_tools, plot_tools # isort: skip # prevent MathJax from bein loaded pio.kaleido.scope.mathjax = None +logger = None + def writeOutput(fig, outfile, extensions=[], postfix=None, args=None, meta_info=None): name, _ = os.path.splitext(outfile) @@ -700,6 +702,17 @@ def parseArgs(): type=str, help="fitresults output hdf5 file from fit", ) + parser.add_argument( + "-v", + "--verbose", + type=int, + default=3, + choices=[0, 1, 2, 3, 4], + help="Set verbosity level with logging, the larger the more verbose", + ) + parser.add_argument( + "--noColorLogger", action="store_true", help="Do not use logging with colors" + ) parser.add_argument( "--result", default=None, @@ -712,14 +725,26 @@ def parseArgs(): default=None, type=str, nargs="+", - help="Print impacts on observables use '-m channel axes' for mapping results.", + help="Plot impacts on observables use '-m channel axes' for mapping results.", + ) + parser.add_argument( + "--channel", + default=None, + type=str, + help="Plot impacts for given channel", ) parser.add_argument( - "--mappingRef", + "--refMapping", default=None, type=str, nargs="+", - help="Print impacts on observables use '-m channel axes' for mapping results for reference.", + help="Plot impacts on observables from mapping results for reference.", + ) + parser.add_argument( + "--refChannel", + default=None, + type=str, + help="Plot impacts for given channel", ) parser.add_argument( "-r", @@ -1228,6 +1253,8 @@ def produce_plots_hist( def main(): args = parseArgs() + global logger + logger = logging.setup_logger(__file__, args.verbose, args.noColorLogger) config = plot_tools.load_config(args.config) @@ -1235,12 +1262,12 @@ def main(): fitresult, meta = io_tools.get_fitresult(args.inputFile, args.result, meta=True) if any( - x is not None for x in [args.referenceFile, args.refResult, args.mappingRef] + x is not None for x in [args.referenceFile, args.refResult, args.refMapping] ): - referenceFile = ( - args.referenceFile if args.referenceFile is not None else args.inputFile - ) - fitresult_ref = io_tools.get_fitresult(referenceFile, args.refResult) + if args.referenceFile is not None: + fitresult_ref = io_tools.get_fitresult(args.referenceFile, args.refResult) + else: + fitresult_ref = fitresult else: fitresult_ref = None @@ -1292,13 +1319,17 @@ def get_mapping_key(result, key): channels = res[key]["channels"] return channels, key else: - keys = [key for key in res.keys() if key.startswith(key)] + keys = [k for k in res.keys() if k.startswith(key)] if len(keys) == 0: raise ValueError( f"Mapping {key} not found, available mappings are: {res.keys()}" ) channels = res[keys[0]]["channels"] + logger.info( + f"Found mapping {keys[0]} with channels {[k for k in channels.keys()]}" + ) + return channels, keys[0] mapping_key = " ".join(args.mapping) @@ -1307,8 +1338,8 @@ def get_mapping_key(result, key): if fitresult_ref: mapping_key_ref = ( - " ".join(args.mappingRef) - if args.mappingRef is not None + " ".join(args.refMapping) + if args.refMapping is not None else mapping_key ) @@ -1317,6 +1348,8 @@ def get_mapping_key(result, key): ) for channel, hists in channels.items(): + if args.channel and channel not in args.channel: + continue modes = ["ungrouped", "group"] if args.mode == "both" else [args.mode] for mode in modes: @@ -1331,7 +1364,9 @@ def get_mapping_key(result, key): hist = hists[key].get() if fitresult_ref: - if channel in channels_ref.keys(): + if args.refChannel: + channel_ref = args.refChannel + elif channel in channels_ref.keys(): channel_ref = channel elif len(channels_ref.keys()) == 1: channel_ref = [v for v in channels_ref.keys()][0] From e675c92065ea402fef16e38d494c718f4c9c92b9 Mon Sep 17 00:00:00 2001 From: davidwalter2 Date: Fri, 9 Jan 2026 13:05:57 -0500 Subject: [PATCH 3/5] Fix poiModel Ones --- rabbit/poi_models/poi_model.py | 1 + 1 file changed, 1 insertion(+) diff --git a/rabbit/poi_models/poi_model.py b/rabbit/poi_models/poi_model.py index f45d997..758e8b8 100644 --- a/rabbit/poi_models/poi_model.py +++ b/rabbit/poi_models/poi_model.py @@ -64,6 +64,7 @@ def __init__(self, indata, **kwargs): self.poidefault = tf.zeros([], dtype=self.indata.dtype) self.allowNegativePOI = False + self.is_linear = True def compute(self, poi): rnorm = tf.ones(self.indata.nproc, dtype=self.indata.dtype) From ccd88407b003a704e28aa5adf2194eb3f9a1b343 Mon Sep 17 00:00:00 2001 From: davidwalter2 Date: Fri, 9 Jan 2026 13:06:46 -0500 Subject: [PATCH 4/5] Use common fitting functions in contour scan --- rabbit/fitter.py | 52 ++++++++++++++++++++++++++++-------------------- 1 file changed, 30 insertions(+), 22 deletions(-) diff --git a/rabbit/fitter.py b/rabbit/fitter.py index faba62a..5d40ad6 100644 --- a/rabbit/fitter.py +++ b/rabbit/fitter.py @@ -2367,6 +2367,9 @@ def scipy_hess(x, v): xdn = xval[idx] - (self.cov[idx, idx] * q) ** 0.5 for i, sign in enumerate(signs): + logger.info( + f"Try to find contour on the {'positive' if sign else 'negative'} side" + ) # Objective function and its derivatives if sign == -1: xval_init[idx] = xdn @@ -2380,7 +2383,7 @@ def scipy_hess(x, v): self.minimize() self.defreeze_params(param) - opt = {} + info_minimize = {} if fun is None: # contour scan on parameter def objective_val_grad(x): @@ -2397,7 +2400,7 @@ def objective_val_grad(x): n_params = len(xval_init) obj_hess = csr_matrix((n_params, n_params)) - opt["hess"] = lambda x: obj_hess + info_minimize["hess"] = lambda x: obj_hess else: # contour scan on observable def objective_val_grad(x): @@ -2432,30 +2435,35 @@ def objective_hessp(x, pval): hessp = t2.gradient(grad, self.x, output_gradients=p) return hessp.__array__() - opt["hessp"] = objective_hessp + info_minimize["hessp"] = objective_hessp - res = scipy.optimize.minimize( - objective_val_grad, - xval_init, - method="trust-constr", - jac=True, - constraints=[nlc], - options={ - "maxiter": 50000, - "xtol": 1e-10, - "gtol": 1e-10, - # "barrier_tol": 1e-10, - }, - **opt, - ) + callback = FitterCallback(self) + + try: + res = scipy.optimize.minimize( + objective_val_grad, + xval_init, + method="trust-constr", + jac=True, + tol=0.0, + constraints=[nlc], + callback=callback, + **info_minimize, + ) + except Exception as ex: + # minimizer could have called the loss or hessp functions with "random" values, so restore the + # state from the end of the last iteration before the exception + xval = callback.xval + logger.debug(ex) + else: + xval = res["x"] + logger.debug(res) logger.info(f"Success: {res.success}") logger.debug(f"Status: {res.status}") - if not res.success: - logger.warning(f"Message: {res.message}") - logger.warning(f"Optimality (gtol): {res.optimality}") - logger.warning(f"Constraint Violation: {res.constr_violation}") - continue + logger.debug(f"Message: {res.message}") + logger.debug(f"Optimality (gtol): {res.optimality}") + logger.debug(f"Constraint Violation: {res.constr_violation}") params_values[i] = res["x"] - xval From 0ae2b92c073cb98ab61394cbd7a89dc6efaf47e9 Mon Sep 17 00:00:00 2001 From: davidwalter2 Date: Fri, 9 Jan 2026 13:34:57 -0500 Subject: [PATCH 5/5] Revert "Use common fitting functions in contour scan" This reverts commit ccd88407b003a704e28aa5adf2194eb3f9a1b343. --- rabbit/fitter.py | 52 ++++++++++++++++++++---------------------------- 1 file changed, 22 insertions(+), 30 deletions(-) diff --git a/rabbit/fitter.py b/rabbit/fitter.py index 5d40ad6..faba62a 100644 --- a/rabbit/fitter.py +++ b/rabbit/fitter.py @@ -2367,9 +2367,6 @@ def scipy_hess(x, v): xdn = xval[idx] - (self.cov[idx, idx] * q) ** 0.5 for i, sign in enumerate(signs): - logger.info( - f"Try to find contour on the {'positive' if sign else 'negative'} side" - ) # Objective function and its derivatives if sign == -1: xval_init[idx] = xdn @@ -2383,7 +2380,7 @@ def scipy_hess(x, v): self.minimize() self.defreeze_params(param) - info_minimize = {} + opt = {} if fun is None: # contour scan on parameter def objective_val_grad(x): @@ -2400,7 +2397,7 @@ def objective_val_grad(x): n_params = len(xval_init) obj_hess = csr_matrix((n_params, n_params)) - info_minimize["hess"] = lambda x: obj_hess + opt["hess"] = lambda x: obj_hess else: # contour scan on observable def objective_val_grad(x): @@ -2435,35 +2432,30 @@ def objective_hessp(x, pval): hessp = t2.gradient(grad, self.x, output_gradients=p) return hessp.__array__() - info_minimize["hessp"] = objective_hessp + opt["hessp"] = objective_hessp - callback = FitterCallback(self) - - try: - res = scipy.optimize.minimize( - objective_val_grad, - xval_init, - method="trust-constr", - jac=True, - tol=0.0, - constraints=[nlc], - callback=callback, - **info_minimize, - ) - except Exception as ex: - # minimizer could have called the loss or hessp functions with "random" values, so restore the - # state from the end of the last iteration before the exception - xval = callback.xval - logger.debug(ex) - else: - xval = res["x"] - logger.debug(res) + res = scipy.optimize.minimize( + objective_val_grad, + xval_init, + method="trust-constr", + jac=True, + constraints=[nlc], + options={ + "maxiter": 50000, + "xtol": 1e-10, + "gtol": 1e-10, + # "barrier_tol": 1e-10, + }, + **opt, + ) logger.info(f"Success: {res.success}") logger.debug(f"Status: {res.status}") - logger.debug(f"Message: {res.message}") - logger.debug(f"Optimality (gtol): {res.optimality}") - logger.debug(f"Constraint Violation: {res.constr_violation}") + if not res.success: + logger.warning(f"Message: {res.message}") + logger.warning(f"Optimality (gtol): {res.optimality}") + logger.warning(f"Constraint Violation: {res.constr_violation}") + continue params_values[i] = res["x"] - xval