diff --git a/math/minuit2/CMakeLists.txt b/math/minuit2/CMakeLists.txt index c204e61b78ae7..8f69a7c3696c9 100644 --- a/math/minuit2/CMakeLists.txt +++ b/math/minuit2/CMakeLists.txt @@ -40,7 +40,6 @@ if(CMAKE_PROJECT_NAME STREQUAL ROOT) Minuit2/FumiliStandardMaximumLikelihoodFCN.h Minuit2/FunctionGradient.h Minuit2/FunctionMinimum.h - Minuit2/GenericFunction.h Minuit2/GradientCalculator.h Minuit2/HessianGradientCalculator.h Minuit2/InitialGradientCalculator.h @@ -150,7 +149,6 @@ if(CMAKE_PROJECT_NAME STREQUAL ROOT) src/MnStrategy.cxx src/MnTiny.cxx src/MnTraceObject.cxx - src/MnUserFcn.cxx src/MnUserParameterState.cxx src/MnUserParameters.cxx src/MnUserTransformation.cxx diff --git a/math/minuit2/inc/LinkDef.h b/math/minuit2/inc/LinkDef.h index f78ff0b5a9630..ddcab062d169b 100644 --- a/math/minuit2/inc/LinkDef.h +++ b/math/minuit2/inc/LinkDef.h @@ -13,9 +13,6 @@ #pragma link off all classes; #pragma link off all functions; -//#pragma link C++ namespace ROOT::Minuit2; - -#pragma link C++ class ROOT::Minuit2::GenericFunction; #pragma link C++ class ROOT::Minuit2::FCNBase; #pragma link C++ class ROOT::Minuit2::FCNGradientBase; #pragma link C++ class ROOT::Minuit2::FumiliFCNBase; diff --git a/math/minuit2/inc/Minuit2/FCNBase.h b/math/minuit2/inc/Minuit2/FCNBase.h index 7ab93ce5c6099..444ad7b8d7609 100644 --- a/math/minuit2/inc/Minuit2/FCNBase.h +++ b/math/minuit2/inc/Minuit2/FCNBase.h @@ -12,8 +12,6 @@ #include "Minuit2/MnConfig.h" -#include "Minuit2/GenericFunction.h" - #include #include @@ -48,10 +46,12 @@ Interface (abstract class) defining the function to be minimized, which has to b */ -class FCNBase : public GenericFunction { +class FCNBase { public: + virtual ~FCNBase() = default; + /** The meaning of the vector of parameters is of course defined by the user, @@ -74,7 +74,7 @@ class FCNBase : public GenericFunction { */ - double operator()(std::vector const &v) const override = 0; + virtual double operator()(std::vector const &v) const = 0; /** diff --git a/math/minuit2/inc/Minuit2/FumiliBuilder.h b/math/minuit2/inc/Minuit2/FumiliBuilder.h index 9b219d7ffc3b9..96b03d0e7706d 100644 --- a/math/minuit2/inc/Minuit2/FumiliBuilder.h +++ b/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 @@ -22,6 +21,8 @@ namespace ROOT { namespace Minuit2 { +class MnFcn; + /** Builds the FunctionMinimum using the Fumili method. diff --git a/math/minuit2/inc/Minuit2/GenericFunction.h b/math/minuit2/inc/Minuit2/GenericFunction.h deleted file mode 100644 index f7c65118c6262..0000000000000 --- a/math/minuit2/inc/Minuit2/GenericFunction.h +++ /dev/null @@ -1,59 +0,0 @@ -// @(#)root/minuit2:$Id$ -// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005 - -/********************************************************************** - * * - * Copyright (c) 2005 LCG ROOT Math team, CERN/PH-SFT * - * * - **********************************************************************/ - -#ifndef ROOT_Minuit2_GenericFunction -#define ROOT_Minuit2_GenericFunction - -#include "Minuit2/MnConfig.h" - -#include - -#include - -namespace ROOT { - -namespace Minuit2 { - -//_____________________________________________________________________ -/** - -Class from which all the other classes, representing functions, -inherit. That is why it defines only one method, the operator(), -which allows to call the function. - -@author Andras Zsenei and Lorenzo Moneta, Creation date: 23 Sep 2004 - -@ingroup Minuit - - */ - -class GenericFunction { - -public: - virtual ~GenericFunction() {} - - /** - - Evaluates the function using the vector containing the input values. - - @param x vector of the coordinates (for example the x coordinate for a - one-dimensional Gaussian) - - @return the result of the evaluation of the function. - - */ - - virtual double operator()(std::vector const& x) const = 0; -}; - -} // namespace Minuit2 - -} // namespace ROOT - -#endif // ROOT_Minuit2_GenericFunction diff --git a/math/minuit2/inc/Minuit2/MnFcn.h b/math/minuit2/inc/Minuit2/MnFcn.h index 2280570c4ac55..1688e303a25ac 100644 --- a/math/minuit2/inc/Minuit2/MnFcn.h +++ b/math/minuit2/inc/Minuit2/MnFcn.h @@ -10,13 +10,17 @@ #ifndef ROOT_Minuit2_MnFcn #define ROOT_Minuit2_MnFcn -#include "Minuit2/MnConfig.h" +#include "Minuit2/FCNBase.h" #include "Minuit2/MnMatrix.h" +#include + namespace ROOT { namespace Minuit2 { +class MnUserTransformation; + class FCNBase; /** Wrapper class to FCNBase interface used internally by Minuit. @@ -30,27 +34,51 @@ class FCNBase; class MnFcn { public: - /// constructor of explicit MnFcn(const FCNBase &fcn, int ncall = 0) : fFCN(fcn), fNumCall(ncall) {} + explicit MnFcn(const FCNBase &fcn, const MnUserTransformation &trafo, int ncall = 0) + : fFCN(fcn), fNumCall(ncall), fTransform(&trafo) + { + } - virtual ~MnFcn(); - - virtual double operator()(const MnAlgebraicVector &) const; unsigned int NumOfCalls() const { return fNumCall; } - // - // forward interface - // - double ErrorDef() const; - double Up() const; + double ErrorDef() const + { + return fFCN.Up(); + } + + double Up() const + { + return fFCN.Up(); + } const FCNBase &Fcn() const { return fFCN; } + // Access the parameter transformations. + // For internal use in the Minuit2 implementation. + const MnUserTransformation *Trafo() const { return fTransform; } + + double CallWithTransformedParams(std::vector const &vpar) const; + double CallWithoutDoingTrafo(const MnAlgebraicVector &) const; + private: const FCNBase &fFCN; - -protected: 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 fLastInput; + std::vector fVpar; }; } // namespace Minuit2 diff --git a/math/minuit2/inc/Minuit2/MnUserFcn.h b/math/minuit2/inc/Minuit2/MnUserFcn.h index 890d24e3e4433..9784c072bcaee 100644 --- a/math/minuit2/inc/Minuit2/MnUserFcn.h +++ b/math/minuit2/inc/Minuit2/MnUserFcn.h @@ -1,51 +1,14 @@ -// @(#)root/minuit2:$Id$ -// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005 - -/********************************************************************** - * * - * Copyright (c) 2005 LCG ROOT Math team, CERN/PH-SFT * - * * - **********************************************************************/ - #ifndef ROOT_Minuit2_MnUserFcn #define ROOT_Minuit2_MnUserFcn -#include "Minuit2/MnFcn.h" - -#include +#include namespace ROOT { - namespace Minuit2 { -class MnUserTransformation; - -/** - Wrapper used by Minuit of FCN interface - containing a reference to the transformation object - */ -class MnUserFcn : public MnFcn { - -public: - MnUserFcn(const FCNBase &fcn, const MnUserTransformation &trafo, int ncall = 0) - : MnFcn(fcn, ncall), fTransform(trafo) - { - } - - double operator()(const MnAlgebraicVector &) const override; - - // Access the parameter transformations. - // For internal use in the Minuit2 implementation. - const MnUserTransformation &transform() const { return fTransform; } - - double callWithTransformedParams(std::vector const &vpar) const; - -private: - const MnUserTransformation &fTransform; -}; - -} // namespace Minuit2 +using MnUserFcn = MnFcn; +} } // namespace ROOT -#endif // ROOT_Minuit2_MnUserFcn +#endif diff --git a/math/minuit2/inc/Minuit2/Numerical2PGradientCalculator.h b/math/minuit2/inc/Minuit2/Numerical2PGradientCalculator.h index 2a9daf062f8aa..9aea9e3e02580 100644 --- a/math/minuit2/inc/Minuit2/Numerical2PGradientCalculator.h +++ b/math/minuit2/inc/Minuit2/Numerical2PGradientCalculator.h @@ -14,17 +14,12 @@ #include "Minuit2/GradientCalculator.h" -#include - -#include - namespace ROOT { namespace Minuit2 { class MnFcn; class MnUserTransformation; -class MnMachinePrecision; class MnStrategy; /** @@ -43,15 +38,6 @@ class Numerical2PGradientCalculator : public GradientCalculator { FunctionGradient operator()(const MinimumParameters &, const FunctionGradient &) const override; - const MnFcn &Fcn() const { return fFcn; } - const MnUserTransformation &Trafo() const { return fTransformation; } - const MnMachinePrecision &Precision() const; - const MnStrategy &Strategy() const { return fStrategy; } - - unsigned int Ncycle() const; - double StepTolerance() const; - double GradTolerance() const; - private: const MnFcn &fFcn; const MnUserTransformation &fTransformation; diff --git a/math/minuit2/src/CMakeLists.txt b/math/minuit2/src/CMakeLists.txt index 630750baf233b..eca980389b832 100644 --- a/math/minuit2/src/CMakeLists.txt +++ b/math/minuit2/src/CMakeLists.txt @@ -30,7 +30,6 @@ set(MINUIT2_HEADERS FumiliStandardMaximumLikelihoodFCN.h FunctionGradient.h FunctionMinimum.h - GenericFunction.h GradientCalculator.h HessianGradientCalculator.h InitialGradientCalculator.h @@ -76,7 +75,6 @@ set(MINUIT2_HEADERS MnTiny.h MnTraceObject.h MnUserCovariance.h - MnUserFcn.h MnUserParameterState.h MnUserParameters.h MnUserTransformation.h @@ -139,7 +137,6 @@ set(MINUIT2_SOURCES MnStrategy.cxx MnTiny.cxx MnTraceObject.cxx - MnUserFcn.cxx MnUserParameterState.cxx MnUserParameters.cxx MnUserTransformation.cxx diff --git a/math/minuit2/src/CombinedMinimumBuilder.cxx b/math/minuit2/src/CombinedMinimumBuilder.cxx index 39c86c2adcf5c..dfefa459eb0b7 100644 --- a/math/minuit2/src/CombinedMinimumBuilder.cxx +++ b/math/minuit2/src/CombinedMinimumBuilder.cxx @@ -7,7 +7,6 @@ * * **********************************************************************/ -#include "Minuit2/AnalyticalGradientCalculator.h" #include "Minuit2/CombinedMinimumBuilder.h" #include "Minuit2/FunctionMinimum.h" #include "Minuit2/MnStrategy.h" @@ -37,10 +36,7 @@ FunctionMinimum CombinedMinimumBuilder::Minimum(const MnFcn &fcn, const Gradient return min1; } - // check if gradient calculator is analytical - auto agc = dynamic_cast(&gc); - MinimumSeed seed1 = (agc) ? - fVMMinimizer.SeedGenerator()(fcn, *agc, min1.UserState(), str) : fVMMinimizer.SeedGenerator()(fcn, gc, min1.UserState(), str); + MinimumSeed seed1 = fVMMinimizer.SeedGenerator()(fcn, gc, min1.UserState(), str); FunctionMinimum min2 = fVMMinimizer.Builder().Minimum(fcn, gc, seed1, str, maxfcn, edmval); if (!min2.IsValid()) { diff --git a/math/minuit2/src/FumiliBuilder.cxx b/math/minuit2/src/FumiliBuilder.cxx index f9bf7f3b93caf..077589ebd4c57 100644 --- a/math/minuit2/src/FumiliBuilder.cxx +++ b/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()); } diff --git a/math/minuit2/src/FumiliGradientCalculator.cxx b/math/minuit2/src/FumiliGradientCalculator.cxx index ff51791dad6ba..b02807bc6f3ef 100644 --- a/math/minuit2/src/FumiliGradientCalculator.cxx +++ b/math/minuit2/src/FumiliGradientCalculator.cxx @@ -17,7 +17,7 @@ #include "Minuit2/MnPrint.h" #include "Minuit2/Numerical2PGradientCalculator.h" #include "Minuit2/MnStrategy.h" -#include "Minuit2/MnUserFcn.h" +#include "Minuit2/MnFcn.h" namespace ROOT { @@ -76,7 +76,7 @@ FunctionGradient FumiliGradientCalculator::operator()(const MinimumParameters &p // 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)); + Numerical2PGradientCalculator gc(MnFcn{fFcn, fTransformation}, fTransformation, MnStrategy(1)); FunctionGradient grd2 = gc(par); os << "Fumili Gradient:" << v << std::endl; os << "Minuit Gradient" << grd2.Vec() << std::endl; diff --git a/math/minuit2/src/FumiliMinimizer.cxx b/math/minuit2/src/FumiliMinimizer.cxx index 2c0b03e8119fa..09d78d1735806 100644 --- a/math/minuit2/src/FumiliMinimizer.cxx +++ b/math/minuit2/src/FumiliMinimizer.cxx @@ -20,9 +20,8 @@ #include "Minuit2/MnUserParameterState.h" #include "Minuit2/MnUserParameters.h" #include "Minuit2/MnUserTransformation.h" -#include "Minuit2/MnUserFcn.h" #include "Minuit2/FumiliFCNBase.h" -#include "Minuit2/FCNBase.h" +#include "Minuit2/MnFcn.h" #include "Minuit2/MnStrategy.h" #include "Minuit2/MnPrint.h" @@ -40,7 +39,7 @@ FunctionMinimum FumiliMinimizer::Minimize(const FCNBase &fcn, const MnUserParame // Minimize using Fumili. Create seed and Fumili gradient calculator. // The FCNBase passed must be a FumiliFCNBase type otherwise method will fail ! - MnUserFcn mfcn(fcn, st.Trafo()); + MnFcn mfcn{fcn, st.Trafo()}; unsigned int npar = st.VariableParameters(); @@ -64,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; diff --git a/math/minuit2/src/HessianGradientCalculator.cxx b/math/minuit2/src/HessianGradientCalculator.cxx index c743369b1be50..c964842a7d1c5 100644 --- a/math/minuit2/src/HessianGradientCalculator.cxx +++ b/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 ??? diff --git a/math/minuit2/src/MnFcn.cxx b/math/minuit2/src/MnFcn.cxx index 0aaae9198116c..6defaf55a2f50 100644 --- a/math/minuit2/src/MnFcn.cxx +++ b/math/minuit2/src/MnFcn.cxx @@ -8,36 +8,59 @@ **********************************************************************/ #include "Minuit2/MnFcn.h" -#include "Minuit2/FCNBase.h" +#include "Minuit2/MnUserTransformation.h" namespace ROOT { namespace Minuit2 { -MnFcn::~MnFcn() -{ - // std::cout<<"Total number of calls to FCN: "<{v.Data(), v.Data() + v.size()}); } -// double MnFcn::operator()(const std::vector& par) const { -// return fFCN(par); -// } +// Calling the underlying function with the transformed parameters. +// For internal use in the Minuit2 implementation. +double MnFcn::CallWithTransformedParams(std::vector const &vpar) const +{ + // call Fcn function transforming from a MnAlgebraicVector of internal values to a std::vector of external ones + fNumCall++; + + return Fcn()(vpar); +} -double MnFcn::ErrorDef() const +MnFcnCaller::MnFcnCaller(const MnFcn &mfcn) : fMfcn{mfcn}, fDoInt2ext{static_cast(mfcn.Trafo())} { - return fFCN.Up(); + 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 MnFcn::Up() const +double MnFcnCaller::operator()(const MnAlgebraicVector &v) { - return fFCN.Up(); + 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 diff --git a/math/minuit2/src/MnFcnCaller.h b/math/minuit2/src/MnFcnCaller.h deleted file mode 100644 index d7df64392eeee..0000000000000 --- a/math/minuit2/src/MnFcnCaller.h +++ /dev/null @@ -1,56 +0,0 @@ -#ifndef ROOT_Minuit2_MnFcnCaller -#define ROOT_Minuit2_MnFcnCaller - -#include "Minuit2/MnUserFcn.h" - -namespace ROOT { - -namespace Minuit2 { - -// Helper class to call the MnFcn, cashing the transformed parameters in case -// it is a MnUserFcn that does the parameter transformation when calling. -class MnFcnCaller { -public: - MnFcnCaller(const MnFcn &mfcn) : fMfcn{mfcn}, fIsUserFcn{static_cast(dynamic_cast(&mfcn))} - { - if (!fIsUserFcn) - return; - - MnUserTransformation const &transform = static_cast(fMfcn).transform(); - - // get first initial values of parameter (in case some one is fixed) - fVpar.assign(transform.InitialParValues().begin(), transform.InitialParValues().end()); - } - - double operator()(const MnAlgebraicVector &v) - { - if (!fIsUserFcn) - return fMfcn(v); - - MnUserTransformation const &transform = static_cast(fMfcn).transform(); - - 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 static_cast(fMfcn).callWithTransformedParams(fVpar); - } - -private: - MnFcn const &fMfcn; - bool fIsUserFcn = false; - std::vector fLastInput; - std::vector fVpar; -}; - -} // namespace Minuit2 -} // namespace ROOT - -#endif diff --git a/math/minuit2/src/MnHesse.cxx b/math/minuit2/src/MnHesse.cxx index 2fe8c01d5096a..eb3b6bb52b46d 100644 --- a/math/minuit2/src/MnHesse.cxx +++ b/math/minuit2/src/MnHesse.cxx @@ -9,7 +9,6 @@ #include "Minuit2/MnHesse.h" #include "Minuit2/MnUserParameterState.h" -#include "Minuit2/MnUserFcn.h" #include "Minuit2/FCNBase.h" #include "Minuit2/MnPosDef.h" #include "Minuit2/HessianGradientCalculator.h" @@ -23,7 +22,7 @@ #include "Minuit2/MnPrint.h" #include "./MPIProcess.h" -#include "./MnFcnCaller.h" +#include "Minuit2/MnFcn.h" #include "Math/Util.h" @@ -50,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(); - MnUserFcn 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(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()) { @@ -79,7 +76,7 @@ void MnHesse::operator()(const FCNBase &fcn, FunctionMinimum &min, unsigned int // interface from FunctionMinimum to be used after minimization // use last state from the minimization without the need to re-create a new state // do not reset function calls and keep updating them - MnUserFcn mfcn(fcn, min.UserState().Trafo(), min.NFcn()); + MnFcn mfcn{fcn, min.UserState().Trafo(), min.NFcn()}; MinimumState st = (*this)(mfcn, min.State(), min.UserState().Trafo(), maxcalls); min.Add(st); } diff --git a/math/minuit2/src/MnLineSearch.cxx b/math/minuit2/src/MnLineSearch.cxx index 2b981255e6817..961e406e9a22b 100644 --- a/math/minuit2/src/MnLineSearch.cxx +++ b/math/minuit2/src/MnLineSearch.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()) { diff --git a/math/minuit2/src/MnSeedGenerator.cxx b/math/minuit2/src/MnSeedGenerator.cxx index e668c6c3d09a3..8161ca5eb0e4f 100644 --- a/math/minuit2/src/MnSeedGenerator.cxx +++ b/math/minuit2/src/MnSeedGenerator.cxx @@ -69,7 +69,7 @@ operator()(const MnFcn &fcn, const GradientCalculator &gc, const MnUserParameter // the line search. auto timingScope = std::make_unique([&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 diff --git a/math/minuit2/src/MnUserFcn.cxx b/math/minuit2/src/MnUserFcn.cxx deleted file mode 100644 index f9a1c9aa6a2b1..0000000000000 --- a/math/minuit2/src/MnUserFcn.cxx +++ /dev/null @@ -1,54 +0,0 @@ -// @(#)root/minuit2:$Id$ -// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005 - -/********************************************************************** - * * - * Copyright (c) 2005 LCG ROOT Math team, CERN/PH-SFT * - * * - **********************************************************************/ - -#include "Minuit2/MnUserFcn.h" -#include "Minuit2/FCNBase.h" -#include "Minuit2/MnUserTransformation.h" - -namespace ROOT { - -namespace Minuit2 { - -double MnUserFcn::operator()(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 vpar(fTransform.InitialParValues().begin(), fTransform.InitialParValues().end()); - - const std::vector ¶meters = fTransform.Parameters(); - unsigned int n = v.size(); - for (unsigned int i = 0; i < n; i++) { - int ext = fTransform.ExtOfInt(i); - if (parameters[ext].HasLimits()) { - vpar[ext] = fTransform.Int2ext(i, v(i)); - } else { - vpar[ext] = v(i); - } - } - - return callWithTransformedParams(vpar); -} - -// Calling the underlying function with the transformed parameters. -// For internal use in the Minuit2 implementation. -double MnUserFcn::callWithTransformedParams(std::vector const &vpar) const -{ - // call Fcn function transforming from a MnAlgebraicVector of internal values to a std::vector of external ones - fNumCall++; - - return Fcn()(vpar); -} - -} // namespace Minuit2 - -} // namespace ROOT diff --git a/math/minuit2/src/ModularFunctionMinimizer.cxx b/math/minuit2/src/ModularFunctionMinimizer.cxx index 4fa546c360760..cd4d0da79cae7 100644 --- a/math/minuit2/src/ModularFunctionMinimizer.cxx +++ b/math/minuit2/src/ModularFunctionMinimizer.cxx @@ -19,12 +19,12 @@ #include "Minuit2/MnUserParameters.h" #include "Minuit2/MnUserCovariance.h" #include "Minuit2/MnUserTransformation.h" -#include "Minuit2/MnUserFcn.h" #include "Minuit2/FCNBase.h" #include "Minuit2/MnStrategy.h" #include "Minuit2/MnLineSearch.h" #include "Minuit2/MnParabolaPoint.h" #include "Minuit2/MnPrint.h" +#include "Minuit2/MnFcn.h" namespace ROOT { @@ -33,20 +33,18 @@ namespace Minuit2 { FunctionMinimum ModularFunctionMinimizer::Minimize(const FCNBase &fcn, const MnUserParameterState &st, const MnStrategy &strategy, unsigned int maxfcn, double toler) const { - // need MnUserFcn for difference int-ext parameters - MnUserFcn mfcn(fcn, st.Trafo()); + // need wrapper for difference int-ext parameters + MnFcn mfcn{fcn, st.Trafo()}; std::unique_ptr gc; if (!fcn.HasGradient()) { // minimize from a FCNBase and a MnUserparameterState - interface used by all the previous ones // based on FCNBase. Create in this case a NumericalGradient calculator - // Create the minuit FCN wrapper (MnUserFcn) containing the transformation (int<->ext) gc = std::make_unique(mfcn, st.Trafo(), strategy); } else if (fcn.gradParameterSpace() == GradientParameterSpace::Internal) { // minimize from a function with gradient and a MnUserParameterState - // interface based on function with gradient (external/analytical gradients) // Create in this case an AnalyticalGradient calculator - // Create the minuit FCN wrapper (MnUserFcn) containing the transformation (int<->ext) gc = std::make_unique(fcn, st.Trafo()); } else { gc = std::make_unique(fcn, st.Trafo()); diff --git a/math/minuit2/src/Numerical2PGradientCalculator.cxx b/math/minuit2/src/Numerical2PGradientCalculator.cxx index a511ccc09c7b6..0530b2dc795f2 100644 --- a/math/minuit2/src/Numerical2PGradientCalculator.cxx +++ b/math/minuit2/src/Numerical2PGradientCalculator.cxx @@ -18,7 +18,6 @@ #include "Minuit2/MnPrint.h" #include "./MPIProcess.h" -#include "./MnFcnCaller.h" #ifdef _OPENMP #include @@ -48,20 +47,20 @@ operator()(const MinimumParameters &par, const FunctionGradient &Gradient) const // calculate numerical gradient from MinimumParameters object // the algorithm takes correctly care when the gradient is approximately zero + MnUserTransformation const &trafo = fTransformation; + // std::cout<<"########### Numerical2PDerivative"< 0.5) step = 0.5; } @@ -130,7 +131,7 @@ operator()(const MinimumParameters &par, const FunctionGradient &Gradient) const step = stpmin; // std::cout<<" "< ParametricFunction::GetGradient(std::vector 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; diff --git a/math/minuit2/src/SimplexBuilder.cxx b/math/minuit2/src/SimplexBuilder.cxx index 490f6b59d0d7e..1e5163adc782e 100644 --- a/math/minuit2/src/SimplexBuilder.cxx +++ b/math/minuit2/src/SimplexBuilder.cxx @@ -28,6 +28,8 @@ FunctionMinimum SimplexBuilder::Minimum(const MnFcn &mfcn, const GradientCalcula // method to find initial simplex is slightly different than in the original Fortran // Minuit since has not been proofed that one to be better + MnFcnCaller mfcnCaller{mfcn}; + // synchronize print levels MnPrint print("SimplexBuilder", PrintLevel()); @@ -41,7 +43,7 @@ FunctionMinimum SimplexBuilder::Minimum(const MnFcn &mfcn, const GradientCalcula // If there are no parameters, we just return a minimum state if (n == 0) { - MinimumState st(MinimumParameters(MnAlgebraicVector{0}, MnAlgebraicVector{0}, mfcn(x)), 0.0, + MinimumState st(MinimumParameters(MnAlgebraicVector{0}, MnAlgebraicVector{0}, mfcnCaller(x)), 0.0, mfcn.NumOfCalls()); return FunctionMinimum(seed, {&st, 1}, mfcn.Up()); } @@ -71,7 +73,7 @@ FunctionMinimum SimplexBuilder::Minimum(const MnFcn &mfcn, const GradientCalcula if (step(i) < dmin) step(i) = dmin; x(i) += step(i); - double tmp = mfcn(x); + double tmp = mfcnCaller(x); if (tmp < amin) { amin = tmp; jl = i + 1; @@ -123,7 +125,7 @@ FunctionMinimum SimplexBuilder::Minimum(const MnFcn &mfcn, const GradientCalcula } MnAlgebraicVector pstar = (1. + alpha) * pbar - alpha * simplex(jh).second; - const double ystar = mfcn(pstar); + const double ystar = mfcnCaller(pstar); print.Debug("pbar", pbar, "pstar", pstar, "f(pstar)", ystar); @@ -134,7 +136,7 @@ FunctionMinimum SimplexBuilder::Minimum(const MnFcn &mfcn, const GradientCalcula continue; } MnAlgebraicVector pstst = beta * simplex(jh).second + (1. - beta) * pbar; - double ystst = mfcn(pstst); + double ystst = mfcnCaller(pstst); print.Debug("Reduced simplex pstst", pstst, "f(pstst)", ystst); @@ -145,7 +147,7 @@ FunctionMinimum SimplexBuilder::Minimum(const MnFcn &mfcn, const GradientCalcula } MnAlgebraicVector pstst = gamma * pstar + (1. - gamma) * pbar; - double ystst = mfcn(pstst); + double ystst = mfcnCaller(pstst); print.Debug("pstst", pstst, "f(pstst)", ystst); @@ -162,7 +164,7 @@ FunctionMinimum SimplexBuilder::Minimum(const MnFcn &mfcn, const GradientCalcula if (rho > rhomax) rho = rhomax; MnAlgebraicVector prho = rho * pbar + (1. - rho) * simplex(jh).second; - double yrho = mfcn(prho); + double yrho = mfcnCaller(prho); print.Debug("prho", prho, "f(prho)", yrho); @@ -183,7 +185,7 @@ FunctionMinimum SimplexBuilder::Minimum(const MnFcn &mfcn, const GradientCalcula } if (ystar > simplex(jh).first) { pstst = beta * simplex(jh).second + (1. - beta) * pbar; - ystst = mfcn(pstst); + ystst = mfcnCaller(pstst); if (ystst > simplex(jh).first) break; simplex.Update(ystst, pstst); @@ -203,7 +205,7 @@ FunctionMinimum SimplexBuilder::Minimum(const MnFcn &mfcn, const GradientCalcula continue; pbar += (wg * simplex(i).second); } - double ybar = mfcn(pbar); + double ybar = mfcnCaller(pbar); if (ybar < amin) simplex.Update(ybar, pbar); else { diff --git a/math/minuit2/src/SimplexSeedGenerator.cxx b/math/minuit2/src/SimplexSeedGenerator.cxx index c477b3fc209f4..3a740daf72763 100644 --- a/math/minuit2/src/SimplexSeedGenerator.cxx +++ b/math/minuit2/src/SimplexSeedGenerator.cxx @@ -28,10 +28,8 @@ operator()(const MnFcn &fcn, const GradientCalculator &, const MnUserParameterSt 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); InitialGradientCalculator igc(fcn, st.Trafo()); FunctionGradient dgrad = igc(pa); diff --git a/math/minuit2/test/MnSim/GaussianModelFunction.h b/math/minuit2/test/MnSim/GaussianModelFunction.h index 88f924a28cfec..9adc1eb24cd50 100644 --- a/math/minuit2/test/MnSim/GaussianModelFunction.h +++ b/math/minuit2/test/MnSim/GaussianModelFunction.h @@ -14,7 +14,6 @@ #include "Minuit2/ParametricFunction.h" -#include "Minuit2/MnFcn.h" #include "Minuit2/MnStrategy.h" #include "Minuit2/MnUserParameterState.h"