[Minuit2] Consistently use MnFcnCaller

Author: guitargeek

math/minuit2/inc/Minuit2/FumiliBuilder.h

@@ -13,7 +13,6 @@
 #include "Minuit2/MinimumBuilder.h"
 #include "Minuit2/VariableMetricEDMEstimator.h"
 #include "Minuit2/FumiliErrorUpdator.h"
-#include "Minuit2/MnFcn.h"
 #include "Minuit2/FunctionMinimum.h"
 
 #include <vector>
@@ -22,6 +21,8 @@ namespace ROOT {
 
 namespace Minuit2 {
 
+class MnFcn;
+
 /**
 
 Builds the FunctionMinimum using the Fumili method.

math/minuit2/inc/Minuit2/MnFcn.h

@@ -40,11 +40,6 @@ class MnFcn {
    {
    }
 
-   inline double operator()(const MnAlgebraicVector &v) const
-   {
-      return fTransform ? CallWithDoingTrafo(v) : CallWithoutDoingTrafo(v);
-   }
-
    unsigned int NumOfCalls() const { return fNumCall; }
 
    double ErrorDef() const
@@ -61,19 +56,31 @@ class MnFcn {
 
    // Access the parameter transformations.
    // For internal use in the Minuit2 implementation.
-   const MnUserTransformation *transform() const { return fTransform; }
+   const MnUserTransformation *Trafo() const { return fTransform; }
 
    double CallWithTransformedParams(std::vector<double> const &vpar) const;
-
-private:
    double CallWithoutDoingTrafo(const MnAlgebraicVector &) const;
-   double CallWithDoingTrafo(const MnAlgebraicVector &) const;
 
+private:
    const FCNBase &fFCN;
    mutable int fNumCall;
    const MnUserTransformation *fTransform = nullptr;
 };
 
+// Helper class to call the MnFcn, caching the transformed parameters if necessary.
+class MnFcnCaller {
+public:
+   MnFcnCaller(const MnFcn &mfcn);
+
+   double operator()(const MnAlgebraicVector &v);
+
+private:
+   MnFcn const &fMfcn;
+   bool fDoInt2ext = false;
+   std::vector<double> fLastInput;
+   std::vector<double> fVpar;
+};
+
 } // namespace Minuit2
 
 } // namespace ROOT

math/minuit2/src/FumiliBuilder.cxx

@@ -10,7 +10,6 @@
 #include "Minuit2/FumiliBuilder.h"
 #include "Minuit2/FumiliStandardMaximumLikelihoodFCN.h"
 #include "Minuit2/GradientCalculator.h"
-//#include "Minuit2/Numerical2PGradientCalculator.h"
 #include "Minuit2/MinimumState.h"
 #include "Minuit2/MinimumError.h"
 #include "Minuit2/FunctionGradient.h"
@@ -237,6 +236,8 @@ FunctionMinimum FumiliBuilder::Minimum(const MnFcn &fcn, const GradientCalculato
 
    double lambda = (doLineSearch) ? 0.001 : 0;
 
+   MnFcnCaller fcnCaller{fcn};
+
    do {
 
       //     const MinimumState& s0 = result.back();
@@ -269,7 +270,7 @@ FunctionMinimum FumiliBuilder::Minimum(const MnFcn &fcn, const GradientCalculato
       // take a full step
 
       //evaluate function only if doing a line search
-      double fval2 = (doLineSearch) ? fcn(s0.Vec() + step) : 0;
+      double fval2 = (doLineSearch) ? fcnCaller(s0.Vec() + step) : 0;
       MinimumParameters p(s0.Vec() + step, fval2);
 
       // check that taking the full step does not deteriorate minimum
@@ -324,7 +325,7 @@ FunctionMinimum FumiliBuilder::Minimum(const MnFcn &fcn, const GradientCalculato
 
          // if the proposed point (newton step) is inside the trust region radius accept it
          if (norm <= delta) {
-            p = MinimumParameters(s0.Vec() + step, fcn(s0.Vec() + step));
+            p = MinimumParameters(s0.Vec() + step, fcnCaller(s0.Vec() + step));
             print.Debug("Accept full Newton step - it is inside TR ",delta);
          } else {
             //step = - (delta/norm) * step;
@@ -403,7 +404,7 @@ FunctionMinimum FumiliBuilder::Minimum(const MnFcn &fcn, const GradientCalculato
             }
             print.Debug("New accepted step is ",step);
 
-            p = MinimumParameters(s0.Vec() + step, fcn(s0.Vec() + step));
+            p = MinimumParameters(s0.Vec() + step, fcnCaller(s0.Vec() + step));
             norm = delta;
             gdel = inner_product(step, s0.Gradient().Grad());
          }

math/minuit2/src/FumiliMinimizer.cxx

@@ -21,7 +21,7 @@
 #include "Minuit2/MnUserParameters.h"
 #include "Minuit2/MnUserTransformation.h"
 #include "Minuit2/FumiliFCNBase.h"
-#include "Minuit2/FCNBase.h"
+#include "Minuit2/MnFcn.h"
 #include "Minuit2/MnStrategy.h"
 #include "Minuit2/MnPrint.h"
 
@@ -63,11 +63,8 @@ FunctionMinimum FumiliMinimizer::Minimize(const FCNBase &fcn, const MnUserParame
    }
 
    // compute initial values;
-   const unsigned int n = st.VariableParameters();
-   MnAlgebraicVector x(n);
-   for (unsigned int i = 0; i < n; i++)
-      x(i) = st.IntParameters()[i];
-   double fcnmin = mfcn(x);
+   MnAlgebraicVector x{st.IntParameters()};
+   double fcnmin = MnFcnCaller{mfcn}(x);
    MinimumParameters pa(x, fcnmin);
    FunctionGradient grad = fgc(pa);
    FumiliErrorUpdator errUpdator;

math/minuit2/src/HessianGradientCalculator.cxx

@@ -84,6 +84,8 @@ HessianGradientCalculator::DeltaGradient(const MinimumParameters &par, const Fun
    unsigned int startElementIndex = mpiproc.StartElementIndex();
    unsigned int endElementIndex = mpiproc.EndElementIndex();
 
+   MnFcnCaller fcnCaller{fFcn};
+
    for (unsigned int i = startElementIndex; i < endElementIndex; i++) {
       double xtf = x(i);
       double dmin = 4. * Precision().Eps2() * (xtf + Precision().Eps2());
@@ -100,9 +102,9 @@ HessianGradientCalculator::DeltaGradient(const MinimumParameters &par, const Fun
       double grdnew = 0.;
       for (unsigned int j = 0; j < Ncycle(); j++) {
          x(i) = xtf + d;
-         double fs1 = Fcn()(x);
+         double fs1 = fcnCaller(x);
          x(i) = xtf - d;
-         double fs2 = Fcn()(x);
+         double fs2 = fcnCaller(x);
          x(i) = xtf;
          //       double sag = 0.5*(fs1+fs2-2.*fcnmin);
          // LM: should I calculate also here second derivatives ???

math/minuit2/src/MnFcn.cxx

@@ -21,23 +21,6 @@ double MnFcn::CallWithoutDoingTrafo(const MnAlgebraicVector &v) const
    return fFCN(std::vector<double>{v.Data(), v.Data() + v.size()});
 }
 
-double MnFcn::CallWithDoingTrafo(const MnAlgebraicVector &v) const
-{
-   // calling fTransform() like here was not thread safe because it was using a cached vector
-   // return Fcn()( fTransform(v) );
-   // make a new thread-safe implementation creating a vector each time
-   // a bit slower few% in stressFit and 10% in Rosenbrock function but it is negligible in big fits
-
-   // get first initial values of parameter (in case some one is fixed)
-   std::vector<double> vpar(fTransform->InitialParValues().begin(), fTransform->InitialParValues().end());
-
-   for (unsigned int i = 0; i < v.size(); i++) {
-      vpar[fTransform->ExtOfInt(i)] = fTransform->Int2ext(i, v(i));
-   }
-
-   return CallWithTransformedParams(vpar);
-}
-
 // Calling the underlying function with the transformed parameters.
 // For internal use in the Minuit2 implementation.
 double MnFcn::CallWithTransformedParams(std::vector<double> const &vpar) const
@@ -48,6 +31,38 @@ double MnFcn::CallWithTransformedParams(std::vector<double> const &vpar) const
    return Fcn()(vpar);
 }
 
+MnFcnCaller::MnFcnCaller(const MnFcn &mfcn) : fMfcn{mfcn}, fDoInt2ext{static_cast<bool>(mfcn.Trafo())}
+{
+   if (!fDoInt2ext)
+      return;
+
+   MnUserTransformation const &transform = *fMfcn.Trafo();
+
+   // get first initial values of parameter (in case some one is fixed)
+   fVpar.assign(transform.InitialParValues().begin(), transform.InitialParValues().end());
+}
+
+double MnFcnCaller::operator()(const MnAlgebraicVector &v)
+{
+   if (!fDoInt2ext)
+      return fMfcn.CallWithoutDoingTrafo(v);
+
+   MnUserTransformation const &transform = *fMfcn.Trafo();
+
+   bool firstCall = fLastInput.size() != v.size();
+
+   fLastInput.resize(v.size());
+
+   for (unsigned int i = 0; i < v.size(); i++) {
+      if (firstCall || fLastInput[i] != v(i)) {
+         fVpar[transform.ExtOfInt(i)] = transform.Int2ext(i, v(i));
+         fLastInput[i] = v(i);
+      }
+   }
+
+   return fMfcn.CallWithTransformedParams(fVpar);
+}
+
 } // namespace Minuit2
 
 } // namespace ROOT

math/minuit2/src/MnFcnCaller.h

@@ -22,7 +22,7 @@
 #include "Minuit2/MnPrint.h"
 
 #include "./MPIProcess.h"
-#include "./MnFcnCaller.h"
+#include "Minuit2/MnFcn.h"
 
 #include "Math/Util.h"
 
@@ -49,11 +49,9 @@ MnHesse::operator()(const FCNBase &fcn, const MnUserParameterState &state, unsig
    // interface from MnUserParameterState
    // create a new Minimum state and use that interface
    unsigned int n = state.VariableParameters();
-   MnFcn mfcn{fcn, state.Trafo(), state.NFcn()};
-   MnAlgebraicVector x(n);
-   for (unsigned int i = 0; i < n; i++)
-      x(i) = state.IntParameters()[i];
-   double amin = mfcn(x);
+   MnFcn mfcn{fcn, state.Trafo(), static_cast<int>(state.NFcn())};
+   MnAlgebraicVector x(state.IntParameters());
+   double amin = MnFcnCaller{mfcn}(x);
    MinimumParameters par(x, amin);
    // check if we can use analytical gradient
    if (fcn.HasGradient()) {

math/minuit2/src/MnHesse.cxx

@@ -69,6 +69,8 @@ MnParabolaPoint MnLineSearch::operator()(const MnFcn &fcn, const MinimumParamete
    // start as in Fortran from 1 and count all the time we evaluate the function
    int niter = 1;
 
+   MnFcnCaller fcnCaller{fcn};
+
    for (unsigned int i = 0; i < step.size(); i++) {
       if (step(i) == 0)
          continue;
@@ -83,7 +85,7 @@ MnParabolaPoint MnLineSearch::operator()(const MnFcn &fcn, const MinimumParamete
    slamin *= prec.Eps2();
 
    double f0 = st.Fval();
-   double f1 = fcn(st.Vec() + step);
+   double f1 = fcnCaller(st.Vec() + step);
    niter++;
    double fvmin = st.Fval();
    double xvmin = 0.;
@@ -171,7 +173,7 @@ MnParabolaPoint MnLineSearch::operator()(const MnFcn &fcn, const MinimumParamete
       //     MnAlgebraicVector tmp = step;
       //     tmp *= slam;
       //     f2 = fcn(st.Vec()+tmp);
-      f2 = fcn(st.Vec() + slam * step);
+      f2 = fcnCaller(st.Vec() + slam * step);
 
       niter++; // do as in Minuit (count all func evalu)
 
@@ -263,7 +265,7 @@ MnParabolaPoint MnLineSearch::operator()(const MnFcn &fcn, const MinimumParamete
          // take the step
          //       MnAlgebraicVector tmp = step;
          //       tmp *= slam;
-         f3 = fcn(st.Vec() + slam * step);
+         f3 = fcnCaller(st.Vec() + slam * step);
          print.Trace("f3", f3, "f3-p(2-0).Y()", f3 - p2.Y(), f3 - p1.Y(), f3 - p0.Y());
          // if latest point worse than all three previous, cut step
          if (f3 > p0.Y() && f3 > p1.Y() && f3 > p2.Y()) {

math/minuit2/src/MnLineSearch.cxx

@@ -69,7 +69,7 @@ operator()(const MnFcn &fcn, const GradientCalculator &gc, const MnUserParameter
    // the line search.
    auto timingScope = std::make_unique<ROOT::Math::Util::TimingScope>([&print](std::string const &s) { print.Info(s); },
                                                                       "Evaluated function and gradient in");
-   MinimumParameters pa(x, fcn(x));
+   MinimumParameters pa(x, MnFcnCaller{fcn}(x));
    FunctionGradient dgrad = gc(pa);
    timingScope.reset();
 
@@ -155,10 +155,8 @@ MnSeedGenerator::CallWithAnalyticalGradientCalculator(const MnFcn &fcn, const An
    const MnMachinePrecision &prec = st.Precision();
 
    // initial starting values
-   MnAlgebraicVector x(n);
-   for (unsigned int i = 0; i < n; i++)
-      x(i) = st.IntParameters()[i];
-   double fcnmin = fcn(x);
+   MnAlgebraicVector x(st.IntParameters());
+   double fcnmin = MnFcnCaller{fcn}(x);
    MinimumParameters pa(x, fcnmin);
 
    // compute function gradient

math/minuit2/src/MnSeedGenerator.cxx

@@ -18,7 +18,6 @@
 #include "Minuit2/MnPrint.h"
 
 #include "./MPIProcess.h"
-#include "./MnFcnCaller.h"
 
 #ifdef _OPENMP
 #include <omp.h>

math/minuit2/src/Numerical2PGradientCalculator.cxx

@@ -36,7 +36,7 @@ std::vector<double> ParametricFunction::GetGradient(std::vector<double> const &x
 
    Numerical2PGradientCalculator gc(mfcn, st.Trafo(), strategy);
    MnAlgebraicVector xVec{x};
-   FunctionGradient g = gc(MinimumParameters{xVec, mfcn(xVec)});
+   FunctionGradient g = gc(MinimumParameters{xVec, MnFcnCaller{mfcn}(xVec)});
    const MnAlgebraicVector &grad = g.Vec();
    assert(grad.size() == x.size());
    //  std::cout << "Param Function Gradient " << grad << std::endl;