[Minuit2] Typesafe strategy handling

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;

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());

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<double> xvalues_ylo(1);
    migrad1.State().SetValue(py, valy + ey.first);

math/minuit2/src/MnFunctionCross.cxx

@@ -104,7 +104,7 @@ MnCross MnFunctionCross::operator()(std::span<const unsigned int> 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:";

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));

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);

math/minuit2/src/MnStrategy.cxx

@@ -9,6 +9,9 @@
 
 #include "Minuit2/MnStrategy.h"
 
+#include <algorithm>
+#include <limits>
+
 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<double>::infinity();
+   if (fStrategy == 1)
+      return 0.05;
+   return -std::numeric_limits<double>::infinity();
+}
+
 } // namespace Minuit2
 
 } // namespace ROOT

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());