fixes for calling hesse and minos without existing function minimum (#…

…668)

doc/changelog.rst

@@ -7,6 +7,13 @@ Changelog
 
 2.7.0
 -----
+API change
+~~~~~~~~~~
+- If ``Minuit.hesse`` is called when ``Minuit.fmin`` is ``None``, an instance
+  ``Minuit.fmin`` is now created. If Hesse fails, the code does not raise an exception
+  anymore, since now the error state can be accessed as usual from the ``Minuit.fmin``
+  object. Users who relied on the exception should check ``Minuit.fmin`` from now on.
+
 New features
 ~~~~~~~~~~~~
 - ``Minuit.scipy`` can be used to minimise with SciPy algorithms. These may succeed
@@ -30,6 +37,10 @@ Fixes
   (patch submitted to ROOT)
 - ``Minuit.hesse`` no longer sets the status of the FunctionMinimum unconditionally
   to valid if it was invalid before
+- Repeated calls to ``Minuit.hesse`` no longer accumulate calls and eventually exhaust
+  the call limit, the call counter is now properly reset
+- Calling ``Minuit.minos`` repeatedly now does not recompute the Hessian and avoids
+  a bug that used to exhaust the call limit before in this case
 
 Documentation
 ~~~~~~~~~~~~~

src/hesse.cpp

@@ -10,7 +10,8 @@ using namespace ROOT::Minuit2;
 
 void update_fmin(MnHesse& self, const FCNBase& fcn, FunctionMinimum& min,
                  unsigned maxcalls, float maxedm) {
-  MnUserFcn mfcn(fcn, min.UserState().Trafo(), min.NFcn());
+  // We reset call counter here in contrast to MnHesse.cxx:83
+  MnUserFcn mfcn(fcn, min.UserState().Trafo(), 0);
 
   // Run Hesse
   MinimumState st = self(mfcn, min.State(), min.UserState().Trafo(), maxcalls);
@@ -32,15 +33,6 @@ void bind_hesse(py::module m) {
   py::class_<MnHesse>(m, "MnHesse")
 
       .def(py::init<const MnStrategy&>())
-      // pybind11 needs help to figure out the return value that's why we use lambdas
-      .def(
-          "__call__",
-          [](MnHesse& self, const FCNBase& fcn, const MnUserParameterState& state,
-             unsigned maxcalls) -> MnUserParameterState {
-            return self(fcn, state, maxcalls);
-          },
-          py::keep_alive<1, 2>())
-
       .def("__call__", update_fmin)
 
       ;

src/iminuit/minuit.py

@@ -1251,37 +1251,39 @@ def hesse(self, ncall: Optional[int] = None) -> "Minuit":
             )
             return self
 
-        hesse = MnHesse(self.strategy)
-
-        if self._fmin is None:
-            # _fmin does not exist or _last_state was modified,
-            # so we cannot just update last _fmin
-            self._last_state = hesse(self._fcn, self._last_state, ncall)
+        if self._fmin is None or self._fmin._src.state is not self._last_state:
+            # _fmin does not exist or last_state was modified, create a seed minimum
+            edm_goal = self._edm_goal(migrad_factor=True)
+            fm = FunctionMinimum(
+                self._fcn,
+                self._last_state,
+                self._strategy,
+                edm_goal,
+            )
+            self._fmin = mutil.FMin(
+                fm,
+                "External",
+                self.nfcn,
+                self.ngrad,
+                self.ndof,
+                edm_goal,
+            )
             self._merrors = mutil.MErrors()
-        else:
-            fm = self._fmin._src
-            if fm.state is self._last_state:
-                # fmin exists and _last_state not modified,
-                # can update _fmin which is more efficient
-                hesse(self._fcn, fm, ncall, self._fmin.edm_goal)
-                self._last_state = fm.state
-                self._fmin = mutil.FMin(
-                    fm,
-                    self._fmin.algorithm,
-                    self.nfcn,
-                    self.ngrad,
-                    self.ndof,
-                    self._fmin.edm_goal,
-                )
-            else:
-                # fmin exists but _last_state was modified
-                self._fmin = None
-                self._last_state = hesse(self._fcn, self._last_state, ncall)
-                self._merrors = mutil.MErrors()
 
-        if self._last_state.has_covariance is False:
-            if not self._fmin:
-                raise RuntimeError("Hesse Failed")
+        fm = self._fmin._src
+
+        # update _fmin with Hesse
+        hesse = MnHesse(self.strategy)
+        hesse(self._fcn, fm, ncall, self._fmin.edm_goal)
+        self._last_state = fm.state
+        self._fmin = mutil.FMin(
+            fm,
+            self._fmin.algorithm,
+            self.nfcn,
+            self.ngrad,
+            self.ndof,
+            self._fmin.edm_goal,
+        )
 
         self._make_covariance()
 
@@ -1335,28 +1337,18 @@ def minos(
         else:
             try:
                 from scipy.stats import chi2
-            except ImportError:  # pragma: no cover
-                raise ImportError("setting cl requires scipy")  # pragma: no cover
+            except ModuleNotFoundError as exc:  # pragma: no cover
+                exc.msg += "\nPlease install scipy to set the cl argument."
+                raise
             factor = chi2(1).ppf(cl)
 
-        if not self._fmin:
-            # create a seed minimum for MnMinos
-            fm = FunctionMinimum(
-                self._fcn,
-                self._last_state,
-                self._strategy,
-                self._edm_goal(migrad_factor=True),
-            )
-            # running MnHesse on seed is necessary for MnMinos to work
-            hesse = MnHesse(self.strategy)
-            hesse(self._fcn, fm, ncall, self._edm_goal(migrad_factor=True))
-            self._last_state = fm.state
-            self._make_covariance()
-        else:
-            fm = self._fmin._src
+        if not self._fmin or self._fmin._src.state is not self._last_state:
+            # _fmin does not exist or last_state was modified
+            self.hesse()  # also creates self._fmin
+        fm = self._fmin._src
 
         if not fm.is_valid:
-            raise RuntimeError("Function minimum is not valid.")
+            raise RuntimeError(f"Function minimum is not valid: {repr(self._fmin)}")
 
         if len(parameters) == 0:
             pars = [par for par in self.parameters if not self.fixed[par]]

tests/test_core.py

@@ -78,6 +78,13 @@ def test_MnUserParameterState():
     assert st[1].lower_limit == 1
     assert st[1].upper_limit == 4
 
+    st2 = MnUserParameterState(st)
+    assert st2 == st
+
+    st2.set_value(0, 1.1)
+
+    assert st2 != st
+
 
 def test_MnMigrad():
     fcn = FCN(fn, None, False, 1)

tests/test_minuit.py

@@ -1081,12 +1081,13 @@ def test_hesse_without_migrad():
     m.values["x"] = 1
     m.hesse()
     assert m.errors["x"] == approx((1.0 / 14.0) ** 0.5, abs=1e-4)
-    assert m.fmin is None
+    assert m.fmin
 
     m = Minuit(lambda x: 0, 0)
     m.errordef = 1
-    with pytest.raises(RuntimeError):
-        m.hesse()
+    m.hesse()
+    assert not m.accurate
+    assert m.fmin.hesse_failed
 
 
 def test_edm_goal():
@@ -1516,3 +1517,24 @@ def test_call_limit_reached_in_hesse():
     m.migrad(ncall=200)
     assert m.fmin.has_reached_call_limit
     assert m.fmin.nfcn < 205
+
+
+def test_issue_643():
+    def fcn(x, y, z):
+        return (x - 2) ** 2 + (y - 3) ** 2 + (z - 4) ** 2
+
+    fcn.errordef = Minuit.LEAST_SQUARES
+
+    m = Minuit(fcn, x=2, y=3, z=4)
+    m.migrad()
+
+    m2 = Minuit(fcn, x=m.values["x"], y=m.values["y"], z=m.values["z"])
+    # used to call MnHesse when it was not needed and quickly exhaust call limit
+    for i in range(10):
+        m2.minos()
+
+    m2.reset()
+    # used to exhaust call limit, because calls to MnHesse did not reset call count
+    for i in range(10):
+        m2.values = m.values
+        m2.minos()