feat: Improve handling of parameter bounds when computing impacts in ranking

src/cabinetry/fit/__init__.py

@@ -551,6 +551,9 @@ def ranking(
     maxiter: Optional[int] = None,
     tolerance: Optional[float] = None,
     custom_fit: bool = False,
+    minos: Optional[Union[List[str], Tuple[str, ...]]] = None,
+    minos_cl: Optional[float] = None,
+    vary_within_bounds: bool = False,
 ) -> RankingResults:
     """Calculates the impact of nuisance parameters on the parameter of interest (POI).
 
@@ -582,6 +585,17 @@ def ranking(
             None (use ``iminuit`` default of 0.1)
         custom_fit (bool, optional): whether to use the ``pyhf.infer`` API or
             ``iminuit``, defaults to False (using ``pyhf.infer``)
+        minos (Optional[Union[List[str], Tuple[str, ...]]], optional): runs the MINOS
+            algorithm for all parameters specified, defaults to None (does not run
+            MINOS). Used only if fit_results is not provided.
+        minos_cl (Optional[float]), optional): confidence level for the MINOS confidence
+            interval, defaults to None (use ``iminuit`` default of 68.27%). Used only
+            if fit_results is not provided.
+        vary_within_bounds (bool, optional): if parameter bounds are not specified,
+            this option specifies that the ``pyhf`` suggested parameter bounds should
+            be used to limit values of the parameters during ranking. Useful if one
+            of the ranking fits is failing. If some parameter bounds are specified,
+            the suggested bound will not be used. Default is False.
 
     Raises:
         ValueError: if no POI is found
@@ -600,6 +614,8 @@ def ranking(
             maxiter=maxiter,
             tolerance=tolerance,
             custom_fit=custom_fit,
+            minos=minos,
+            minos_cl=minos_cl,
         )
 
     labels = model.config.par_names
@@ -618,6 +634,27 @@ def ranking(
     init_pars = init_pars or model.config.suggested_init()
     fix_pars = fix_pars or model.config.suggested_fixed()
 
+    if par_bounds is None and vary_within_bounds:
+        log.info(
+            "No parameter bounds specified. Falling back to suggested bounds from pyhf."
+        )
+        par_bounds = model.config.suggested_bounds()
+    elif par_bounds is None and not vary_within_bounds:
+        log.warning(
+            "No parameter bounds specified and suggested bounds are disabled."
+            + " Ranking fits might be unstable."
+        )
+    else:
+        if (shape_par := np.asarray(par_bounds).shape) != (
+            shape_suggested := np.asarray(model.config.suggested_bounds()).shape
+        ):
+            raise ValueError(
+                f"Number of parameter bounds provided ({shape_par}) "
+                + f"does not match number of parameters ({shape_suggested})."
+            )
+        else:
+            log.info("Parameter bounds are specified and will be used for ranking.")
+
     all_impacts = []
     for i_par, label in enumerate(labels):
         if i_par == poi_index:
@@ -628,21 +665,42 @@ def ranking(
         fix_pars_ranking = fix_pars.copy()
         fix_pars_ranking[i_par] = True
 
+        # Get post-fit uncertainty from MINOS if available
+        postfit_unc_up = fit_results.uncertainty[i_par]
+        postfit_unc_do = -1 * fit_results.uncertainty[i_par]
+        # walrus assignment to avoid multiple lookups
+        if minos_unc := fit_results.minos_uncertainty.get(label):
+            postfit_unc_up = minos_unc[1]
+            postfit_unc_do = minos_unc[0]
+
         parameter_impacts = []
         # calculate impacts: pre-fit up, pre-fit down, post-fit up, post-fit down
-        for np_val in [
-            fit_results.bestfit[i_par] + prefit_unc[i_par],
-            fit_results.bestfit[i_par] - prefit_unc[i_par],
-            fit_results.bestfit[i_par] + fit_results.uncertainty[i_par],
-            fit_results.bestfit[i_par] - fit_results.uncertainty[i_par],
-        ]:
+        fit_labels = ["pre-fit up", "pre-fit down", "post-fit up", "post-fit down"]
+        for fit_idx, np_val in enumerate(
+            [
+                fit_results.bestfit[i_par] + prefit_unc[i_par],
+                fit_results.bestfit[i_par] - prefit_unc[i_par],
+                fit_results.bestfit[i_par] + postfit_unc_up,
+                fit_results.bestfit[i_par] + postfit_unc_do,
+            ]
+        ):
+            fit_label = fit_labels[fit_idx]
             # can skip pre-fit calculation for unconstrained parameters (their
             # pre-fit uncertainty is set to 0), and pre- and post-fit calculation
             # for fixed parameters (both uncertainties set to 0 as well)
             if np_val == fit_results.bestfit[i_par]:
                 log.debug(f"impact of {label} is zero, skipping fit")
                 parameter_impacts.append(0.0)
             else:
+                if par_bounds is not None:
+                    if not par_bounds[i_par][0] <= np_val <= par_bounds[i_par][1]:
+                        log.info(
+                            f"Parameter {label} value is out of bounds in the"
+                            + f" {fit_label} fit. Fixing it to its boundary."
+                        )
+                        np_val = min(
+                            max(np_val, par_bounds[i_par][0]), par_bounds[i_par][1]
+                        )
                 init_pars_ranking = init_pars.copy()
                 init_pars_ranking[i_par] = np_val  # value of current nuisance parameter
                 fit_results_ranking = _fit_model(

tests/fit/test_fit.py

@@ -622,9 +622,63 @@ def test_fit(mock_fit, mock_print, mock_gof):
         fit.FitResults(
             np.asarray([0.7, 0.9]), np.asarray([0.1, 0.1]), ["a", "b"], np.empty(0), 0.0
         ),
+        # for fourth ranking call without reference results with param bounds and
+        # vary_within_bounds used together
+        fit.FitResults(
+            np.asarray([1.3, 0.9]), np.asarray([0.1, 0.1]), ["a", "b"], np.empty(0), 0.0
+        ),
+        fit.FitResults(
+            np.asarray([0.7, 0.9]), np.asarray([0.1, 0.1]), ["a", "b"], np.empty(0), 0.0
+        ),
+        fit.FitResults(
+            np.asarray([1.2, 0.9]), np.asarray([0.1, 0.1]), ["a", "b"], np.empty(0), 0.0
+        ),
+        fit.FitResults(
+            np.asarray([0.8, 0.9]), np.asarray([0.1, 0.1]), ["a", "b"], np.empty(0), 0.0
+        ),
+        # for fifth ranking call with reference result including minos
+        fit.FitResults(
+            np.asarray([1.3, 0.9]), np.asarray([0.1, 0.1]), ["a", "b"], np.empty(0), 0.0
+        ),
+        fit.FitResults(
+            np.asarray([0.7, 0.9]), np.asarray([0.1, 0.1]), ["a", "b"], np.empty(0), 0.0
+        ),
+        fit.FitResults(
+            np.asarray([1.2, 0.9]), np.asarray([0.1, 0.1]), ["a", "b"], np.empty(0), 0.0
+        ),
+        fit.FitResults(
+            np.asarray([0.8, 0.9]), np.asarray([0.1, 0.1]), ["a", "b"], np.empty(0), 0.0
+        ),
+        # for sixth ranking call with reference result including minos
+        fit.FitResults(
+            np.asarray([1.3, 0.9]), np.asarray([0.1, 0.1]), ["a", "b"], np.empty(0), 0.0
+        ),
+        fit.FitResults(
+            np.asarray([0.7, 0.9]), np.asarray([0.1, 0.1]), ["a", "b"], np.empty(0), 0.0
+        ),
+        fit.FitResults(
+            np.asarray([1.2, 0.9]), np.asarray([0.1, 0.1]), ["a", "b"], np.empty(0), 0.0
+        ),
+        fit.FitResults(
+            np.asarray([0.8, 0.9]), np.asarray([0.1, 0.1]), ["a", "b"], np.empty(0), 0.0
+        ),
+        # for seventh ranking call with reference result including minos
+        fit.FitResults(
+            np.asarray([1.3, 0.9]), np.asarray([0.1, 0.1]), ["a", "b"], np.empty(0), 0.0
+        ),
+        fit.FitResults(
+            np.asarray([0.7, 0.9]), np.asarray([0.1, 0.1]), ["a", "b"], np.empty(0), 0.0
+        ),
+        fit.FitResults(
+            np.asarray([1.2, 0.9]), np.asarray([0.1, 0.1]), ["a", "b"], np.empty(0), 0.0
+        ),
+        fit.FitResults(
+            np.asarray([0.8, 0.9]), np.asarray([0.1, 0.1]), ["a", "b"], np.empty(0), 0.0
+        ),
     ],
 )
-def test_ranking(mock_fit, example_spec):
+def test_ranking(mock_fit, example_spec, caplog):
+    caplog.set_level(logging.DEBUG)
     example_spec["measurements"][0]["config"]["parameters"][0]["fixed"] = False
     bestfit = np.asarray([1.0, 0.9])
     uncertainty = np.asarray([0.1, 0.02])
@@ -649,6 +703,13 @@ def test_ranking(mock_fit, example_spec):
         assert mock_fit.call_args_list[i][1]["tolerance"] is None
         assert mock_fit.call_args_list[i][1]["custom_fit"] is False
 
+    assert (
+        "No parameter bounds specified and suggested bounds are disabled."
+        + " Ranking fits might be unstable."
+        in [rec.message for rec in caplog.records]
+    )
+    caplog.clear()
+
     # POI removed from fit results
     assert np.allclose(ranking_results.bestfit, [0.9])
     assert np.allclose(ranking_results.uncertainty, [0.02])
@@ -680,7 +741,7 @@ def test_ranking(mock_fit, example_spec):
     assert np.allclose(ranking_results.postfit_up, [0.2])
     assert np.allclose(ranking_results.postfit_down, [-0.2])
 
-    # no reference results, init/fixed pars, par bounds, strategy/maxiter/tolerance
+    # no fit results, init/fixed_pars, bounds, strategy/maxiter/tolerance,minos
     ranking_results = fit.ranking(
         model,
         data,
@@ -692,6 +753,8 @@ def test_ranking(mock_fit, example_spec):
         maxiter=100,
         tolerance=0.01,
         custom_fit=True,
+        minos=["Signal strength", "staterror_Signal-Region[0]"],
+        minos_cl=0.95,
     )
     assert mock_fit.call_count == 9
     # reference fit
@@ -705,6 +768,8 @@ def test_ranking(mock_fit, example_spec):
             "maxiter": 100,
             "tolerance": 0.01,
             "custom_fit": True,
+            "minos": ["Signal strength", "staterror_Signal-Region[0]"],
+            "minos_cl": 0.95,
         },
     )
     # fits for impact (comparing each option separately since init_pars needs allclose)
@@ -734,6 +799,103 @@ def test_ranking(mock_fit, example_spec):
     with pytest.raises(ValueError, match="no POI specified, cannot calculate ranking"):
         fit.ranking(model, data, fit_results=fit_results)
 
+    # test parameter bounding
+    # parameter bounds fully specified
+    caplog.clear()
+    example_spec["measurements"][0]["config"]["parameters"][0]["fixed"] = False
+    example_spec["measurements"][0]["config"]["poi"] = "Signal strength"
+    model, data = model_utils.model_and_data(example_spec)
+    ranking_results = fit.ranking(
+        model,
+        data,
+        fit_results=fit_results,
+        par_bounds=[(0, 5), (0.1, 10)],
+    )
+
+    assert "Parameter bounds are specified and will be used for ranking." in [
+        rec.message for rec in caplog.records
+    ]
+    caplog.clear()
+    assert mock_fit.call_count == 13
+
+    # parameter bounds partially specified and vary_within_bounds used
+    with pytest.raises(
+        ValueError,
+        match=r"Number of parameter bounds provided \(\(1, 2\)\) "
+        + r"does not match number of parameters \(\(2, 2\)\).",
+    ):
+        fit.ranking(
+            model,
+            data,
+            fit_results=fit_results,
+            par_bounds=[(0, 5)],
+            vary_within_bounds=True,
+        )
+
+    caplog.clear()
+    assert mock_fit.call_count == 13
+
+    # parameter bounds not specified and vary_within_bounds is used
+    ranking_results = fit.ranking(
+        model,
+        data,
+        fit_results=fit_results,
+        vary_within_bounds=True,
+    )
+    assert (
+        "No parameter bounds specified. Falling back to suggested bounds from pyhf."
+        in [rec.message for rec in caplog.records]
+    )
+    caplog.clear()
+    assert mock_fit.call_count == 17
+
+    # parameter bounds specified but np_value for a ranking fit is out of bound
+    fit_results = fit.FitResults(
+        np.asarray([1.0, 0.9]), np.asarray([0.1, 1.1]), ["a", "b"], np.empty(0), 0.0
+    )
+    ranking_results = fit.ranking(
+        model,
+        data,
+        fit_results=fit_results,
+        vary_within_bounds=True,
+    )
+    assert (
+        "Parameter staterror_Signal-Region[0] value is out of bounds in the"
+        + " post-fit down fit. Fixing it to its boundary."
+        in [rec.message for rec in caplog.records]
+    )
+    caplog.clear()
+    assert mock_fit.call_count == 21
+
+    fit_results = fit.FitResults(
+        bestfit,
+        uncertainty,
+        labels,
+        np.empty(0),
+        0.0,
+        0.0,
+        {"Signal strength": (-0.28, 0.32), "staterror_Signal-Region[0]": (-0.18, 0.22)},
+    )
+    ranking_results = fit.ranking(model, data, fit_results=fit_results)
+    # post-fit down, POI (stays same because we only vary NPs)
+    assert np.allclose(
+        mock_fit.call_args_list[-1][1]["init_pars"][0], (fit_results.bestfit[0])
+    )
+    # post-fit up, POI (stays same because we only vary NPs)
+    assert np.allclose(
+        mock_fit.call_args_list[-2][1]["init_pars"][0], (fit_results.bestfit[0])
+    )
+    # post-fit down, staterror
+    assert np.allclose(
+        mock_fit.call_args_list[-1][1]["init_pars"][1], (fit_results.bestfit[1] + -0.18)
+    )
+    # post-fit up, staterror
+    assert np.allclose(
+        mock_fit.call_args_list[-2][1]["init_pars"][1], (fit_results.bestfit[1] + 0.22)
+    )
+    caplog.clear()
+    assert mock_fit.call_count == 25
+
 
 @mock.patch(
     "cabinetry.fit._fit_model",