diff --git a/math/minuit2/inc/Minuit2/MnStrategy.h b/math/minuit2/inc/Minuit2/MnStrategy.h index a1cb48df85844..8ff00fb609af0 100644 --- a/math/minuit2/inc/Minuit2/MnStrategy.h +++ b/math/minuit2/inc/Minuit2/MnStrategy.h @@ -33,8 +33,6 @@ class MnStrategy { // user defined strategy (0, 1, 2, >=3) explicit MnStrategy(unsigned int); - unsigned int Strategy() const { return fStrategy; } - unsigned int GradientNCycles() const { return fGradNCyc; } double GradientStepTolerance() const { return fGradTlrStp; } double GradientTolerance() const { return fGradTlr; } @@ -48,15 +46,11 @@ class MnStrategy { int StorageLevel() const { return fStoreLevel; } - bool IsLow() const { return fStrategy == 0; } - bool IsMedium() const { return fStrategy == 1; } - bool IsHigh() const { return fStrategy == 2; } - bool IsVeryHigh() const { return fStrategy >= 3; } + bool RefineGradientInHessian() const { return fStrategy > 0; } - void SetLowStrategy(); - void SetMediumStrategy(); - void SetHighStrategy(); - void SetVeryHighStrategy(); + bool ComputeInitialHessian() const { return fStrategy == 2; } + + double HessianRecomputeThreshold() const; void SetGradientNCycles(unsigned int n) { fGradNCyc = n; } void SetGradientStepTolerance(double stp) { fGradTlrStp = stp; } @@ -80,6 +74,15 @@ class MnStrategy { void SetStorageLevel(unsigned int level) { fStoreLevel = level; } private: + friend class MnFunctionCross; + friend class MnContours; + MnStrategy NextLower() const; + + void SetLowStrategy(); + void SetMediumStrategy(); + void SetHighStrategy(); + void SetVeryHighStrategy(); + unsigned int fStrategy; unsigned int fGradNCyc; diff --git a/math/minuit2/src/FumiliBuilder.cxx b/math/minuit2/src/FumiliBuilder.cxx index 34182056ad6bf..196e8f6c1106e 100644 --- a/math/minuit2/src/FumiliBuilder.cxx +++ b/math/minuit2/src/FumiliBuilder.cxx @@ -97,7 +97,7 @@ FunctionMinimum FumiliBuilder::Minimum(const MnFcn &fcn, const GradientCalculato // call always Hesse (error matrix from Fumili is never accurate since is approximate) - if (strategy.Strategy() > 0) { + if (min.Error().Dcovar() > strategy.HessianRecomputeThreshold()) { print.Debug("FumiliBuilder will verify convergence and Error matrix; " "dcov", min.Error().Dcovar()); diff --git a/math/minuit2/src/MnContours.cxx b/math/minuit2/src/MnContours.cxx index 04b2e5f7af6d5..2fde939a3d736 100644 --- a/math/minuit2/src/MnContours.cxx +++ b/math/minuit2/src/MnContours.cxx @@ -84,7 +84,7 @@ ContoursError MnContours::Contour(unsigned int px, unsigned int py, unsigned int print.Debug("Minos error for p0: ",ey.first,ey.second); // 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)))); + MnMigrad migrad0(fFCN, fMinimum.UserState(), fStrategy.NextLower()); // function cross for a single parameter - copy the state MnUserParameterState ustate = fMinimum.UserState(); @@ -190,7 +190,7 @@ ContoursError MnContours::Contour(unsigned int px, unsigned int py, unsigned int // now look for x values when y is fixed - MnMigrad migrad1(fFCN, fMinimum.UserState(), MnStrategy(std::max(0, int(fStrategy.Strategy() - 1)))); + MnMigrad migrad1(fFCN, fMinimum.UserState(), fStrategy.NextLower()); migrad1.State().Fix(py); std::vector xvalues_ylo(1); migrad1.State().SetValue(py, valy + ey.first); diff --git a/math/minuit2/src/MnFunctionCross.cxx b/math/minuit2/src/MnFunctionCross.cxx index 2c2a17ef17fcc..e20738a301e70 100644 --- a/math/minuit2/src/MnFunctionCross.cxx +++ b/math/minuit2/src/MnFunctionCross.cxx @@ -104,7 +104,7 @@ MnCross MnFunctionCross::operator()(std::span par, std::span if (aulim < aopt + tla) limset = true; - MnMigrad migrad(fFCN, fState, MnStrategy(std::max(0, int(fStrategy.Strategy() - 1)))); + MnMigrad migrad(fFCN, fState, fStrategy.NextLower()); print.Info([&](std::ostream &os) { os << "Run Migrad with fixed parameters:"; diff --git a/math/minuit2/src/MnHesse.cxx b/math/minuit2/src/MnHesse.cxx index eb3b6bb52b46d..a5fa41edc3c61 100644 --- a/math/minuit2/src/MnHesse.cxx +++ b/math/minuit2/src/MnHesse.cxx @@ -310,7 +310,7 @@ MinimumState ComputeNumerical(const MnFcn &mfcn, const MinimumState &st, const M print.Debug("Second derivatives", g2); - if (strat.Strategy() > 0) { + if (strat.RefineGradientInHessian()) { // refine first derivative HessianGradientCalculator hgc(mfcn, trafo, strat); FunctionGradient gr = hgc(st.Parameters(), FunctionGradient(grd, g2, gst)); diff --git a/math/minuit2/src/MnSeedGenerator.cxx b/math/minuit2/src/MnSeedGenerator.cxx index 6e7a0d6129c77..a40506fb049d2 100644 --- a/math/minuit2/src/MnSeedGenerator.cxx +++ b/math/minuit2/src/MnSeedGenerator.cxx @@ -106,7 +106,7 @@ operator()(const MnFcn &fcn, const GradientCalculator &gc, const MnUserParameter } } - if (stra.Strategy() == 2 && !st.HasCovariance()) { + if (stra.ComputeInitialHessian() && !st.HasCovariance()) { // calculate full 2nd derivative print.Debug("calling MnHesse"); @@ -234,7 +234,7 @@ MnSeedGenerator::CallWithAnalyticalGradientCalculator(const MnFcn &fcn, const An } // compute Hessian above will not have posdef check as it is done if we call MnHesse - if (stra.Strategy() == 2 && !st.HasCovariance() && !computedHessian) { + if (stra.ComputeInitialHessian() && !st.HasCovariance() && !computedHessian) { // can calculate full 2nd derivative MinimumState tmpState = MnHesse(stra)(fcn, state, st.Trafo()); print.Info("Compute full Hessian: Initial seeding state is ",tmpState); diff --git a/math/minuit2/src/MnStrategy.cxx b/math/minuit2/src/MnStrategy.cxx index 3cd36a939c03a..e44727380d34c 100644 --- a/math/minuit2/src/MnStrategy.cxx +++ b/math/minuit2/src/MnStrategy.cxx @@ -9,6 +9,9 @@ #include "Minuit2/MnStrategy.h" +#include +#include + namespace ROOT { namespace Minuit2 { @@ -89,6 +92,20 @@ void MnStrategy::SetVeryHighStrategy() SetHessianForcePosDef(0); } +MnStrategy MnStrategy::NextLower() const +{ + return MnStrategy(std::max(0, int(fStrategy - 1))); +} + +double MnStrategy::HessianRecomputeThreshold() const +{ + if (fStrategy == 0) + return std::numeric_limits::infinity(); + if (fStrategy == 1) + return 0.05; + return -std::numeric_limits::infinity(); +} + } // namespace Minuit2 } // namespace ROOT diff --git a/math/minuit2/src/VariableMetricBuilder.cxx b/math/minuit2/src/VariableMetricBuilder.cxx index d2b8fcc506c63..84bdac9fede22 100644 --- a/math/minuit2/src/VariableMetricBuilder.cxx +++ b/math/minuit2/src/VariableMetricBuilder.cxx @@ -130,7 +130,7 @@ FunctionMinimum VariableMetricBuilder::Minimum(const MnFcn &fcn, const GradientC edm = result.back().Edm(); // need to correct again for Dcovar: edm *= (1. + 3. * e.Dcovar()) ??? - if ((strategy.Strategy() >= 2) || (strategy.Strategy() == 1 && min.Error().Dcovar() > 0.05)) { + if (min.Error().Dcovar() > strategy.HessianRecomputeThreshold()) { print.Debug("MnMigrad will verify convergence and Error matrix; dcov =", min.Error().Dcovar());