feat!: Implement covariance-based impact calculations

src/cabinetry/cli/__init__.py

@@ -154,8 +154,17 @@ def fit(
     default="figures",
     help='folder to save figures to (default: "figures")',
 )
+@click.option(
+    "--impacts_method",
+    default="covariance",
+    help="The method to be used for computing impacts",
+)
 def ranking(
-    ws_spec: io.TextIOWrapper, asimov: bool, max_pars: int, figfolder: str
+    ws_spec: io.TextIOWrapper,
+    asimov: bool,
+    max_pars: int,
+    figfolder: str,
+    impacts_method: str,
 ) -> None:
     """Ranks nuisance parameters and visualizes the result.
 
@@ -165,9 +174,13 @@ def ranking(
     ws = json.load(ws_spec)
     model, data = cabinetry_model_utils.model_and_data(ws, asimov=asimov)
     fit_results = cabinetry_fit.fit(model, data)
-    ranking_results = cabinetry_fit.ranking(model, data, fit_results=fit_results)
+    ranking_results = cabinetry_fit.ranking(
+        model, data, fit_results=fit_results, impacts_method=impacts_method
+    )
     cabinetry_visualize.ranking(
-        ranking_results, figure_folder=figfolder, max_pars=max_pars
+        ranking_results,
+        figure_folder=figfolder,
+        max_pars=max_pars,
     )
 
 

src/cabinetry/fit/results_containers.py

@@ -43,6 +43,7 @@ class RankingResults(NamedTuple):
         prefit_down (np.ndarray): pre-fit impact in "down" direction
         postfit_up (np.ndarray): post-fit impact in "up" direction
         postfit_down (np.ndarray): post-fit impact in "down" direction
+        impacts_method (str): the method used to compute parameter impacts
     """
 
     bestfit: np.ndarray
@@ -52,6 +53,7 @@ class RankingResults(NamedTuple):
     prefit_down: np.ndarray
     postfit_up: np.ndarray
     postfit_down: np.ndarray
+    impacts_method: str
 
 
 class ScanResults(NamedTuple):

src/cabinetry/visualize/__init__.py

@@ -566,7 +566,11 @@ def ranking(
         matplotlib.figure.Figure: the ranking figure
     """
     # path is None if figure should not be saved
-    figure_path = pathlib.Path(figure_folder) / "ranking.pdf" if save_figure else None
+    figure_path = (
+        pathlib.Path(figure_folder) / f"ranking_{ranking_results.impacts_method}.pdf"
+        if save_figure
+        else None
+    )
 
     # sort parameters by decreasing maximum post-fit impact
     max_postfit_impact = np.maximum(
@@ -603,6 +607,7 @@ def ranking(
         postfit_down,
         figure_path=figure_path,
         close_figure=close_figure,
+        impacts_method=ranking_results.impacts_method,
     )
     return fig
 

src/cabinetry/visualize/plot_result.py

@@ -128,6 +128,7 @@ def ranking(
     *,
     figure_path: Optional[pathlib.Path] = None,
     close_figure: bool = False,
+    impacts_method: str = "covariance",
 ) -> mpl.figure.Figure:
     """Draws a ranking plot.
 
@@ -146,6 +147,8 @@ def ranking(
         close_figure (bool, optional): whether to close each figure immediately after
             saving it, defaults to False (enable when producing many figures to avoid
             memory issues, prevents rendering in notebooks)
+        impacts_method (str, optional): the method used to compute parameter impacts.
+            Options are np_shift, covariance, auxdata_shift
 
     Returns:
         matplotlib.figure.Figure: the ranking figure
@@ -163,6 +166,22 @@ def ranking(
     else:
         layout = None  # pragma: no cover  # layout set after figure creation instead
 
+    if impacts_method not in ["np_shift", "covariance", "auxdata_shift"]:
+        raise ValueError(
+            f"The impacts method {impacts_method} provided is not supported."
+            + " Valid options are (np_shift, covariance, auxdata_shift)"
+        )
+    if impacts_method == "auxdata_shift":
+        raise NotImplementedError(
+            "Plotting impacts computed by shifting auxiliary data is not supported yet."
+        )
+
+    impacts_color_map = {
+        "np_shift": ["C0", "C5"],
+        "covariance": ["#2CA02C", "#98DF8A"],
+        "auxdata_shift": ["#9467BD", "#C5B0D5"],
+    }
+
     mpl.style.use(MPL_STYLE)
     fig, ax_pulls = plt.subplots(
         figsize=(8, 2.5 + num_pars * 0.45), dpi=100, layout=layout
@@ -197,18 +216,33 @@ def ranking(
 
     y_pos = np.arange(num_pars)[::-1]
 
-    # pre-fit up
-    pre_up = ax_impact.barh(
-        y_pos, impact_prefit_up, fill=False, linewidth=1, edgecolor="C0"
-    )
-    # pre-fit down
-    pre_down = ax_impact.barh(
-        y_pos, impact_prefit_down, fill=False, linewidth=1, edgecolor="C5"
-    )
+    pre_up, pre_down = None, None
+    if impacts_method == "np_shift":
+        # pre-fit up
+        pre_up = ax_impact.barh(
+            y_pos,
+            impact_prefit_up,
+            fill=False,
+            linewidth=1,
+            edgecolor=impacts_color_map[impacts_method][0],
+        )
+        # pre-fit down
+        pre_down = ax_impact.barh(
+            y_pos,
+            impact_prefit_down,
+            fill=False,
+            linewidth=1,
+            edgecolor=impacts_color_map[impacts_method][1],
+        )
+
     # post-fit up
-    post_up = ax_impact.barh(y_pos, impact_postfit_up, color="C0")
+    post_up = ax_impact.barh(
+        y_pos, impact_postfit_up, color=impacts_color_map[impacts_method][0]
+    )
     # post-fit down
-    post_down = ax_impact.barh(y_pos, impact_postfit_down, color="C5")
+    post_down = ax_impact.barh(
+        y_pos, impact_postfit_down, color=impacts_color_map[impacts_method][1]
+    )
     # nuisance parameter pulls
     pulls = ax_pulls.errorbar(bestfit, y_pos, xerr=uncertainty, fmt="o", color="k")
 
@@ -244,20 +278,44 @@ def ranking(
     ax_pulls.tick_params(direction="in", which="both")
     ax_impact.tick_params(direction="in", which="both")
 
-    fig.legend(
-        (pre_up, pre_down, post_up, post_down, pulls),
-        (
-            r"pre-fit impact: $\theta = \hat{\theta} + \Delta \theta$",
-            r"pre-fit impact: $\theta = \hat{\theta} - \Delta \theta$",
-            r"post-fit impact: $\theta = \hat{\theta} + \Delta \hat{\theta}$",
-            r"post-fit impact: $\theta = \hat{\theta} - \Delta \hat{\theta}$",
-            "pulls",
-        ),
-        frameon=False,
-        loc="upper left",
-        ncol=3,
-        fontsize="large",
+    leg_handlers = (
+        (pre_up, pre_down, post_up, post_down, pulls)
+        if impacts_method == "np_shift"
+        else (post_up, post_down, pulls)
     )
+    leg_settings = {
+        "frameon": False,
+        "loc": "upper left",
+        "ncol": 3,
+        "fontsize": "large",
+    }
+
+    if impacts_method == "np_shift":
+        fig.legend(
+            leg_handlers,
+            (
+                r"pre-fit impact: $\theta = \hat{\theta} + \Delta \theta$",
+                r"pre-fit impact: $\theta = \hat{\theta} - \Delta \theta$",
+                r"post-fit impact: $\theta = \hat{\theta} + \Delta \hat{\theta}$",
+                r"post-fit impact: $\theta = \hat{\theta} - \Delta \hat{\theta}$",
+                "pulls",
+            ),
+            **leg_settings,
+        )
+    elif impacts_method == "covariance":
+        fig.legend(
+            leg_handlers,
+            (
+                r"post-fit impact: $\sigma_{\text{POI}}\,$"
+                + r"$\rho(\text{POI},\theta)\,\sigma_{\theta}$",
+                r"post-fit impact: $-\sigma_{\text{POI}}\,$"
+                + r"$\rho(\text{POI},\theta)\,\sigma_{\theta}$",
+                "pulls",
+            ),
+            **leg_settings,
+        )
+    else:
+        pass  # to be implemented
 
     if not MPL_GEQ_36:
         fig.tight_layout(rect=[0, 0, 1.0, 1 - leg_space])  # pragma: no cover

tests/cli/test_cli.py

@@ -203,6 +203,7 @@ def test_fit(mock_utils, mock_fit, mock_pulls, mock_corrmat, tmp_path):
         np.asarray([[0.8]]),
         np.asarray([[1.1]]),
         np.asarray([[0.9]]),
+        impacts_method="np_shift",
     ),
     autospec=True,
 )
@@ -240,7 +241,10 @@ def test_ranking(mock_utils, mock_fit, mock_rank, mock_vis, tmp_path):
     assert mock_utils.call_args_list == [((workspace,), {"asimov": False})]
     assert mock_fit.call_args_list == [(("model", "data"), {})]
     assert mock_rank.call_args_list == [
-        (("model", "data"), {"fit_results": fit_results})
+        (
+            ("model", "data"),
+            {"fit_results": fit_results, "impacts_method": "covariance"},
+        )
     ]
     assert mock_vis.call_count == 1
     assert np.allclose(mock_vis.call_args[0][0].prefit_up, [[1.2]])
@@ -252,12 +256,24 @@ def test_ranking(mock_utils, mock_fit, mock_rank, mock_vis, tmp_path):
     # Asimov, maximum amount of parameters, custom folder
     result = runner.invoke(
         cli.ranking,
-        ["--asimov", "--max_pars", 3, "--figfolder", "folder", workspace_path],
+        [
+            "--impacts_method",
+            "np_shift",
+            "--asimov",
+            "--max_pars",
+            3,
+            "--figfolder",
+            "folder",
+            workspace_path,
+        ],
     )
     assert result.exit_code == 0
     assert mock_utils.call_args == ((workspace,), {"asimov": True})
     assert mock_fit.call_args == (("model", "data"), {})
-    assert mock_rank.call_args == (("model", "data"), {"fit_results": fit_results})
+    assert mock_rank.call_args == (
+        ("model", "data"),
+        {"fit_results": fit_results, "impacts_method": "np_shift"},
+    )
     assert mock_vis.call_args[1] == {"figure_folder": "folder", "max_pars": 3}