[math][minuit2] Apply small fixes and improvements on new Fumili version

Fix also fit option V, VV and VVV

hist/hist/src/HFitImpl.cxx

@@ -787,10 +787,12 @@ void ROOT::Fit::FitOptionsMake(EFitObjectType type, const char *option, Foption_
          fitOption.User = 0;
       }
    }
-   if (opt.Contains("Q")) fitOption.Quiet   = 1;
-   if (opt.Contains("V")) {fitOption.Verbose = 1; fitOption.Quiet   = 0;}
-   if (opt.Contains("VV")) {fitOption.Verbose = 2; }
-    if (opt.Contains("VVV") || opt.Contains("DEBUG")) {fitOption.Verbose = 3; }
+
+   // in case of Q and V options V has precedence
+   if (opt.Contains("VVV") || opt.Contains("DEBUG")) { fitOption.Verbose = 3; }
+   else if (opt.Contains("VV")) {fitOption.Verbose = 2; }
+   else if (opt.Contains("V")) {fitOption.Verbose = 1; }
+   else if (opt.Contains("Q")) {fitOption.Quiet   = 1; }
 
 
    if (opt.Contains("E")) fitOption.Errors  = 1;

math/mathmore/src/GSLNLSMinimizer.cxx

@@ -275,13 +275,8 @@ GSLNLSMinimizer::GSLNLSMinimizer( int type )
 }
 
 GSLNLSMinimizer::~GSLNLSMinimizer () {
-<<<<<<< HEAD
-   assert(fGSLMultiFit != nullptr);
-   delete fGSLMultiFit;
-=======
    if (fGSLMultiFit) delete fGSLMultiFit;
    if (fGSLMultiFit2) delete fGSLMultiFit2;
->>>>>>> c1c95febd88 ([math][mathmore] Add new fitter algorithm from GSL under GSLMultiFit)
 }
 
 
@@ -466,13 +461,8 @@ bool GSLNLSMinimizer::DoMinimize(const Func & fitFunc, FitterType * fitter) {
    }
 
    // save state with values and function value
-<<<<<<< HEAD
-   const double * x = fGSLMultiFit->X();
-   if (x == nullptr) return false;
-=======
    const double * x = fitter->X();
-   if (x == 0) return false;
->>>>>>> c1c95febd88 ([math][mathmore] Add new fitter algorithm from GSL under GSLMultiFit)
+   if (x == nullptr) return false;
    // apply transformation outside SetFinalValues(..)
    // because trFunc is not a MinimTransformFunction but a FitTransFormFunction
    if (trFunc)  x = trFunc->Transformation(x);

math/minuit2/inc/Minuit2/FumiliFCNAdapter.h

@@ -54,13 +54,9 @@ class FumiliFCNAdapter : public FumiliFCNBase {
 
    void SetErrorDef(double up) override { fUp = up; }
 
-   // virtual std::vector<double> Gradient(std::vector<double> const &) const;
-
-   // forward interface
-   // virtual double operator()(int npar, double* params,int iflag = 4) const;
 
    /**
-       evaluate gradient hessian and function value needed by fumili
+       evaluate gradient hessian and function value needed by Fumili
      */
    void EvaluateAll(std::vector<double> const &v) override;
 
@@ -130,7 +126,7 @@ void FumiliFCNAdapter<Function>::EvaluateAll(std::vector<double> const &v)
    } else if (fFunc.Type() == Function::kPoissonLikelihood) {
       print.Debug("Poisson Likelihood FCN: Evaluate gradient and Hessian");
       // for Poisson need Hessian computed in DataElement since one needs the bin expected value ad bin observed value
-     for (unsigned int i = 0; i < ndata; ++i) {
+      for (unsigned int i = 0; i < ndata; ++i) {
          // calculate data element and gradient
          fFunc.DataElement(&v.front(), i, gf.data(), h.data());
          for (size_t j = 0; j < npar; ++j) {
@@ -140,7 +136,7 @@ void FumiliFCNAdapter<Function>::EvaluateAll(std::vector<double> const &v)
                hess[idx] += h[idx];
             }
          }
-     }
+      }
    } else {
       print.Error("Type of fit method is not supported, it must be chi2 or log-likelihood or Poisson Likelihood");
    }

math/minuit2/inc/Minuit2/FumiliFCNBase.h

@@ -70,7 +70,6 @@ class FumiliFCNBase : public FCNBase {
    {
    }
 
-   //   FumiliFCNBase(const ParametricFunction& modelFCN) { fModelFunction = &modelFCN; }
 
    ~FumiliFCNBase() override {}
 

math/minuit2/src/FumiliBuilder.cxx

@@ -48,16 +48,19 @@ FunctionMinimum FumiliBuilder::Minimum(const MnFcn &fcn, const GradientCalculato
    print.Debug("Convergence when edm <", edmval);
 
    if (seed.Parameters().Vec().size() == 0) {
+      print.Warn("No variable parameters are defined! - Return current function value ");
       return FunctionMinimum(seed, fcn.Up());
    }
 
-   //   double edm = Estimator().Estimate(seed.Gradient(), seed.Error());
-   double edm = seed.State().Edm();
+   // estimate initial edm value
+   double edm = Estimator().Estimate(seed.Gradient(), seed.Error());
+   print.Debug("initial edm is ", edm);
+   //double edm = seed.State().Edm();
 
    FunctionMinimum min(seed, fcn.Up());
 
    if (edm < 0.) {
-      print.Warn("Initial matrix not pos.def.");
+      print.Error("Initial matrix not positive defined, edm = ",edm,"\nExit minimization ");
       return min;
    }
 

math/minuit2/src/FumiliErrorUpdator.cxx

@@ -49,11 +49,13 @@ MinimumError FumiliErrorUpdator::Update(const MinimumState &s0, const MinimumPar
    // we apply also the Marquard lambda factor to increase weight of diagonal term
    // as suggester in Numerical Receipt for Marquard method
 
+   MnPrint print("FumiliErrorUpdator");
+   print.Debug("Compute covariance matrix using Fumili method");
+
    // need to downcast to FumiliGradientCalculator
    FumiliGradientCalculator *fgc = dynamic_cast<FumiliGradientCalculator *>(const_cast<GradientCalculator *>(&gc));
    assert(fgc != nullptr);
 
-   MnPrint print("FumiliErrorUpdator");
 
    // get Hessian from Gradient calculator
 
@@ -83,6 +85,9 @@ MinimumError FumiliErrorUpdator::Update(const MinimumState &s0, const MinimumPar
       for (unsigned int i = 0; i < cov.Nrow(); i++) {
          cov(i, i) = 1. / cov(i, i);
       }
+
+      // shiould we return the Hessian in addition to cov in this case?
+      return MinimumError(cov, MinimumError::MnInvertFailed);
    }
 
    double dcov = -1;

math/minuit2/src/FumiliGradientCalculator.cxx

@@ -39,13 +39,14 @@ FunctionGradient FumiliGradientCalculator::operator()(const MinimumParameters &p
    // Need to apply internal to external for parameters and the external to int transformation
    // for the return gradient and Hessian
 
+   MnPrint print("FumiliGradientCalculator");
+
    int nvar = par.Vec().size();
    std::vector<double> extParam = fTransformation(par.Vec());
 
    // eval Gradient
    FumiliFCNBase &fcn = const_cast<FumiliFCNBase &>(fFcn);
 
-   // evaluate gradient and Hessian
    fcn.EvaluateAll(extParam);
 
    MnAlgebraicVector v(nvar);
@@ -70,14 +71,16 @@ FunctionGradient FumiliGradientCalculator::operator()(const MinimumParameters &p
       }
    }
 
-   MnPrint print("FumiliGradientCalculator");
+
    print.Debug([&](std::ostream &os) {
       // compare Fumili with Minuit gradient
+      os << "Comparison of Fumili Gradient and standard (numerical) Minuit Gradient (done only when debugging enabled)" << std::endl;
       int plevel = MnPrint::SetGlobalLevel(MnPrint::GlobalLevel()-1);
       Numerical2PGradientCalculator gc(MnUserFcn(fFcn, fTransformation), fTransformation, MnStrategy(1));
       FunctionGradient grd2 = gc(par);
-      os << "Fumili Gradient:" << v << "\nMinuit Gradient" << grd2.Vec();
-      os << "\nFumili Hessian:  " << h << std::endl;
+      os << "Fumili Gradient:" << v << std::endl;
+      os << "Minuit Gradient" << grd2.Vec() << std::endl;
+      os << "Fumili Hessian:  " << h << std::endl;
       os << "Numerical g2 " << grd2.G2() << std::endl;
       MnPrint::SetGlobalLevel(plevel);
    });

math/minuit2/src/FumiliMinimizer.cxx

@@ -41,19 +41,15 @@ FunctionMinimum FumiliMinimizer::Minimize(const FCNBase &fcn, const MnUserParame
    // The FCNBase passed must be a FumiliFCNBase type otherwise method will fail !
 
    MnUserFcn mfcn(fcn, st.Trafo());
-   //  Numerical2PGradientCalculator gc(mfcn, st.Trafo(), strategy);
+
 
    unsigned int npar = st.VariableParameters();
    if (maxfcn == 0)
       maxfcn = 200 + 100 * npar + 5 * npar * npar;
    // FUMILI needs much less function calls
    maxfcn = int(0.1 * maxfcn);
 
-   // MinimumSeed mnseeds = SeedGenerator()(mfcn, gc, st, strategy);
-
-   // std::cout << "FCN type " << typeid(&fcn).Name() << std::endl;
-
-   // Minimize using Fumili. Case of interface is a function with gradient.
+   // Minimize using Fumili - function interface must be a FumiliFCNBase type
    FumiliFCNBase *fumiliFcn = dynamic_cast<FumiliFCNBase *>(const_cast<FCNBase *>(&fcn));
    if (!fumiliFcn) {
       print.Error("Wrong FCN type; try to use default minimizer");
@@ -77,7 +73,8 @@ FunctionMinimum FumiliMinimizer::Minimize(const FCNBase &fcn, const MnUserParame
    FunctionGradient grad = fgc(pa);
    FumiliErrorUpdator errUpdator;
    MinimumError err = errUpdator.Update(MinimumState(0), pa, fgc, 0.);
-   MinimumSeed mnseeds(MinimumState(pa, err, grad, -1., 1), st.Trafo());
+   // set an arbitrary large initial edm (1.E10)
+   MinimumSeed mnseeds(MinimumState(pa, err, grad, 1.E10, 1), st.Trafo());
 
    return ModularFunctionMinimizer::Minimize(mfcn, fgc, mnseeds, strategy, maxfcn, toler);
 }

math/minuit2/src/InitialGradientCalculator.cxx

@@ -71,7 +71,7 @@ FunctionGradient InitialGradientCalculator::operator()(const MinimumParameters &
       gr2(i) = g2;
       gst(i) = gstep;
 
-      print.Debug("Computed initial gradient for parameter", Trafo().Name(exOfIn), "value", var, "[", vmin, ",", vplu,
+      print.Trace("Computed initial gradient for parameter", Trafo().Name(exOfIn), "value", var, "[", vmin, ",", vplu,
                   "]", "dirin", dirin, "grd", grd, "g2", g2);
    }
 

math/minuit2/src/Minuit2Minimizer.cxx

@@ -487,6 +487,7 @@ bool Minuit2Minimizer::Minimize()
    fMinuitFCN->SetErrorDef(ErrorDef());
 
    const int printLevel = PrintLevel();
+   print.Debug("Minuit print level is", printLevel);
    if (PrintLevel() >= 1) {
       // print the real number of maxfcn used (defined in ModularFunctionMinimizer)
       int maxfcn_used = maxfcn;