[minuit2] Fixes in Minuit2 for finding contour points

math/minuit2/inc/Minuit2/MinosError.h

@@ -55,17 +55,30 @@ class MinosError {
    {
       if (AtLowerLimit())
          return LowerState().Parameter(Parameter()).LowerLimit() - fMinParValue;
-      if (LowerValid())
-         return -1. * LowerState().Error(Parameter()) * (1. + fLower.Value());
+      if (LowerValid()) {
+         // Minos error is   value - error - aopt * error  where aopt is MnCross.Value()
+         // If value - error is below the limit, error must be truncated at limit
+         double err = LowerState().Error(Parameter());
+         // error is truncated if over the limit
+         if (LowerState().Parameter(Parameter()).HasLowerLimit())
+            err = std::min(err, fMinParValue - LowerState().Parameter(Parameter()).LowerLimit());
+         return -1. * err * (1. + fLower.Value());
+      }
       // return Hessian Error in case is invalid
       return -LowerState().Error(Parameter());
    }
    double Upper() const
    {
       if (AtUpperLimit())
          return UpperState().Parameter(Parameter()).UpperLimit() - fMinParValue;
-      if (UpperValid())
-         return UpperState().Error(Parameter()) * (1. + fUpper.Value());
+      if (UpperValid()) {
+         // Minos error is   value + error + aopt * error  where aopt is MnCross.Value()
+         // If value + error is over the limit,  err must be truncated at limit
+         double err = UpperState().Error(Parameter());
+         if (UpperState().Parameter(Parameter()).HasUpperLimit())
+            err = std::min(err, UpperState().Parameter(Parameter()).UpperLimit() - fMinParValue);
+         return err * (1. + fUpper.Value());
+      }
       // return Hessian Error in case is invalid
       return UpperState().Error(Parameter());
    }

math/minuit2/src/MnContours.cxx

@@ -40,75 +40,97 @@ ContoursError MnContours::Contour(unsigned int px, unsigned int py, unsigned int
    unsigned int nfcn = 0;
 
    MnPrint print("MnContours");
+   print.Debug("MnContours: finding ",npoints," contours points for ",px,py," at level ",fFCN.Up()," from value ",fMinimum.Fval());
 
    std::vector<std::pair<double, double>> result;
    result.reserve(npoints);
    std::vector<MnUserParameterState> states;
    //   double edmmax = 0.5*0.05*fFCN.Up()*1.e-3;
 
    // double toler = 0.05;
-   double toler = 0.1; // use same default value as in Minos
-
-   // get first four points
-   //   std::cout<<"MnContours: get first 4 params."<<std::endl;
+   double toler = 0.1; // use same defaut value as in Minos
+
+   // get first four points running Minos separately on the two parameters
+   // and then finding the corresponding minimum in the other
+   // P1( exlow, ymin1)  where ymin1 is the parameter value (y) at the minimum of f when x is fixed to exlow
+   // P2(xmin1, eylow)  where  xmin1 is the  the parameter value (x) at the minimum of f when y is fixed to eylow
+   // P3(exup, ymin2)
+   // P4(xmin2, eyup)
    MnMinos minos(fFCN, fMinimum, fStrategy);
 
    double valx = fMinimum.UserState().Value(px);
    double valy = fMinimum.UserState().Value(py);
 
-   MinosError meX = minos.Minos(px);
-   nfcn += meX.NFcn();
-   if (!meX.IsValid()) {
+   print.Debug("Run Minos to find first 4 contour points. Current minimum is : ",valx,valy);
+
+   MinosError mnex = minos.Minos(px);
+   nfcn += mnex.NFcn();
+   if (!mnex.IsValid()) {
       print.Error("unable to find first two points");
-      return ContoursError(px, py, result, meX, meX, nfcn);
+      return ContoursError(px, py, result, mnex, mnex, nfcn);
    }
-   std::pair<double, double> ex = meX();
+   std::pair<double, double> ex = mnex();
+
+   print.Debug("Minos error for p0:  ",ex.first,ex.second);
 
-   MinosError meY = minos.Minos(py);
-   nfcn += meY.NFcn();
-   if (!meY.IsValid()) {
+   MinosError mney = minos.Minos(py);
+   nfcn += mney.NFcn();
+   if (!mney.IsValid()) {
       print.Error("unable to find second two points");
-      return ContoursError(px, py, result, meX, meY, nfcn);
+      return ContoursError(px, py, result, mnex, mney, nfcn);
    }
-   std::pair<double, double> ey = meY();
+   std::pair<double, double> ey = mney();
 
-   MnMigrad migrad(fFCN, fMinimum.UserState(), MnStrategy(std::max(0, int(fStrategy.Strategy() - 1))));
+   print.Debug("Minos error for p0:  ",ey.first,ey.second);
 
-   migrad.Fix(px);
-   migrad.SetValue(px, valx + ex.second);
-   FunctionMinimum exy_up = migrad();
-   nfcn += exy_up.NFcn();
-   if (!exy_up.IsValid()) {
-      print.Error("unable to find Upper y Value for x Parameter", px);
-      return ContoursError(px, py, result, meX, meY, nfcn);
-   }
+   // if Minos is not at limits we can use migrad to find the other corresponding point coordinate
+   MnMigrad migrad0(fFCN, fMinimum.UserState(), MnStrategy(std::max(0, int(fStrategy.Strategy() - 1))));
 
-   migrad.SetValue(px, valx + ex.first);
-   FunctionMinimum exy_lo = migrad();
+
+   // start from minimizing in p1 and fixing p0 to Minos value
+   migrad0.Fix(px);
+   migrad0.SetValue(px, valx + ex.first);
+   FunctionMinimum exy_lo = migrad0();
    nfcn += exy_lo.NFcn();
    if (!exy_lo.IsValid()) {
       print.Error("unable to find Lower y Value for x Parameter", px);
-      return ContoursError(px, py, result, meX, meY, nfcn);
+      return ContoursError(px, py, result, mnex, mney, nfcn);
    }
 
+   print.Debug("Minimum p1 found for p0 set to ",migrad0.Value(px)," is ",exy_lo.UserState().Value(py),"fcn = ",exy_lo.Fval());
+
+   migrad0.SetValue(px, valx + ex.second);
+   FunctionMinimum exy_up = migrad0();
+   nfcn += exy_up.NFcn();
+   if (!exy_up.IsValid()) {
+      print.Error("unable to find Upper y Value for x Parameter", px);
+      return ContoursError(px, py, result, mnex, mney, nfcn);
+   }
+   print.Debug("Minimum p1 found for p0 set to ",migrad0.Value(px)," is ",exy_up.UserState().Value(py),"fcn = ",exy_up.Fval());
+
+
    MnMigrad migrad1(fFCN, fMinimum.UserState(), MnStrategy(std::max(0, int(fStrategy.Strategy() - 1))));
    migrad1.Fix(py);
    migrad1.SetValue(py, valy + ey.second);
    FunctionMinimum eyx_up = migrad1();
    nfcn += eyx_up.NFcn();
    if (!eyx_up.IsValid()) {
       print.Error("unable to find Upper x Value for y Parameter", py);
-      return ContoursError(px, py, result, meX, meY, nfcn);
+      return ContoursError(px, py, result, mnex, mney, nfcn);
    }
+   print.Debug("Minimum p0 found for p1 set to ",migrad1.Value(py)," is ",eyx_up.UserState().Value(px),"fcn = ",eyx_up.Fval());
 
    migrad1.SetValue(py, valy + ey.first);
    FunctionMinimum eyx_lo = migrad1();
    nfcn += eyx_lo.NFcn();
    if (!eyx_lo.IsValid()) {
       print.Error("unable to find Lower x Value for y Parameter", py);
-      return ContoursError(px, py, result, meX, meY, nfcn);
+      return ContoursError(px, py, result, mnex, mney, nfcn);
    }
 
+   print.Debug("Minimum p0 found for p1 set to ",migrad1.Value(py)," is ",eyx_lo.UserState().Value(px),"fcn = ",eyx_lo.Fval());
+
+
    double scalx = 1. / (ex.second - ex.first);
    double scaly = 1. / (ey.second - ey.first);
 
@@ -119,7 +141,7 @@ ContoursError MnContours::Contour(unsigned int px, unsigned int py, unsigned int
 
    MnUserParameterState upar = fMinimum.UserState();
 
-   print.Info("List of found points", '\n', "  Parameter x is", upar.Name(px), '\n', "  Parameter y is", upar.Name(py),
+   print.Debug("List of first 4 found contour points", '\n', "  Parameter x is", upar.Name(px), '\n', "  Parameter y is", upar.Name(py),
               '\n', result[0], '\n', result[1], '\n', result[2], '\n', result[3]);
 
    upar.Fix(px);
@@ -130,15 +152,17 @@ ContoursError MnContours::Contour(unsigned int px, unsigned int py, unsigned int
    par[1] = py;
    MnFunctionCross cross(fFCN, upar, fMinimum.Fval(), fStrategy);
 
+   // find the remaining points of the contour
    for (unsigned int i = 4; i < npoints; i++) {
 
+      //  find the two neighbouring points with largest separation
       auto idist1 = result.end() - 1;
       auto idist2 = result.begin();
       double dx = idist1->first - (idist2)->first;
       double dy = idist1->second - (idist2)->second;
       double bigdis = scalx * scalx * dx * dx + scaly * scaly * dy * dy;
 
-      for (auto ipair = result.begin(); ipair != result.end() - 1; ++ipair) {
+      for (auto  ipair = result.begin(); ipair != result.end() - 1; ++ipair) {
          double distx = ipair->first - (ipair + 1)->first;
          double disty = ipair->second - (ipair + 1)->second;
          double dist = scalx * scalx * distx * distx + scaly * scaly * disty * disty;
@@ -157,11 +181,17 @@ ContoursError MnContours::Contour(unsigned int px, unsigned int py, unsigned int
 
       if (nfcn > maxcalls) {
          print.Error("maximum number of function calls exhausted");
-         return ContoursError(px, py, result, meX, meY, nfcn);
+         return ContoursError(px, py, result, mnex, mney, nfcn);
       }
 
+      print.Debug("Find new contour point between points with max sep:  (",idist1->first,", ",idist1->second,") and (",
+                                                                           idist2->first,", ",idist2->second,")  with weights ",a1,a2);
+      // find next point between the found 2 with max separation
+      // start from point situated at the middle (a1,a2=0.5)
+      // and direction
       double xmidcr = a1 * idist1->first + a2 * (idist2)->first;
       double ymidcr = a1 * idist1->second + a2 * (idist2)->second;
+      // direction is the perpendicular one
       double xdir = (idist2)->second - idist1->second;
       double ydir = idist1->first - (idist2)->first;
       double scalfac = sca * std::max(std::fabs(xdir * scalx), std::fabs(ydir * scaly));
@@ -177,30 +207,39 @@ ContoursError MnContours::Contour(unsigned int px, unsigned int py, unsigned int
       MnCross opt = cross(par, pmid, pdir, toler, maxcalls);
       nfcn += opt.NFcn();
       if (!opt.IsValid()) {
-         //       if(a1 > 0.5) {
-         if (sca < 0.) {
+         if(a1 > 0.5) {
+         // LM 20/10/23 : remove switch of direction and look instead closer (this is what is done in TMinuit)
+         // should we try again closer to P2 (e.g. a1=0.25, a2 = 0.75) if failing?
+         //if (sca < 0.) {
             print.Error("unable to find point on Contour", i + 1, '\n', "found only", i, "points");
-            return ContoursError(px, py, result, meX, meY, nfcn);
+            return ContoursError(px, py, result, mnex, mney, nfcn);
          }
-         //       a1 = 0.75;
-         //       a2 = 0.25;
-         //       std::cout<<"*****switch direction"<<std::endl;
-         sca = -1.;
+         a1 = 0.75;
+         a2 = 0.25;
+         print.Debug("Unable to find point, try closer to p1 with weight values",a1,a2);
+         //std::cout<<"*****switch direction"<<std::endl;
+         //sca = -1.;
          goto L300;
       }
       double aopt = opt.Value();
+      int pos = result.size();
       if (idist2 == result.begin()) {
          result.emplace_back(xmidcr + (aopt)*xdircr, ymidcr + (aopt)*ydircr);
          print.Info(result.back());
       } else {
          print.Info(*idist2);
-         result.insert(idist2, {xmidcr + (aopt)*xdircr, ymidcr + (aopt)*ydircr});
+         auto itr = result.insert(idist2, {xmidcr + (aopt)*xdircr, ymidcr + (aopt)*ydircr});
+         pos = std::distance(result.begin(),itr);
       }
+      print.Info("Found new contour point - pos: ",pos,result[pos]);
    }
 
    print.Info("Number of contour points =", result.size());
+   print.Debug("List of contour points");
+   for (size_t i = 0; i < result.size(); i++)
+      print.Debug("point ",i,result[i]);
 
-   return ContoursError(px, py, result, meX, meY, nfcn);
+   return ContoursError(px, py, result, mnex, mney, nfcn);
 }
 
 } // namespace Minuit2

math/minuit2/src/MnSeedGenerator.cxx

@@ -69,8 +69,8 @@ operator()(const MnFcn &fcn, const GradientCalculator &gc, const MnUserParameter
       dcovar = 0.;
    } else {
       for (unsigned int i = 0; i < n; i++)
-         mat(i, i) = (std::fabs(dgrad.G2()(i)) > prec.Eps2() ? 1. / dgrad.G2()(i) :
-          (dgrad.G2()(i) >= 0) ? 1./prec.Eps2() : -1./prec.Eps2());
+        // if G2 is small better using an arbitrary value (e.g. 1)
+         mat(i, i) = std::fabs(dgrad.G2()(i)) > prec.Eps() ? 1. / dgrad.G2()(i) : 1.0;
    }
    MinimumError err(mat, dcovar);
 
@@ -177,6 +177,7 @@ MinimumSeed MnSeedGenerator::operator()(const MnFcn &fcn, const AnalyticalGradie
    // do this only when we have not computed the Hessian or always ?
    if (!computedHessian) {
       // check if minimum state has covariance - if not use computed G2
+      // should maybe this an option, sometimes is not good to re-use existing covariance
       if (st.HasCovariance()) {
          print.Info("Using existing covariance matrix");
          for (unsigned int i = 0; i < n; i++)
@@ -185,10 +186,9 @@ MinimumSeed MnSeedGenerator::operator()(const MnFcn &fcn, const AnalyticalGradie
          dcovar = 0.;
       } else {
          for (unsigned int i = 0; i < n; i++) {
-            // should not use a cut-off here like 1./prec.Eps()
-            mat(i, i) = (std::fabs(grad.G2()(i)) > prec.Eps2() ? 1. / grad.G2()(i)
-                         : (grad.G2()(i) >= 0)                 ? 1. / prec.Eps2()
-                                                               : -1. / prec.Eps2());
+            // if G2 is very small, better using an arbitrary value (e.g. 1.)
+            mat(i, i) = std::fabs(grad.G2()(i)) > prec.Eps() ? 1. / grad.G2()(i)
+                        : 1.0;
          }
          dcovar = 1.;
       }

math/minuit2/src/NegativeG2LineSearch.cxx

@@ -121,8 +121,8 @@ MinimumState NegativeG2LineSearch::operator()(const MnFcn &fcn, const MinimumSta
    print.Debug("Approximate new covariance after NegativeG2LS using only G2");
    MnAlgebraicSymMatrix mat(n);
    for (unsigned int i = 0; i < n; i++) {
-      mat(i, i) = (std::fabs(dgrad.G2()(i)) > prec.Eps2() ? 1. / dgrad.G2()(i) :
-                  (dgrad.G2()(i) >= 0) ? 1./prec.Eps2() : -1./prec.Eps2());
+      mat(i, i) = std::fabs(dgrad.G2()(i)) > prec.Eps() ? 1. / dgrad.G2()(i) :
+                  1; // use an arbitrary value (e.g. 1)
    }
 
    MinimumError err(mat, 1.);

tutorials/roofit/rf601_intminuit.C

@@ -54,9 +54,6 @@ void rf601_intminuit()
    // Create MINUIT interface object
    RooMinimizer m(*nll);
 
-   // The contour algorithm in Minuit2 does not work for this tutorial
-   m.setMinimizerType("Minuit");
-
    // Activate verbose logging of MINUIT parameter space stepping
    m.setVerbose(true);
 

tutorials/roofit/rf601_intminuit.py

@@ -46,9 +46,6 @@
 # Create MINUIT interface object
 m = ROOT.RooMinimizer(nll)
 
-# The contour algorithm in Minuit2 does not work for this tutorial
-m.setMinimizerType("Minuit")
-
 # Activate verbose logging of MINUIT parameter space stepping
 m.setVerbose(True)