From a9a4c16391bad2dc614b4509cc5d3352d3baf0ef Mon Sep 17 00:00:00 2001 From: moneta Date: Wed, 29 Jan 2025 11:45:32 +0100 Subject: [PATCH] [math][minuit2] Apply small fixes and improvements on new Fumili version Fix also fit option V, VV and VVV --- hist/hist/src/HFitImpl.cxx | 10 ++++++---- math/minuit2/inc/Minuit2/FumiliFCNAdapter.h | 10 +++------- math/minuit2/inc/Minuit2/FumiliFCNBase.h | 1 - math/minuit2/src/FumiliBuilder.cxx | 9 ++++++--- math/minuit2/src/FumiliErrorUpdator.cxx | 7 ++++++- math/minuit2/src/FumiliGradientCalculator.cxx | 11 +++++++---- math/minuit2/src/FumiliMinimizer.cxx | 11 ++++------- math/minuit2/src/InitialGradientCalculator.cxx | 2 +- math/minuit2/src/Minuit2Minimizer.cxx | 1 + 9 files changed, 34 insertions(+), 28 deletions(-) diff --git a/hist/hist/src/HFitImpl.cxx b/hist/hist/src/HFitImpl.cxx index 99ccf9b383574..f1a46467b8aed 100644 --- a/hist/hist/src/HFitImpl.cxx +++ b/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; diff --git a/math/minuit2/inc/Minuit2/FumiliFCNAdapter.h b/math/minuit2/inc/Minuit2/FumiliFCNAdapter.h index e5c6995fc06f8..8060c08f938c4 100644 --- a/math/minuit2/inc/Minuit2/FumiliFCNAdapter.h +++ b/math/minuit2/inc/Minuit2/FumiliFCNAdapter.h @@ -54,13 +54,9 @@ class FumiliFCNAdapter : public FumiliFCNBase { void SetErrorDef(double up) override { fUp = up; } - // virtual std::vector Gradient(std::vector 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 const &v) override; @@ -130,7 +126,7 @@ void FumiliFCNAdapter::EvaluateAll(std::vector 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::EvaluateAll(std::vector 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"); } diff --git a/math/minuit2/inc/Minuit2/FumiliFCNBase.h b/math/minuit2/inc/Minuit2/FumiliFCNBase.h index 1e754b9c03b91..5fa2b15a04a93 100644 --- a/math/minuit2/inc/Minuit2/FumiliFCNBase.h +++ b/math/minuit2/inc/Minuit2/FumiliFCNBase.h @@ -70,7 +70,6 @@ class FumiliFCNBase : public FCNBase { { } - // FumiliFCNBase(const ParametricFunction& modelFCN) { fModelFunction = &modelFCN; } ~FumiliFCNBase() override {} diff --git a/math/minuit2/src/FumiliBuilder.cxx b/math/minuit2/src/FumiliBuilder.cxx index 9200e8f0f8d2f..7371c7f6f19f9 100644 --- a/math/minuit2/src/FumiliBuilder.cxx +++ b/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; } diff --git a/math/minuit2/src/FumiliErrorUpdator.cxx b/math/minuit2/src/FumiliErrorUpdator.cxx index 6c2ec7dcf887f..800a0a593590f 100644 --- a/math/minuit2/src/FumiliErrorUpdator.cxx +++ b/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(const_cast(&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; diff --git a/math/minuit2/src/FumiliGradientCalculator.cxx b/math/minuit2/src/FumiliGradientCalculator.cxx index 09c65f3b6acbf..ff51791dad6ba 100644 --- a/math/minuit2/src/FumiliGradientCalculator.cxx +++ b/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 extParam = fTransformation(par.Vec()); // eval Gradient FumiliFCNBase &fcn = const_cast(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); }); diff --git a/math/minuit2/src/FumiliMinimizer.cxx b/math/minuit2/src/FumiliMinimizer.cxx index 8fb2ddc03bb92..a265615397135 100644 --- a/math/minuit2/src/FumiliMinimizer.cxx +++ b/math/minuit2/src/FumiliMinimizer.cxx @@ -41,7 +41,7 @@ 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) @@ -49,11 +49,7 @@ FunctionMinimum FumiliMinimizer::Minimize(const FCNBase &fcn, const MnUserParame // 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(const_cast(&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); } diff --git a/math/minuit2/src/InitialGradientCalculator.cxx b/math/minuit2/src/InitialGradientCalculator.cxx index 91823526022ea..fb179f2f41dae 100644 --- a/math/minuit2/src/InitialGradientCalculator.cxx +++ b/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); } diff --git a/math/minuit2/src/Minuit2Minimizer.cxx b/math/minuit2/src/Minuit2Minimizer.cxx index 1f148a72cd061..37c33804aa74e 100644 --- a/math/minuit2/src/Minuit2Minimizer.cxx +++ b/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;