fix: handle negative model predictions in visualizations (#394)

* catch negative total model prediction in visualize.plot_model.data_mc
* handle negative nominal template prediction in visualize.plot_model.templates

---------

Co-authored-by: Roman Mueller <[email protected]>
Co-authored-by: Alexander Held <[email protected]>

Author: 3 people

src/cabinetry/visualize/plot_model.py

@@ -52,6 +52,9 @@ def data_mc(
             saving it, defaults to False (enable when producing many figures to avoid
             memory issues, prevents rendering in notebooks)
 
+    Raises:
+        ValueError: when total model yield is negative in any bin
+
     Returns:
         matplotlib.figure.Figure: the data/MC figure
     """
@@ -149,6 +152,11 @@ def data_mc(
         )
     nonzero_model_yield = total_yield != 0.0
 
+    if np.any(total_yield < 0.0):
+        raise ValueError(
+            f"{label} total model yield has negative bin(s): {total_yield.tolist()}"
+        )
+
     # add uncertainty band around y=1
     rel_mc_unc = total_model_unc / total_yield
     # do not show band in bins where total model yield is 0
@@ -312,6 +320,15 @@ def templates(
     # x positions for lines drawn showing the template distributions
     line_x = [y for y in bin_edges for _ in range(2)][1:-1]
 
+    neg_nom_bin = False  # negative bin(s) present in nominal histogram
+    if np.any(nominal_histo["yields"] < 0.0):
+        neg_nom_bin = True
+        log.warning(
+            f"{label} nominal histogram yield has negative bin(s): "
+            f"{nominal_histo['yields'].tolist()}, taking absolute value for "
+            "ratio plot uncertainty"
+        )
+
     # draw templates
     for template, color, linestyle, template_label in zip(
         all_templates, colors, linestyles, template_labels
@@ -354,7 +371,7 @@ def templates(
             ax2.errorbar(
                 bin_centers,
                 template_ratio_plot,
-                yerr=template["stdev"] / nominal_histo["yields"],
+                yerr=template["stdev"] / np.abs(nominal_histo["yields"]),
                 fmt="none",
                 color=color,
             )
@@ -395,7 +412,7 @@ def templates(
     ax2.set_xlim([bin_edges[0], bin_edges[-1]])
     ax2.set_ylim([0.5, 1.5])
     ax2.set_xlabel(variable)
-    ax2.set_ylabel("variation / nominal")
+    ax2.set_ylabel(f"variation / {'nominal' if not neg_nom_bin else 'abs(nominal)'}")
     ax2.set_yticks([0.5, 0.75, 1.0, 1.25, 1.5])
     ax2.set_yticklabels([0.5, 0.75, 1.0, 1.25, ""])
     ax2.tick_params(axis="both", which="major", pad=8)

tests/visualize/test_visualize_plot_model.py

@@ -125,10 +125,20 @@ def test_data_mc(tmp_path, caplog):
     # expect three RuntimeWarnings from numpy due to division by zero
     assert sum("divide by zero" in str(m.message) for m in warn_record) == 3
 
+    # negative bin yield
+    histo_dict_list[0]["yields"] = [-50, -50]
+    with pytest.raises(
+        ValueError,
+        match=r"abc total model yield has negative bin\(s\): \[-50, -45\]",
+    ):
+        plot_model.data_mc(
+            histo_dict_list, total_model_unc_log, bin_edges_log, label="abc"
+        )
+
     plt.close("all")
 
 
-def test_templates(tmp_path):
+def test_templates(tmp_path, caplog):
     fname = tmp_path / "fig.png"
     nominal_histo = {
         "yields": np.asarray([1.0, 1.2]),
@@ -172,9 +182,12 @@ def test_templates(tmp_path):
     assert (
         compare_images("tests/visualize/reference/templates.png", str(fname), 0) is None
     )
+    caplog.clear()
 
     # do not save figure, but close it
     # only single variation specified
+    # negative bin present in nominal histogram
+    nominal_histo["yields"][0] = -1
     with mock.patch("cabinetry.visualize.utils._save_and_close") as mock_close_safe:
         fig = plot_model.templates(
             nominal_histo,
@@ -188,6 +201,12 @@ def test_templates(tmp_path):
             close_figure=True,
         )
         assert mock_close_safe.call_args_list == [((fig, None, True), {})]
+    assert (
+        "region: Signal region\nsample: Signal\nsystematic: Modeling nominal histogram "
+        "yield has negative bin(s): [-1.0, 1.2], taking absolute value for ratio plot "
+        "uncertainty"
+    ) in [rec.message for rec in caplog.records]
+    caplog.clear()
 
     plt.close("all")