Allow running Minos without Migrad first (#590)

Author: HDembinski

doc/changelog.rst

@@ -12,6 +12,8 @@ New features
 ~~~~~~~~~~~~
 - Minuit object is now pickle-able and copy-able
 - More efficient internal conversion between Python objects and ``std::vector<double>``
+- ``Minuit.minos`` can now be called without calling ``Minuit.migrad`` first, which allows
+  one to use an external minimiser to find a minimum and then compute Minos errors for it
 
 Bug-fixes
 ~~~~~~~~~

src/functionminimum.cpp

@@ -18,7 +18,7 @@ using namespace ROOT::Minuit2;
 MinimumSeed make_seed(const FCN& fcn, const MnUserFcn& mfcn,
                       const MnUserParameterState& st, const MnStrategy& str);
 
-FunctionMinimum init(const FCN& fcn, const MnUserParameterState& st, double fval,
+FunctionMinimum init(const FCN& fcn, const MnUserParameterState& st,
                      const MnStrategy& str, double edm_goal);
 
 py::tuple fmin_getstate(const FunctionMinimum&);

src/functionminimum_extra.cpp

@@ -26,7 +26,7 @@ MinimumSeed make_seed(const FCN& fcn, const MnUserFcn& mfcn,
   return gen(mfcn, gc, st, str);
 }
 
-FunctionMinimum init(const FCN& fcn, const MnUserParameterState& st, double fval,
+FunctionMinimum init(const FCN& fcn, const MnUserParameterState& st,
                      const MnStrategy& str, double edm_goal) {
   MnUserFcn mfcn(fcn, st.Trafo());
   MinimumSeed seed = make_seed(fcn, mfcn, st, str);
@@ -39,7 +39,7 @@ FunctionMinimum init(const FCN& fcn, const MnUserParameterState& st, double fval
     err(i) = std::sqrt(2. * mfcn.Up() * seed.Error().InvHessian()(i, i));
   }
 
-  MinimumParameters minp(val, err, fval);
+  MinimumParameters minp(val, err, seed.Fval());
   std::vector<MinimumState> minstv(1, MinimumState(minp, seed.Edm(), fcn.nfcn_));
   if (minstv.back().Edm() < edm_goal) return FunctionMinimum(seed, minstv, fcn.Up());
   return FunctionMinimum(seed, minstv, fcn.Up(), FunctionMinimum::MnAboveMaxEdm());

src/iminuit/_minuit.py

@@ -4,6 +4,7 @@
     FCN,
     MnContours,
     MnHesse,
+    MnMachinePrecision,
     MnMigrad,
     MnMinos,
     MnPrint,
@@ -550,13 +551,7 @@ def migrad(self, ncall=None, iterate=5):
 
         self._last_state = fm.state
 
-        # EDM goal
-        # - taken from the source code, see VariableMeticBuilder::Minimum and
-        #   ModularFunctionMinimizer::Minimize
-        # - goal is used to detect convergence but violations by 10x are also accepted;
-        #   see VariableMetricBuilder.cxx:425
-        edm_goal = 2e-3 * max(self._tolerance * fm.errordef, migrad.precision.eps2)
-
+        edm_goal = self._migrad_edm_goal()
         self._fmin = mutil.FMin(fm, self.nfcn, self.ngrad, edm_goal)
         self._make_covariance()
 
@@ -605,6 +600,8 @@ def simplex(self, ncall=None):
 
         edm_goal = max(self._tolerance * fm.errordef, simplex.precision.eps2)
         self._fmin = mutil.FMin(fm, self.nfcn, self.ngrad, edm_goal)
+        self._covariance = None
+        self._merrors = mutil.MErrors()
 
         return self  # return self for method chaining and to autodisplay current state
 
@@ -694,11 +691,11 @@ def run(ipar):
         run(0)
 
         edm_goal = self._tolerance * self._fcn._errordef
-        fm = FunctionMinimum(
-            self._fcn, self._last_state, x[self.npar], self.strategy, edm_goal
-        )
+        fm = FunctionMinimum(self._fcn, self._last_state, self.strategy, edm_goal)
         self._last_state = fm.state
         self._fmin = mutil.FMin(fm, self.nfcn, self.ngrad, edm_goal)
+        self._covariance = None
+        self._merrors = mutil.MErrors()
 
         return self  # return self for method chaining and to autodisplay current state
 
@@ -806,13 +803,16 @@ def minos(self, *parameters, cl=None, ncall=None):
             factor = chi2(1).ppf(cl)
 
         if not self._fmin:
-            raise RuntimeError(
-                "MINOS require function to be at the minimum. Run MIGRAD first."
+            # create a FunctionMinimum for MnMinos
+            fm = FunctionMinimum(
+                self._fcn, self._last_state, self._strategy, self._tolerance
             )
+        else:
+            fm = self._fmin._src
 
         ncall = 0 if ncall is None else int(ncall)
 
-        if not self._fmin.is_valid:
+        if not fm.is_valid:
             raise RuntimeError(
                 "Function minimum is not valid. Make sure MIGRAD converged."
             )
@@ -833,7 +833,7 @@ def minos(self, *parameters, cl=None, ncall=None):
                     pars.append(par)
 
         with TemporaryErrordef(self._fcn, factor):
-            minos = MnMinos(self._fcn, self._fmin._src, self.strategy)
+            minos = MnMinos(self._fcn, fm, self.strategy)
             for par in pars:
                 me = minos(self._var2pos[par], ncall, self._tolerance)
                 self._merrors[par] = mutil.MError(
@@ -854,8 +854,9 @@ def minos(self, *parameters, cl=None, ncall=None):
                     me.min,
                 )
 
-        self._fmin._nfcn = self.nfcn
-        self._fmin._ngrad = self.ngrad
+        if self._fmin:
+            self._fmin._nfcn = self.nfcn
+            self._fmin._ngrad = self.ngrad
 
         return self  # return self for method chaining and to autodisplay current state
 
@@ -1353,6 +1354,17 @@ def _make_covariance(self):
         else:
             self._covariance = None
 
+    def _migrad_edm_goal(self):
+        pr = MnMachinePrecision()
+        if self.precision is not None:
+            pr.eps = self.precision
+        # EDM goal
+        # - taken from the source code, see VariableMeticBuilder::Minimum and
+        #   ModularFunctionMinimizer::Minimize
+        # - goal is used to detect convergence but violations by 10x are also accepted;
+        #   see VariableMetricBuilder.cxx:425
+        return 2e-3 * max(self.tol * self.errordef, pr.eps2)
+
     def __repr__(self):
         s = []
         if self.fmin is not None:

src/machineprecision.cpp

@@ -18,6 +18,8 @@ using namespace ROOT::Minuit2;
 void bind_machineprecision(py::module m) {
   py::class_<MnMachinePrecision>(m, "MnMachinePrecision")
 
+      .def(py::init<>())
+
       .def_property("eps", &MnMachinePrecision::Eps, &MnMachinePrecision::SetPrecision)
       .def_property_readonly("eps2", &MnMachinePrecision::Eps2)
 

tests/test_core.py

@@ -197,43 +197,19 @@ def test_FunctionMinimum():
     st = MnUserParameterState()
     st.add("x", 0.01, 5)
     str = MnStrategy(1)
-    fm1 = FunctionMinimum(fcn, st, 0.1, str, 0.2)
+    fm1 = FunctionMinimum(fcn, st, str, 0.2)
     assert fm1.is_valid
     assert len(fm1.state) == 1
-    assert fm1.fval == 0.1
-    fm2 = FunctionMinimum(fcn, st, 0.1, str, 0)
+    assert fm1.fval == 10.0001
+    fm2 = FunctionMinimum(fcn, st, str, 0)
     assert not fm2.is_valid
 
 
-def test_MnUserTransformation_pickle():
-    tr = MnUserTransformation()
-
-    pkl = pickle.dumps(tr)
-    tr2 = pickle.loads(pkl)
-
-    assert len(tr2) == len(tr)
-
-
-def test_MinimumState_pickle():
-    st = MinimumState(3)
-
-    pkl = pickle.dumps(st)
-    st2 = pickle.loads(pkl)
-
-    assert st.vec == st2.vec
-    assert st.fval == st2.fval
-    assert st.edm == st2.edm
-    assert st.nfcn == st2.nfcn
-    assert st.is_valid == st2.is_valid
-    assert st.has_parameters == st2.has_parameters
-    assert st.has_covariance == st2.has_covariance
-
-
 def test_FunctionMinimum_pickle():
     st = MnUserParameterState()
     st.add("x", 1, 0.1)
     st.add("y", 2, 0.1, 1, 3)
-    fm = FunctionMinimum(FCN(fn, None, False, 1), st, 123, 1, 0.1)
+    fm = FunctionMinimum(FCN(fn, None, False, 1), st, 1, 0.1)
 
     pkl = pickle.dumps(fm)
     fm2 = pickle.loads(pkl)
@@ -252,3 +228,27 @@ def test_FunctionMinimum_pickle():
     assert fm.is_above_max_edm == fm2.is_above_max_edm
     assert fm.has_reached_call_limit == fm2.has_reached_call_limit
     assert fm.errordef == fm2.errordef
+
+
+def test_MnUserTransformation_pickle():
+    tr = MnUserTransformation()
+
+    pkl = pickle.dumps(tr)
+    tr2 = pickle.loads(pkl)
+
+    assert len(tr2) == len(tr)
+
+
+def test_MinimumState_pickle():
+    st = MinimumState(3)
+
+    pkl = pickle.dumps(st)
+    st2 = pickle.loads(pkl)
+
+    assert st.vec == st2.vec
+    assert st.fval == st2.fval
+    assert st.edm == st2.edm
+    assert st.nfcn == st2.nfcn
+    assert st.is_valid == st2.is_valid
+    assert st.has_parameters == st2.has_parameters
+    assert st.has_covariance == st2.has_covariance

tests/test_minuit.py

@@ -1459,3 +1459,12 @@ def test_minos_new_min():
     # ...but interval is correct
     assert m.merrors["x"].lower == approx(-0.9, abs=1e-2)
     assert m.merrors["x"].upper == approx(1.1, abs=1e-2)
+
+
+@pytest.mark.parametrize("x", (0.9, 1.0, 1.1))
+def test_minos_without_migrad(x):
+    m = Minuit(lsq(lambda x: (x - 1) ** 2), x)
+    m.minos()
+    me = m.merrors["x"]
+    assert x + me.lower == approx(0, abs=5e-3)
+    assert x + me.upper == approx(2, abs=5e-3)