Fix a bug in computation of Minos errors when a parameter is bounded (r…

…oot-project#5899) (root-project#6060)

* Fix computation of Minos errrors in Minuit2 when parameters are bounded.
A bug was introduced when the conditions "new minimum found" was added in Minos.
(ROOT-10854)

Add some fixes also the message reporting in Minuit2. Restore usage of ROOT message reporting system when MInuit2 is built within ROOT

* Add a protection in TMinuit to avoid a 0/0 division. Fix ROOT-10835
This can happen when parameter limit is 0 and parameter value at minimum is 0
and a Minos error analysis is performed for that parameter.

* Fix running Minos error when a new minimum is found.
Fix the handling of the situation of Minos error in Minuit2 when a new minimum is found.
Add a function in the Minimizer to get the status code of Minos to flag when a new minimum is found
Fix function returning minimum values in Minuit2Minimizer in case of multiple runs. The minimum values are now copied in the Mininuit2Minimizer local vector when the minimization is run

* Update Fitter::CalculateMinosError when a new minimum is found in Minos
Recompute the minos error for the previous run parameter in case a new minimum was found

* Use fabs instead of abs. Fix problem reported by CI builds

* small fixes in documentation as suggested by Stephan's review

* Revert changes in the Minuit2 CMAKE file for enabling Minuit2 error and warning messages.
These changes will be done in a separate commit and  a new PR.

math/mathcore/inc/Math/Minimizer.h

@@ -344,7 +344,6 @@ class Minimizer {
       minos error for variable i, return false if Minos failed or not supported
       and the lower and upper errors are returned in errLow and errUp
       An extra flag  specifies if only the lower (option=-1) or the upper (option=+1) error calculation is run
-      (This feature is not yet implemented)
    */
    virtual bool GetMinosError(unsigned int ivar , double & errLow, double & errUp, int option = 0) {
       MATH_ERROR_MSG("Minimizer::GetMinosError","Minos Error not implemented");
@@ -429,6 +428,9 @@ class Minimizer {
    /// status code of minimizer
    int Status() const { return fStatus; }
 
+   /// status code of Minos (to be re-implemented by the minimizers supporting Minos)
+   virtual int MinosStatus() const { return -1; }
+
    /// return the statistical scale used for calculate the error
    /// is typically 1 for Chi2 and 0.5 for likelihood minimization
    double ErrorDef() const { return fOptions.ErrorDef(); }

math/mathcore/src/Fitter.cxx

@@ -764,15 +764,34 @@ bool Fitter::CalculateMinosErrors() {
    const std::vector<unsigned int> & ipars = fConfig.MinosParams();
    unsigned int n = (ipars.size() > 0) ? ipars.size() : fResult->Parameters().size();
    bool ok = false;
-   for (unsigned int i = 0; i < n; ++i) {
-      double elow, eup;
-      unsigned int index = (ipars.size() > 0) ? ipars[i] : i;
-      bool ret = fMinimizer->GetMinosError(index, elow, eup);
-      if (ret) fResult->SetMinosError(index, elow, eup);
-      ok |= ret;
+
+   int iparNewMin = 0; 
+   int iparMax = n; 
+   int iter = 0; 
+   // rerun minos for the parameters run before a new Minimum has been found
+   do {
+      if (iparNewMin > 0) 
+         MATH_INFO_MSG("Fitter::CalculateMinosErrors","Run again Minos for some parameters because a new Minimum has been found");
+      iparNewMin = 0; 
+      for (int i = 0; i < iparMax; ++i) {
+         double elow, eup;
+         unsigned int index = (ipars.size() > 0) ? ipars[i] : i;
+         bool ret = fMinimizer->GetMinosError(index, elow, eup);
+         // flags case when a new minimum has been found
+         if ((fMinimizer->MinosStatus() & 8) != 0) {
+            iparNewMin = i;
+         }
+         if (ret)
+            fResult->SetMinosError(index, elow, eup);
+         ok |= ret;
+      }
+
+      iparMax = iparNewMin;
+      iter++;  // to avoid infinite looping 
    }
+   while( iparNewMin > 0 && iter < 10);
    if (!ok) {
-       MATH_ERROR_MSG("Fitter::CalculateMinosErrors","Minos error calculation failed for all parameters");
+       MATH_ERROR_MSG("Fitter::CalculateMinosErrors","Minos error calculation failed for all the selected parameters");
    }
 
    // re-give a minimizer instance in case it has been changed

math/minuit/inc/TMinuitMinimizer.h

@@ -192,6 +192,12 @@ class TMinuitMinimizer  : public ROOT::Math::Minimizer {
    /// minos error for variable i, return false if Minos failed
    virtual bool GetMinosError(unsigned int i, double & errLow, double & errUp, int = 0);
 
+   /// minos status code of last Minos run
+   /// minos status = -1   : Minos is not run
+   ///              =  0   : last MINOS run was succesfull
+   ///              >  0   : some problems encountered when running MINOS
+   virtual int MinosStatus() const { return fMinosStatus; }
+
    /**
       perform a full calculation of the Hessian matrix for error calculation
     */
@@ -267,6 +273,7 @@ class TMinuitMinimizer  : public ROOT::Math::Minimizer {
    bool fUsed;
    bool fMinosRun;
    unsigned int fDim;
+   int fMinosStatus = -1;         // Minos status code
    std::vector<double> fParams;   // vector of output values
    std::vector<double> fErrors;   // vector of output errors
    std::vector<double> fCovar;    // vector storing the covariance matrix

math/minuit/src/TMinuit.cxx

@@ -5555,7 +5555,9 @@ void TMinuit::mnmnot(Int_t ilax, Int_t ilax2, Double_t &val2pl, Double_t &val2mi
       fU[ilax-1] = TMath::Max(fU[ilax-1],fAlim[ilax-1]);
       delu = fU[ilax-1] - ut;
 //        stop if already at limit with negligible step size
-      if (TMath::Abs(delu) / (TMath::Abs(ut) + TMath::Abs(fU[ilax-1])) < fEpsmac) goto L440;
+//        add also a check if both numerator and denominarot are not zero (ROOT-10835)(LM)
+      if ( (delu == 0 && ut == 0) ||
+         (TMath::Abs(delu) / (TMath::Abs(ut) + TMath::Abs(fU[ilax-1])) < fEpsmac)) goto L440;
       fac = delu / fMNOTw[it-1];
       for (i = 1; i <= fNpar; ++i) {
          fX[i-1] = fXt[i-1] + fac*fMNOTxdev[i-1];

math/minuit/src/TMinuitMinimizer.cxx

@@ -804,7 +804,7 @@ bool TMinuitMinimizer::GetMinosError(unsigned int i, double & errLow, double & e
    }
 
 
-   // syntax of MINOS is MINOS [maxcalls] [parno]
+   // syntax of MINOS command is:  MINOS [maxcalls] [parno]
    // if parno = 0 all parameters are done
    arglist[0] = MaxFunctionCalls();
    arglist[1] = i+1;  // par number starts from 1 in TMInuit
@@ -815,7 +815,7 @@ bool TMinuitMinimizer::GetMinosError(unsigned int i, double & errLow, double & e
    // check also the status from fCstatu
    if (isValid && fMinuit->fCstatu != "SUCCESSFUL") {
       if (fMinuit->fCstatu == "FAILURE" ) {
-         // in this case MINOS failed on all prameter, so it is not valid !
+         // in this case MINOS failed on all parameters, so it is not valid !
          ierr = 5;
          isValid = false;
       }
@@ -824,6 +824,7 @@ bool TMinuitMinimizer::GetMinosError(unsigned int i, double & errLow, double & e
    }
 
    fStatus += 10*ierr;
+   fMinosStatus = ierr; 
 
    fMinosRun = true;
 

math/minuit2/inc/Minuit2/MinosError.h

@@ -26,19 +26,19 @@ class MinosError {
 
 public:
 
-   MinosError() : fParameter(0), fMinValue(0.), fUpper(MnCross()), fLower(MnCross()) {}
+   MinosError() : fParameter(0), fMinParValue(0.), fUpper(MnCross()), fLower(MnCross()) {}
 
-   MinosError(unsigned int par, double min, const MnCross& low, const MnCross& up) : fParameter(par), fMinValue(min), fUpper(up), fLower(low) {}
+   MinosError(unsigned int par, double value, const MnCross& low, const MnCross& up) : fParameter(par), fMinParValue(value), fUpper(up), fLower(low) {}
 
    ~MinosError() {}
 
-   MinosError(const MinosError& err) : fParameter(err.fParameter), fMinValue(err.fMinValue), fUpper(err.fUpper),  fLower(err.fLower) {}
+   MinosError(const MinosError& err) : fParameter(err.fParameter), fMinParValue(err.fMinParValue), fUpper(err.fUpper),  fLower(err.fLower) {}
 
    MinosError& operator=(const MinosError& ) = default;
 
    MinosError& operator()(const MinosError& err) {
       fParameter = err.fParameter;
-      fMinValue = err.fMinValue;
+      fMinParValue = err.fMinParValue;
       fUpper = err.fUpper;
       fLower = err.fLower;
       return *this;
@@ -48,12 +48,20 @@ class MinosError {
       return std::pair<double,double>(Lower(), Upper());
    }
    double Lower() const {
-      if ( AtLowerLimit() ) return LowerState().Parameter( Parameter() ).LowerLimit() -  fMinValue;
-      return -1.*LowerState().Error(Parameter())*(1. + fLower.Value());
+      if (AtLowerLimit())
+         return LowerState().Parameter(Parameter()).LowerLimit() - fMinParValue;
+      if (LowerValid() )
+         return -1. * LowerState().Error(Parameter()) * (1. + fLower.Value());
+      // return Hessian Error in case is invalid
+      return - LowerState().Error(Parameter());
    }
    double Upper() const {
-      if ( AtUpperLimit() ) return UpperState().Parameter( Parameter() ).UpperLimit() -  fMinValue;
-      return UpperState().Error(Parameter())*(1. + fUpper.Value());
+      if (AtUpperLimit())
+         return UpperState().Parameter(Parameter()).UpperLimit() - fMinParValue;
+      if (UpperValid())
+         return UpperState().Error(Parameter()) * (1. + fUpper.Value());
+      // return Hessian Error in case is invalid
+      return  UpperState().Error(Parameter());
    }
    unsigned int Parameter() const {return fParameter;}
    const MnUserParameterState& LowerState() const {return fLower.State();}
@@ -68,12 +76,13 @@ class MinosError {
    bool LowerNewMin() const {return fLower.NewMinimum();}
    bool UpperNewMin() const {return fUpper.NewMinimum();}
    unsigned int NFcn() const {return fUpper.NFcn() + fLower.NFcn();}
-   double Min() const {return fMinValue;}
+   // return parameter value at the minimum
+   double Min() const {return fMinParValue;}
 
 private:
 
    unsigned int fParameter;
-   double fMinValue;
+   double fMinParValue;
    MnCross fUpper;
    MnCross fLower;
 };

math/minuit2/inc/Minuit2/Minuit2Minimizer.h

@@ -157,7 +157,7 @@ class Minuit2Minimizer : public ROOT::Math::Minimizer {
    virtual double Edm() const { return fState.Edm(); }
 
    /// return  pointer to X values at the minimum
-   virtual const double *  X() const;
+   virtual const double *  X() const { return  &fValues.front();}
 
    /// return pointer to gradient values at the minimum
    virtual const double *  MinGradient() const { return 0; } // not available in Minuit2
@@ -237,16 +237,21 @@ class Minuit2Minimizer : public ROOT::Math::Minimizer {
       get the minos error for parameter i, return false if Minos failed
       A minimizaiton must be performed befre, return false if no minimization has been done
       In case of Minos failed the status error is updated as following
-      status += 10 * minosStatus where the minos status is:
-       status = 1    : maximum number of function calls exceeded when running for lower error
-       status = 2    : maximum number of function calls exceeded when running for upper error
-       status = 3    : new minimum found when running for lower error
-       status = 4    : new minimum found when running for upper error
-       status = 5    : any other failure
-
+      status += 10 * minosStatus. 
+      The Minos status of last Minos run can also be retrieved by calling MinosStatus()
    */
    virtual bool GetMinosError(unsigned int i, double & errLow, double & errUp, int = 0);
 
+   /** 
+      MINOS status code of last Minos run
+       `status & 1 > 0`  : invalid lower error
+       `status & 2 > 0`  : invalid upper error
+       `status & 4 > 0`  : invalid because maximum number of function calls exceeded
+       `status & 8 > 0`  : a new minimum has been found
+       `status & 16 > 0` : error is truncated because parameter is at lower/upper limit
+   */
+   virtual int MinosStatus() const { return fMinosStatus; }
+
    /**
       scan a parameter i around the minimum. A minimization must have been done before,
       return false if it is not the case
@@ -305,10 +310,14 @@ class Minuit2Minimizer : public ROOT::Math::Minimizer {
    /// examine the minimum result
    bool ExamineMinimum(const ROOT::Minuit2::FunctionMinimum & min);
 
+   // internal function to compute Minos errors
+   int RunMinosError(unsigned int i, double & errLow, double & errUp, int runopt);
+
 private:
 
    unsigned int fDim;       // dimension of the function to be minimized
    bool fUseFumili;
+   int fMinosStatus = -1;  // Minos status code
 
    ROOT::Minuit2::MnUserParameterState fState;
    // std::vector<ROOT::Minuit2::MinosError> fMinosErrors;