Work on limits

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())
 

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 <mapping> channel axes' for mapping results.",
+        help="Plot impacts on observables use '-m <mapping> 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 <mapping> 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,19 +1253,21 @@ 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)
 
     translate_label = getattr(config, "impact_labels", {})
 
     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]

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)