Work on impacts

.github/workflows/main.yml

@@ -348,10 +348,10 @@ jobs:
         run: rabbit_print_impacts.py $RABBIT_OUTDIR/fitresults.hdf5 -s
 
       - name: print global impacts
-        run: rabbit_print_impacts.py $RABBIT_OUTDIR/fitresults.hdf5 -s --globalImpacts
+        run: rabbit_print_impacts.py $RABBIT_OUTDIR/fitresults.hdf5 -s --impactType global
 
       - name: print global impacts of observable
-        run: rabbit_print_impacts.py $RABBIT_OUTDIR/fitresults.hdf5 --globalImpacts -m Project ch0_masked --ungroup --sort
+        run: rabbit_print_impacts.py $RABBIT_OUTDIR/fitresults.hdf5 --impactType global -m Project ch0_masked --ungroup --sort
 
       - name: plot impacts
         run: >-
@@ -362,14 +362,14 @@ jobs:
         run: >-
           rabbit_plot_pulls_and_impacts.py $RABBIT_OUTDIR/fitresults.hdf5 -o $WEB_DIR/$PLOT_DIR
           -r $RABBIT_OUTDIR/fitresults_pseudodata.hdf5 --refResult original --refName "Pseudo data" -s absimpact
-          --otherExtensions pdf png -n 50  --config tests/style_config.py --oneSidedImpacts --globalImpacts --diffPullAsym --showNumbers
+          --otherExtensions pdf png -n 50  --config tests/style_config.py --oneSidedImpacts --impactType global --diffPullAsym --showNumbers
           --grouping max --subtitle Preliminary --postfix data_vs_pseudodata
 
       - name: plot global impacts on observable
         run: >-
           rabbit_plot_pulls_and_impacts.py $RABBIT_OUTDIR/fitresults.hdf5 -o $WEB_DIR/$PLOT_DIR
           -s absimpact -m Project ch0_masked
-          --otherExtensions pdf png -n 50  --config tests/style_config.py --oneSidedImpacts --globalImpacts --showNumbers
+          --otherExtensions pdf png -n 50  --config tests/style_config.py --oneSidedImpacts --impactType global --showNumbers
           --grouping max --subtitle Preliminary --postfix ch0_masked
 
       - name: plot prefit distributions

bin/rabbit_fit.py

@@ -125,6 +125,12 @@ def make_parser():
         action="store_true",
         help="propagate global impacts on histogram bins (inclusive in processes)",
     )
+    parser.add_argument(
+        "--computeHistGaussianImpacts",
+        default=False,
+        action="store_true",
+        help="propagate gaussian global impacts on histogram bins (inclusive in processes)",
+    )
     parser.add_argument(
         "--computeVariations",
         default=False,
@@ -165,7 +171,13 @@ def make_parser():
         "--globalImpacts",
         default=False,
         action="store_true",
-        help="compute impacts in terms of variations of global observables (as opposed to nuisance parameters directly)",
+        help="compute impacts in terms of variations from the likelihood of global observables (as opposed to nuisance parameters directly)",
+    )
+    parser.add_argument(
+        "--gaussianGlobalImpacts",
+        default=False,
+        action="store_true",
+        help="compute impacts in terms of variations of global observables in the fully gaussian approximation (as opposed to nuisance parameters directly)",
     )
     parser.add_argument(
         "--globalImpactsDisableJVP",
@@ -219,6 +231,7 @@ def save_hists(args, mappings, fitter, ws, prefit=True, profile=False):
                 compute_cov=args.computeHistCov,
                 compute_chi2=not args.noChi2 and mapping.has_data,
                 compute_global_impacts=args.computeHistImpacts,
+                compute_gaussian_global_impacts=args.computeHistGaussianImpacts,
                 profile=profile,
             )
 
@@ -229,11 +242,13 @@ def save_hists(args, mappings, fitter, ws, prefit=True, profile=False):
                 cov=aux[1],
                 impacts=aux[2],
                 impacts_grouped=aux[3],
+                gaussian_impacts=aux[4],
+                gaussian_impacts_grouped=aux[5],
                 prefit=prefit,
             )
 
-            if aux[4] is not None:
-                ws.add_chi2(aux[4], aux[5], prefit, mapping)
+            if aux[-2] is not None:
+                ws.add_chi2(aux[-2], aux[-1], prefit, mapping)
 
         if args.saveHistsPerProcess and not mapping.skip_per_process:
             logger.info(f"Save processes histogram for {mapping.key}")
@@ -299,13 +314,13 @@ def fit(args, fitter, ws, dofit=True):
 
     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)
+        _, grad, hess = fitter.loss_val_grad_hess()
+        edmval, cov = fitter.edmval_cov(grad, hess)
         logger.info(f"edmval: {edmval}")
 
-        fitter.cov.assign(cov)
+        ws.add_cov_hist(cov)
 
+        fitter.cov.assign(cov)
         del cov
 
         if fitter.binByBinStat and fitter.diagnostics:
@@ -314,7 +329,7 @@ def fit(args, fitter, ws, dofit=True):
             # 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)
+            edmvalbeta = edmval_cov(gradbeta, hessbeta)
             logger.info(f"edmvalbeta: {edmvalbeta}")
 
         if args.doImpacts:
@@ -323,7 +338,18 @@ def fit(args, fitter, ws, dofit=True):
         del hess
 
         if args.globalImpacts:
-            ws.add_global_impacts_hists(*fitter.global_impacts_parms())
+            ws.add_impacts_hists(
+                *fitter.global_impacts_parms(),
+                base_name="global_impacts",
+                global_impacts=True,
+            )
+
+        if args.gaussianGlobalImpacts:
+            ws.add_impacts_hists(
+                *fitter.gaussian_global_impacts_parms(),
+                base_name="gaussian_global_impacts",
+                global_impacts=True,
+            )
 
     nllvalreduced = fitter.reduced_nll().numpy()
 
@@ -359,12 +385,11 @@ def fit(args, fitter, ws, dofit=True):
         hist_name="parms",
     )
 
-    if not args.noHessian:
-        ws.add_cov_hist(fitter.cov)
-
     if args.nonProfiledImpacts:
         # TODO: based on covariance
-        ws.add_nonprofiled_impacts_hist(*fitter.nonprofiled_impacts_parms())
+        ws.add_impacts_asym_hist(
+            *fitter.nonprofiled_impacts_parms(), base_name="nonprofiled_impacts_asym"
+        )
 
     # Likelihood scans
     if args.scan is not None:

bin/rabbit_plot_hists_uncertainties.py

@@ -222,6 +222,12 @@ def parseArgs():
         action="store_true",
         help="Don't allow scientific notation for y axis",
     )
+    parser.add_argument(
+        "--impactType",
+        default="global",
+        choices=["gaussian_global", "global"],
+        help="Impact definition",
+    )
     args = parser.parse_args()
 
     return args
@@ -381,11 +387,13 @@ def make_plots(
     args=None,
     channel="",
     lumi=1,
+    impacts_name="global_impacts_grouped",
     *opts,
     **kwopts,
 ):
-    if "hist_postfit_inclusive_global_impacts_grouped" in result.keys():
-        hist_impacts = result["hist_postfit_inclusive_global_impacts_grouped"].get()
+    hist_impacts_name = f"hist_postfit_inclusive_{impacts_name}"
+    if hist_impacts_name in result.keys():
+        hist_impacts = result[hist_impacts_name].get()
     else:
         hist_impacts = None
 
@@ -559,12 +567,15 @@ def main():
                 ]:
                     suffix = suffix.replace(sign, rpl)
 
+                impacts_name = f"{args.impactType}_impacts_grouped"
+
                 make_plots(
                     outdir,
                     result,
                     config,
                     channel=suffix,
                     lumi=info.get("lumi", None),
+                    impacts_name=impacts_name,
                     **opts,
                 )
 

bin/rabbit_plot_pulls_and_impacts.py

@@ -455,7 +455,7 @@ def readFitInfoFromFile(
     fitresult,
     poi,
     group=False,
-    impact_type=False,
+    impact_type=None,
     grouping=None,
     asym=False,
     filters=[],
@@ -474,6 +474,7 @@ def readFitInfoFromFile(
             impact_type=impact_type,
             add_total=group and not impact_type == "nonprofiled",
         )
+
         if group:
             impacts, labels = out
             if normalize:
@@ -554,16 +555,19 @@ def readFitInfoFromFile(
         df["abspull"] = np.abs(df["pull"])
 
         if asym:
-            df["constraint_down"] = -constraints[..., 1]
-            df["constraint_up"] = constraints[..., 0]
+            if impact_type == "nonprofiled":
+                df["constraint_down"] = -constraints
+                df["constraint_up"] = constraints
+            else:
+                df["constraint_down"] = -constraints[..., 1]
+                df["constraint_up"] = constraints[..., 0]
         else:
             df["constraint"] = constraints
             valid = (1 - constraints**2) > 0
             df["newpull"] = 999.0
             df.loc[valid, "newpull"] = df.loc[valid]["pull"] / np.sqrt(
                 1 - df.loc[valid]["constraint"] ** 2
             )
-
     if poi:
         df = df.drop(df.loc[df["label"] == poi].index)
 
@@ -850,12 +854,19 @@ def parseArgs():
         type=str,
         help="Title for impacts",
     )
-    parser.add_argument("--noImpacts", action="store_true", help="Don't show impacts")
-    parser.add_argument(
-        "--globalImpacts", action="store_true", help="Print global impacts"
-    )
     parser.add_argument(
-        "--nonprofiledImpacts", action="store_true", help="Print non-profiled impacts"
+        "--impactType",
+        type=str,
+        default="traditional",
+        choices=[
+            None,
+            "none",
+            "traditional",
+            "global",
+            "gaussian_global",
+            "nonprofiled",
+        ],
+        help="Impact definition",
     )
     parser.add_argument(
         "--asymImpacts",
@@ -937,6 +948,7 @@ def make_plots(
     pullrange=None,
     meta=None,
     postfix=None,
+    impacts=False,
     impact_title=None,
 ):
 
@@ -980,16 +992,14 @@ def make_plots(
         pullrange=pullrange,
         title=args.title,
         subtitle=args.subtitle,
-        impacts=not args.noImpacts,
+        impacts=impacts,
         asym=asym,
         asym_pulls=args.diffPullAsym,
         include_ref=include_ref,
         ref_name=args.refName,
         name=args.name,
         show_numbers=args.showNumbers,
-        show_legend=(not group and not args.noImpacts)
-        or include_ref
-        or args.name is not None,
+        show_legend=(not group and impacts) or include_ref or args.name is not None,
         group=group,
         diff_pulls=not args.pullsNoDiff,
     )
@@ -1019,15 +1029,9 @@ def load_dataframe_parms(
     normalize=False,
     fitresult_ref=None,
     grouping=None,
+    impact_type=None,
     translate_label={},
 ):
-    if args.globalImpacts:
-        impact_type = "global"
-    elif args.nonprofiledImpacts:
-        impact_type = "nonprofiled"
-    else:
-        impact_type = "traditional"
-
     if not group:
         df = readFitInfoFromFile(
             fitresult,
@@ -1176,6 +1180,7 @@ def produce_plots_parms(
     pullrange=None,
     meta=None,
     postfix=None,
+    impact_type=None,
     impact_title=None,
     grouping=None,
     translate_label={},
@@ -1190,6 +1195,7 @@ def produce_plots_parms(
         normalize=normalize,
         fitresult_ref=fitresult_ref,
         grouping=grouping,
+        impact_type=impact_type,
         translate_label=translate_label,
     )
 
@@ -1206,6 +1212,7 @@ def produce_plots_parms(
         pullrange=pullrange,
         meta=meta,
         postfix=postfix,
+        impacts=impact_type not in [None, "none"],
         impact_title=impact_title,
     )
 
@@ -1260,6 +1267,7 @@ def produce_plots_hist(
         pullrange=pullrange,
         meta=meta,
         postfix=postfix,
+        impacts=hist_impacts is not None,
         impact_title=impact_title,
     )
 
@@ -1300,17 +1308,15 @@ def main():
         translate_label=translate_label,
     )
 
-    if args.noImpacts:
+    if args.impactType in [None, "none"]:
         # do one pulls plot, ungrouped
-        produce_plots_parms(args, fitresult, outdir, outfile="pulls.html", **kwargs)
+        produce_plots_parms(
+            args, fitresult, outdir, outfile="pulls.html", impact_type=None, **kwargs
+        )
     else:
         kwargs.update(dict(normalize=args.normalize, impact_title=args.impactTitle))
 
-        impacts_name = "impacts"
-        if args.globalImpacts:
-            impacts_name = f"global_{impacts_name}"
-        elif args.nonprofiledImpacts:
-            impacts_name = f"nonprofiled_{impacts_name}"
+        impacts_name = f"{args.impactType}_impacts"
 
         grouping = None
         if args.grouping is not None:
@@ -1321,9 +1327,9 @@ def main():
                 raise NotImplementedError(
                     "Asymetric impacts on observables is not yet implemented"
                 )
-            if not args.globalImpacts:
+            if args.impactType not in ["global", "gaussian_global"]:
                 raise NotImplementedError(
-                    "Only global impacts on observables is implemented (use --globalImpacts)"
+                    "Only global impacts on observables is implemented (use '--impactType' with 'global' or 'gaussian_global')"
                 )
 
             def get_mapping_key(result, key):
@@ -1369,7 +1375,7 @@ def get_mapping_key(result, key):
                 for mode in modes:
                     group = mode == "group"
 
-                    key = "hist_postfit_inclusive_global_impacts"
+                    key = f"hist_postfit_inclusive_{impacts_name}"
                     if group:
                         key += "_grouped"
 
@@ -1436,7 +1442,13 @@ def get_mapping_key(result, key):
                     if not args.noPulls:
                         name = f"pulls_and_{name}"
                     produce_plots_parms(
-                        args, fitresult, outdir, outfile=name, poi=poi, **kwargs
+                        args,
+                        fitresult,
+                        outdir,
+                        outfile=name,
+                        poi=poi,
+                        impact_type=args.impactType,
+                        **kwargs,
                     )
                 if args.mode in ["both", "group"]:
                     produce_plots_parms(
@@ -1447,6 +1459,7 @@ def get_mapping_key(result, key):
                         poi=poi,
                         group=True,
                         grouping=grouping,
+                        impact_type=args.impactType,
                         **kwargs,
                     )