Repository: scikit-hep/cabinetry

Issue:
Add Minuit instance to `FitResults`
Adding the `iminuit.Minuit` instance to the `FitResults` container would allow interacting with that directly to access more of the `iminuit` functionality. This could be useful if users want to run other steps after a fit and would work similarly to `return_result_obj=True` in `pyhf.infer.mle.fit`.

Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/scikit-hep/cabinetry/pull/504?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=scikit-hep) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Project coverage is 100.00%. Comparing base [(`ffa5ca2`)](https://app.codecov.io/gh/scikit-hep/cabinetry/commit/ffa5ca2a1626da0f56569c420977f8f97de368fe?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=scikit-hep) to head [(`afb84db`)](https://app.codecov.io/gh/scikit-hep/cabinetry/commit/afb84db0c430eca76e60fd0179db1ff2ba8eb801?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=scikit-hep).

<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff            @@
##            master      #504   +/-   ##
=========================================
  Coverage   100.00%   100.00%           
=========================================
  Files           22        22           
  Lines         2094      2096    +2     
  Branches       347       347           
=========================================
+ Hits          2094      2096    +2     
```

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/scikit-hep/cabinetry/pull/504?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=scikit-hep).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=scikit-hep).

MoAly98: @alexander-held any thoughts on whether or not additional tests are needed to validate the fit results container? I am not sure what kind of test would be useful to validate in the cases where we set the FitResults container manually. 
MoAly98: I implemented this as an optional argument, just so that it does not break usage of the `FitResults` container for users who directly access this API. 
MoAly98: @alexander-held I think these comments are now resolved and this is not a breaking change -- any other issues with it or we update and merge?

Patch:
diff --git a/src/cabinetry/fit/__init__.py b/src/cabinetry/fit/__init__.py
index caa91d2..d8fb79c 100644
--- a/src/cabinetry/fit/__init__.py
+++ b/src/cabinetry/fit/__init__.py
@@ -132,6 +132,7 @@ def _fit_model_pyhf(
         corr_mat,
         best_twice_nll,
         minos_uncertainty=minos_results,
+        minuit_obj=result_obj.minuit,
     )
     pyhf.set_backend(pyhf.tensorlib, initial_optimizer)  # restore optimizer settings
     return fit_results
@@ -251,6 +252,7 @@ def _fit_model_custom(
         corr_mat,
         best_twice_nll,
         minos_uncertainty=minos_results,
+        minuit_obj=m,
     )
     pyhf.set_backend(pyhf.tensorlib, initial_optimizer)  # restore optimizer settings
     return fit_results
diff --git a/src/cabinetry/fit/results_containers.py b/src/cabinetry/fit/results_containers.py
index f35f0a8..cafa405 100644
--- a/src/cabinetry/fit/results_containers.py
+++ b/src/cabinetry/fit/results_containers.py
@@ -1,7 +1,8 @@
 """Provides containers for inference results."""
 
-from typing import Dict, List, NamedTuple, Tuple
+from typing import Dict, List, NamedTuple, Optional, Tuple
 
+import iminuit
 import numpy as np
 
 
@@ -18,6 +19,7 @@ class FitResults(NamedTuple):
         goodess_of_fit (float, optional): goodness-of-fit p-value, defaults to -1
         minos_uncertainty (Dict[str, Tuple[float, float]]): uncertainties of best-fit
             parameter results indexed by parameter name, calculated with MINOS
+        minuit_obj (iminuit.Minuit, optional): underlying minimizer
     """
 
     bestfit: np.ndarray
@@ -27,6 +29,7 @@ class FitResults(NamedTuple):
     best_twice_nll: float
     goodness_of_fit: float = -1
     minos_uncertainty: Dict[str, Tuple[float, float]] = {}
+    minuit_obj: Optional[iminuit.Minuit] = None
 
 
 class RankingResults(NamedTuple):


Test patch:
diff --git a/tests/fit/test_fit.py b/tests/fit/test_fit.py
index edbdee3..8c23c0f 100644
--- a/tests/fit/test_fit.py
+++ b/tests/fit/test_fit.py
@@ -42,23 +42,52 @@ def test__fit_model_pyhf(mock_minos, example_spec, example_spec_multibin):
     pyhf.set_backend("numpy", "scipy")
     model, data = model_utils.model_and_data(example_spec)
     fit_results = fit._fit_model_pyhf(model, data)
+    assert fit_results.minuit_obj is not None
     assert np.allclose(fit_results.bestfit, [8.33624084, 1.1])
+    assert np.allclose(fit_results.minuit_obj.values, [8.33624084, 1.1])
     assert np.allclose(fit_results.uncertainty, [0.38182003, 0.0])
+    assert np.allclose(
+        np.where(fit_results.minuit_obj.fixed, 0.0, fit_results.minuit_obj.errors),
+        [0.38182003, 0.0],
+    )
     assert fit_results.labels == ["Signal strength", "staterror_Signal-Region[0]"]
+    assert fit_results.minuit_obj.parameters == (
+        "Signal strength",
+        "staterror_Signal-Region[0]",
+    )
     assert np.allclose(fit_results.best_twice_nll, 7.82495235)
+    assert np.allclose(fit_results.minuit_obj.fval, 7.82495235)
     assert np.allclose(fit_results.corr_mat, [[1.0, 0.0], [0.0, 0.0]])
+    assert np.allclose(
+        fit_results.minuit_obj.covariance.correlation(), [[1.0, 0.0], [0.0, 0.0]]
+    )
+
     assert pyhf.get_backend()[1].name == "scipy"  # optimizer was reset
 
     # Asimov fit, with fixed gamma (fixed not to Asimov MLE)
     model, data = model_utils.model_and_data(example_spec, asimov=True)
     fit_results = fit._fit_model_pyhf(model, data)
+    assert fit_results.minuit_obj is not None
     # the gamma factor is multiplicative and fixed to 1.1, so the
     # signal strength needs to be 1/1.1 to compensate
     assert np.allclose(fit_results.bestfit, [0.90917877, 1.1])
+    assert np.allclose(fit_results.minuit_obj.values, [0.90917877, 1.1])
     assert np.allclose(fit_results.uncertainty, [0.12628017, 0.0])
+    assert np.allclose(
+        np.where(fit_results.minuit_obj.fixed, 0.0, fit_results.minuit_obj.errors),
+        [0.12628017, 0.0],
+    )
     assert fit_results.labels == ["Signal strength", "staterror_Signal-Region[0]"]
+    assert fit_results.minuit_obj.parameters == (
+        "Signal strength",
+        "staterror_Signal-Region[0]",
+    )
     assert np.allclose(fit_results.best_twice_nll, 5.61189476)
+    assert np.allclose(fit_results.minuit_obj.fval, 5.61189476)
     assert np.allclose(fit_results.corr_mat, [[1.0, 0.0], [0.0, 0.0]])
+    assert np.allclose(
+        fit_results.minuit_obj.covariance.correlation(), [[1.0, 0.0], [0.0, 0.0]]
+    )
 
     # parameters held constant via keyword argument
     model, data = model_utils.model_and_data(example_spec_multibin)
@@ -71,15 +100,26 @@ def test__fit_model_pyhf(mock_minos, example_spec, example_spec_multibin):
     fit_results = fit._fit_model_pyhf(
         model, data, init_pars=init_pars, fix_pars=fix_pars
     )
+    assert fit_results.minuit_obj is not None
     assert np.allclose(fit_results.bestfit, [1.48041923, 0.9, 1.1, 0.97511112])
+    assert np.allclose(
+        fit_results.minuit_obj.values, [1.48041923, 0.9, 1.1, 0.97511112]
+    )
     assert np.allclose(fit_results.uncertainty, [0.20694409, 0.0, 0.0, 0.11792805])
+    assert np.allclose(
+        np.where(fit_results.minuit_obj.fixed, 0.0, fit_results.minuit_obj.errors),
+        [0.20694409, 0.0, 0.0, 0.11792805],
+    )
     assert np.allclose(fit_results.best_twice_nll, 10.4531891)
+    assert np.allclose(fit_results.minuit_obj.fval, 10.4531891)
 
     # custom parameter bounds
     model, data = model_utils.model_and_data(example_spec)
     par_bounds = [(0, 5), (0.1, 2)]
     fit_results = fit._fit_model_pyhf(model, data, par_bounds=par_bounds)
+    assert fit_results.minuit_obj is not None
     assert np.allclose(fit_results.bestfit, [5.0, 1.1])
+    assert np.allclose(fit_results.minuit_obj.values, [5.0, 1.1])
 
     # propagation of strategy, max iterations, tolerance
     model, data = model_utils.model_and_data(example_spec)
@@ -122,23 +162,51 @@ def test__fit_model_custom(mock_minos, example_spec, example_spec_multibin):
     pyhf.set_backend("numpy", "scipy")
     model, data = model_utils.model_and_data(example_spec)
     fit_results = fit._fit_model_custom(model, data)
+    assert fit_results.minuit_obj is not None
     assert np.allclose(fit_results.bestfit, [8.33625071, 1.1])
+    assert np.allclose(fit_results.minuit_obj.values, [8.33625071, 1.1])
     assert np.allclose(fit_results.uncertainty, [0.38182151, 0.0])
+    assert np.allclose(
+        np.where(fit_results.minuit_obj.fixed, 0.0, fit_results.minuit_obj.errors),
+        [0.38182151, 0.0],
+    )
     assert fit_results.labels == ["Signal strength", "staterror_Signal-Region[0]"]
+    assert fit_results.minuit_obj.parameters == (
+        "Signal strength",
+        "staterror_Signal-Region[0]",
+    )
     assert np.allclose(fit_results.best_twice_nll, 7.82495235)
+    assert np.allclose(fit_results.minuit_obj.fval, 7.82495235)
     assert np.allclose(fit_results.corr_mat, [[1.0, 0.0], [0.0, 0.0]])
+    assert np.allclose(
+        fit_results.minuit_obj.covariance.correlation(), [[1.0, 0.0], [0.0, 0.0]]
+    )
     assert pyhf.get_backend()[1].name == "scipy"  # optimizer was reset
 
     # Asimov fit, with fixed gamma (fixed not to Asimov MLE)
     model, data = model_utils.model_and_data(example_spec, asimov=True)
     fit_results = fit._fit_model_custom(model, data)
+    assert fit_results.minuit_obj is not None
     # the gamma factor is multiplicative and fixed to 1.1, so the
     # signal strength needs to be 1/1.1 to compensate
     assert np.allclose(fit_results.bestfit, [0.90917877, 1.1])
+    assert np.allclose(fit_results.minuit_obj.values, [0.90917877, 1.1])
     assert np.allclose(fit_results.uncertainty, [0.12628023, 0.0])
+    assert np.allclose(
+        np.where(fit_results.minuit_obj.fixed, 0.0, fit_results.minuit_obj.errors),
+        [0.12628023, 0.0],
+    )
     assert fit_results.labels == ["Signal strength", "staterror_Signal-Region[0]"]
+    assert fit_results.minuit_obj.parameters == (
+        "Signal strength",
+        "staterror_Signal-Region[0]",
+    )
     assert np.allclose(fit_results.best_twice_nll, 5.61189476)
+    assert np.allclose(fit_results.minuit_obj.fval, 5.61189476)
     assert np.allclose(fit_results.corr_mat, [[1.0, 0.0], [0.0, 0.0]])
+    assert np.allclose(
+        fit_results.minuit_obj.covariance.correlation(), [[1.0, 0.0], [0.0, 0.0]]
+    )
 
     # parameters held constant via keyword argument
     model, data = model_utils.model_and_data(example_spec_multibin)
@@ -151,15 +219,26 @@ def test__fit_model_custom(mock_minos, example_spec, example_spec_multibin):
     fit_results = fit._fit_model_custom(
         model, data, init_pars=init_pars, fix_pars=fix_pars
     )
+    assert fit_results.minuit_obj is not None
     assert np.allclose(fit_results.bestfit, [1.48041923, 0.9, 1.1, 0.97511112])
+    assert np.allclose(
+        fit_results.minuit_obj.values, [1.48041923, 0.9, 1.1, 0.97511112]
+    )
     assert np.allclose(fit_results.uncertainty, [0.20694409, 0.0, 0.0, 0.11792805])
+    assert np.allclose(
+        np.where(fit_results.minuit_obj.fixed, 0.0, fit_results.minuit_obj.errors),
+        [0.20694409, 0.0, 0.0, 0.11792805],
+    )
     assert np.allclose(fit_results.best_twice_nll, 10.45318909)
+    assert np.allclose(fit_results.minuit_obj.fval, 10.45318909)
 
     # custom parameter bounds
     model, data = model_utils.model_and_data(example_spec)
     par_bounds = [(0, 5), (0.1, 2)]
     fit_results = fit._fit_model_custom(model, data, par_bounds=par_bounds)
+    assert fit_results.minuit_obj is not None
     assert np.allclose(fit_results.bestfit, [5.0, 1.1])
+    assert np.allclose(fit_results.minuit_obj.values, [5.0, 1.1])
 
     # propagation of strategy, max iterations, tolerance
     model, data = model_utils.model_and_data(example_spec)
@@ -436,7 +515,7 @@ def test_fit(mock_fit, mock_print, mock_gof):
     assert mock_print.call_args[0][0].uncertainty == [0.1]
     assert mock_print.call_args[0][0].labels == ["par"]
     assert mock_gof.call_count == 0
-    assert fit_results.bestfit == [1.0]
+    assert np.allclose(fit_results.bestfit, [1.0])
 
     # custom fit, init/fix pars, par bounds, strategy/maxiter/tolerance
     init_pars = [2.0]
@@ -470,7 +549,7 @@ def test_fit(mock_fit, mock_print, mock_gof):
     )
     assert mock_print.call_args[0][0].bestfit == [1.0]
     assert mock_print.call_args[0][0].uncertainty == [0.1]
-    assert fit_results.bestfit == [1.0]
+    assert np.allclose(fit_results.bestfit, [1.0])
 
     # parameters for MINOS
     fit_results = fit.fit(model, data, minos=["abc"], minos_cl=0.95)
@@ -486,7 +565,9 @@ def test_fit(mock_fit, mock_print, mock_gof):
         "tolerance": None,
         "custom_fit": False,
     }
-    assert fit_results.bestfit == [1.0]
+    assert np.allclose(fit_results.bestfit, [1.0])
+
+    # custom fit
     fit_results = fit.fit(model, data, minos="abc", minos_cl=0.95, custom_fit=True)
     assert mock_fit.call_count == 4
     assert mock_fit.call_args[1] == {
@@ -500,7 +581,7 @@ def test_fit(mock_fit, mock_print, mock_gof):
         "tolerance": None,
         "custom_fit": True,
     }
-    assert fit_results.bestfit == [1.0]
+    assert np.allclose(fit_results.bestfit, [1.0])
 
     # goodness-of-fit test
     fit_results_gof = fit.fit(model, data, goodness_of_fit=True, fix_pars=fix_pars)
diff --git a/tests/fit/test_fit_results_containers.py b/tests/fit/test_fit_results_containers.py
index 0739663..e630262 100644
--- a/tests/fit/test_fit_results_containers.py
+++ b/tests/fit/test_fit_results_containers.py
@@ -1,5 +1,6 @@
 import logging
 
+import iminuit
 import numpy as np
 
 from cabinetry import fit
@@ -11,7 +12,10 @@ def test_FitResults():
     labels = ["par_a"]
     corr_mat = np.asarray([[1.0]])
     best_twice_nll = 2.0
-    fit_results = fit.FitResults(bestfit, uncertainty, labels, corr_mat, best_twice_nll)
+    minuit_obj = iminuit.Minuit(lambda x: x, x=1)
+    fit_results = fit.FitResults(
+        bestfit, uncertainty, labels, corr_mat, best_twice_nll, minuit_obj=minuit_obj
+    )
     assert np.allclose(fit_results.bestfit, bestfit)
     assert np.allclose(fit_results.uncertainty, uncertainty)
     assert fit_results.labels == labels
@@ -19,6 +23,7 @@ def test_FitResults():
     assert fit_results.best_twice_nll == best_twice_nll
     assert fit_results.goodness_of_fit == -1
     assert fit_results.minos_uncertainty == {}
+    assert fit_results.minuit_obj == minuit_obj
 
 
 def test_RankingResults():