From caac0876d8a8191a2162027e0e012edec79f7f8e Mon Sep 17 00:00:00 2001 From: "E. G. Patrick Bos" Date: Fri, 10 Nov 2017 13:13:04 +0100 Subject: [PATCH] Working on extending Minuit2 Analytical gradient implementation --- math/minuit2/inc/Minuit2/FCNGradAdapter.h | 103 ++- math/minuit2/inc/Minuit2/FCNGradientBase.h | 28 +- .../inc/Minuit2/MnUserTransformation.h | 239 +++--- .../inc/Minuit2/SinParameterTransformation.h | 30 +- .../Minuit2/SqrtLowParameterTransformation.h | 34 +- .../Minuit2/SqrtUpParameterTransformation.h | 36 +- .../src/AnalyticalGradientCalculator.cxx | 13 +- .../minuit2/src/InitialGradientCalculator.cxx | 2 + math/minuit2/src/MnSeedGenerator.cxx | 56 +- math/minuit2/src/MnUserTransformation.cxx | 758 +++++++++--------- math/minuit2/src/NegativeG2LineSearch.cxx | 14 +- .../src/Numerical2PGradientCalculator.cxx | 12 +- .../src/SinParameterTransformation.cxx | 72 +- .../src/SqrtLowParameterTransformation.cxx | 55 +- .../src/SqrtUpParameterTransformation.cxx | 56 +- math/minuit2/src/VariableMetricBuilder.cxx | 1 + .../inc/NumericalDerivatorMinuit2.h | 135 ++-- roofit/roofitcore/inc/RooGradMinimizerFcn.h | 103 ++- .../src/NumericalDerivatorMinuit2.cxx | 576 +++++++------ roofit/roofitcore/src/RooGradMinimizerFcn.cxx | 60 +- 20 files changed, 1340 insertions(+), 1043 deletions(-) diff --git a/math/minuit2/inc/Minuit2/FCNGradAdapter.h b/math/minuit2/inc/Minuit2/FCNGradAdapter.h index cc121d7496eea..2ed0f0df77290 100644 --- a/math/minuit2/inc/Minuit2/FCNGradAdapter.h +++ b/math/minuit2/inc/Minuit2/FCNGradAdapter.h @@ -1,9 +1,10 @@ // @(#)root/minuit2:$Id$ -// Author: L. Moneta 10/2006 +// Authors: L. Moneta, E.G.P. Bos 2006-2017 /********************************************************************** * * * Copyright (c) 2006 ROOT Foundation, CERN/PH-SFT * + * Copyright (c) 2017 Patrick Bos, Netherlands eScience Center * * * **********************************************************************/ @@ -19,65 +20,85 @@ namespace ROOT { - namespace Minuit2 { + namespace Minuit2 { -/** + /** + template wrapped class for adapting to FCNBase signature a IGradFunction -template wrapped class for adapting to FCNBase signature a IGradFunction + @author Lorenzo Moneta, Patrick Bos -@author Lorenzo Moneta + @ingroup Minuit -@ingroup Minuit + */ -*/ + template< class Function> + class FCNGradAdapter : public FCNGradientBase { -template< class Function> -class FCNGradAdapter : public FCNGradientBase { + public: -public: + FCNGradAdapter(const Function & f, double up = 1.) : + fFunc(f) , + fUp (up) , + fGrad(std::vector(fFunc.NDim() ) ), + fG2(fFunc.hasG2ndDerivative() ? std::vector(fFunc.NDim()) : std::vector(0) ), + fGStep(fFunc.hasGStepSize() ? std::vector(fFunc.NDim()) : std::vector(0) ) + {} - FCNGradAdapter(const Function & f, double up = 1.) : - fFunc(f) , - fUp (up) , - fGrad(std::vector(fFunc.NDim() ) ) + ~FCNGradAdapter() {} - {} - ~FCNGradAdapter() {} + double operator()(const std::vector& v) const { + return fFunc.operator()(&v[0]); + } + double operator()(const double * v) const { + return fFunc.operator()(v); + } + double Up() const {return fUp;} - double operator()(const std::vector& v) const { - return fFunc.operator()(&v[0]); - } - double operator()(const double * v) const { - return fFunc.operator()(v); - } - - double Up() const {return fUp;} - - std::vector Gradient(const std::vector& v) const { - fFunc.Gradient(&v[0], &fGrad[0]); + virtual std::vector Gradient(const std::vector& v) const { + fFunc.Gradient(v.data(), fGrad.data()); #ifdef DEBUG - std::cout << " gradient in FCNAdapter = { " ; + std::cout << " gradient in FCNAdapter = { " ; for (unsigned int i = 0; i < fGrad.size(); ++i) std::cout << fGrad[i] << "\t"; std::cout << "}" << std::endl; #endif - return fGrad; - } - // forward interface - //virtual double operator()(int npar, double* params,int iflag = 4) const; - bool CheckGradient() const { return false; } - -private: - const Function & fFunc; - double fUp; - mutable std::vector fGrad; -}; - - } // end namespace Minuit2 + return fGrad; + } + // forward interface + //virtual double operator()(int npar, double* params,int iflag = 4) const; + bool CheckGradient() const { return false; } + + virtual std::vector G2ndDerivative(const std::vector& v) const { + fFunc.G2ndDerivative(v.data(), fG2.data()); + return fG2; + }; + + virtual std::vector GStepSize(const std::vector& v) const { + fFunc.GStepSize(v.data(), fGStep.data()); + return fGStep; + }; + + virtual bool hasG2ndDerivative() const { + return fFunc.hasG2ndDerivative(); + } + + virtual bool hasGStepSize() const { + return fFunc.hasGStepSize(); + } + + private: + const Function & fFunc; + double fUp; + mutable std::vector fGrad; + mutable std::vector fG2; + mutable std::vector fGStep; + }; + + } // end namespace Minuit2 } // end namespace ROOT diff --git a/math/minuit2/inc/Minuit2/FCNGradientBase.h b/math/minuit2/inc/Minuit2/FCNGradientBase.h index ba9a7fe6232ea..aaf660e771ba8 100644 --- a/math/minuit2/inc/Minuit2/FCNGradientBase.h +++ b/math/minuit2/inc/Minuit2/FCNGradientBase.h @@ -1,9 +1,10 @@ // @(#)root/minuit2:$Id$ -// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005 +// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei, E.G.P. Bos 2003-2017 /********************************************************************** * * * Copyright (c) 2005 LCG ROOT Math team, CERN/PH-SFT * + * Copyright (c) 2017 Patrick Bos, Netherlands eScience Center * * * **********************************************************************/ @@ -14,7 +15,7 @@ namespace ROOT { - namespace Minuit2 { + namespace Minuit2 { //________________________________________________________________________ /** Extension of the FCNBase for providing the analytical Gradient of the @@ -29,17 +30,28 @@ namespace ROOT { "false". */ -class FCNGradientBase : public FCNBase { + class FCNGradientBase : public FCNBase { -public: + public: - virtual ~FCNGradientBase() {} + virtual ~FCNGradientBase() {} - virtual std::vector Gradient(const std::vector&) const = 0; + virtual std::vector Gradient(const std::vector&) const = 0; - virtual bool CheckGradient() const {return true;} + virtual std::vector G2ndDerivative(const std::vector&) const = 0; + virtual std::vector GStepSize(const std::vector&) const = 0; -}; + virtual bool hasG2ndDerivative() const { + return false; + } + + virtual bool hasGStepSize() const { + return false; + } + + virtual bool CheckGradient() const {return true;} + + }; } // namespace Minuit2 diff --git a/math/minuit2/inc/Minuit2/MnUserTransformation.h b/math/minuit2/inc/Minuit2/MnUserTransformation.h index b6d349bcb4d55..3c3931f27e0d0 100644 --- a/math/minuit2/inc/Minuit2/MnUserTransformation.h +++ b/math/minuit2/inc/Minuit2/MnUserTransformation.h @@ -1,9 +1,10 @@ // @(#)root/minuit2:$Id$ -// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005 +// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei, E.G.P. Bos 2003-2017 /********************************************************************** * * * Copyright (c) 2005 LCG ROOT Math team, CERN/PH-SFT * + * Copyright (c) 2017 Patrick Bos, Netherlands eScience Center * * * **********************************************************************/ @@ -23,10 +24,10 @@ namespace ROOT { - namespace Minuit2 { + namespace Minuit2 { -class MnUserCovariance; + class MnUserCovariance; // class MnMachinePrecision; @@ -35,162 +36,164 @@ class MnUserCovariance; internal parameters used for minimization */ -class MnUserTransformation { - -public: - - MnUserTransformation() : fPrecision(MnMachinePrecision()), - fParameters(std::vector()), - fExtOfInt(std::vector()), - fDoubleLimTrafo(SinParameterTransformation()), - fUpperLimTrafo(SqrtUpParameterTransformation()), - fLowerLimTrafo(SqrtLowParameterTransformation()), - fCache(std::vector()) {} - - MnUserTransformation(const std::vector&, const std::vector&); - - ~MnUserTransformation() {} - - MnUserTransformation(const MnUserTransformation& trafo) : - fPrecision(trafo.fPrecision), - fParameters(trafo.fParameters),fExtOfInt(trafo.fExtOfInt), - fDoubleLimTrafo(trafo.fDoubleLimTrafo), - fUpperLimTrafo(trafo.fUpperLimTrafo), - fLowerLimTrafo(trafo.fLowerLimTrafo), fCache(trafo.fCache) {} - - MnUserTransformation& operator=(const MnUserTransformation& trafo) { - if(this != &trafo) { - fPrecision = trafo.fPrecision; - fParameters = trafo.fParameters; - fExtOfInt = trafo.fExtOfInt; - fDoubleLimTrafo = trafo.fDoubleLimTrafo; - fUpperLimTrafo = trafo.fUpperLimTrafo; - fLowerLimTrafo = trafo.fLowerLimTrafo; - fCache = trafo.fCache; + class MnUserTransformation { + + public: + + MnUserTransformation() : fPrecision(MnMachinePrecision()), + fParameters(std::vector()), + fExtOfInt(std::vector()), + fDoubleLimTrafo(SinParameterTransformation()), + fUpperLimTrafo(SqrtUpParameterTransformation()), + fLowerLimTrafo(SqrtLowParameterTransformation()), + fCache(std::vector()) {} + + MnUserTransformation(const std::vector&, const std::vector&); + + ~MnUserTransformation() {} + + MnUserTransformation(const MnUserTransformation& trafo) : + fPrecision(trafo.fPrecision), + fParameters(trafo.fParameters),fExtOfInt(trafo.fExtOfInt), + fDoubleLimTrafo(trafo.fDoubleLimTrafo), + fUpperLimTrafo(trafo.fUpperLimTrafo), + fLowerLimTrafo(trafo.fLowerLimTrafo), fCache(trafo.fCache) {} + + MnUserTransformation& operator=(const MnUserTransformation& trafo) { + if(this != &trafo) { + fPrecision = trafo.fPrecision; + fParameters = trafo.fParameters; + fExtOfInt = trafo.fExtOfInt; + fDoubleLimTrafo = trafo.fDoubleLimTrafo; + fUpperLimTrafo = trafo.fUpperLimTrafo; + fLowerLimTrafo = trafo.fLowerLimTrafo; + fCache = trafo.fCache; + } + return *this; } - return *this; - } //#ifdef MINUIT2_THREAD_SAFE - // thread-safe version (do not use cache) - std::vector operator()(const MnAlgebraicVector&) const; + // thread-safe version (do not use cache) + std::vector operator()(const MnAlgebraicVector&) const; //#else // not thread safe // const std::vector& operator()(const MnAlgebraicVector&) const; //#endif - // Index = internal Parameter - double Int2ext(unsigned int, double) const; + // Index = internal Parameter + double Int2ext(unsigned int, double) const; - // Index = internal Parameter - double Int2extError(unsigned int, double, double) const; + // Index = internal Parameter + double Int2extError(unsigned int, double, double) const; - MnUserCovariance Int2extCovariance(const MnAlgebraicVector&, const MnAlgebraicSymMatrix&) const; + MnUserCovariance Int2extCovariance(const MnAlgebraicVector&, const MnAlgebraicSymMatrix&) const; - // Index = external Parameter - double Ext2int(unsigned int, double) const; + // Index = external Parameter + double Ext2int(unsigned int, double) const; - // Index = internal Parameter - double DInt2Ext(unsigned int, double) const; + // Index = internal Parameter + double DInt2Ext(unsigned int, double) const; + double D2Int2Ext(unsigned int, double) const; + double GStepInt2Ext(unsigned int, double) const; // // Index = external Parameter // double dExt2Int(unsigned int, double) const; - // Index = external Parameter - unsigned int IntOfExt(unsigned int) const; + // Index = external Parameter + unsigned int IntOfExt(unsigned int) const; - // Index = internal Parameter - unsigned int ExtOfInt(unsigned int internal) const { - assert(internal < fExtOfInt.size()); - return fExtOfInt[internal]; - } + // Index = internal Parameter + unsigned int ExtOfInt(unsigned int internal) const { + assert(internal < fExtOfInt.size()); + return fExtOfInt[internal]; + } - const std::vector& Parameters() const { - return fParameters; - } + const std::vector& Parameters() const { + return fParameters; + } - unsigned int VariableParameters() const {return static_cast ( fExtOfInt.size() );} + unsigned int VariableParameters() const {return static_cast ( fExtOfInt.size() );} - // return initial parameter values (useful especially to get fixed parameter values) - const std::vector & InitialParValues() const { - return fCache; - } + // return initial parameter values (useful especially to get fixed parameter values) + const std::vector & InitialParValues() const { + return fCache; + } - /** forwarded interface */ + /** forwarded interface */ - const MnMachinePrecision& Precision() const {return fPrecision;} - void SetPrecision(double eps) {fPrecision.SetPrecision(eps);} + const MnMachinePrecision& Precision() const {return fPrecision;} + void SetPrecision(double eps) {fPrecision.SetPrecision(eps);} - /// access to parameters and errors in column-wise representation + /// access to parameters and errors in column-wise representation - std::vector Params() const; - std::vector Errors() const; + std::vector Params() const; + std::vector Errors() const; - //access to single Parameter - const MinuitParameter& Parameter(unsigned int) const; + //access to single Parameter + const MinuitParameter& Parameter(unsigned int) const; - //add free Parameter - bool Add(const std::string &, double, double); - //add limited Parameter - bool Add(const std::string &, double, double, double, double); - //add const Parameter - bool Add(const std::string &, double); + //add free Parameter + bool Add(const std::string &, double, double); + //add limited Parameter + bool Add(const std::string &, double, double, double, double); + //add const Parameter + bool Add(const std::string &, double); - //interaction via external number of Parameter - void Fix(unsigned int); - void Release(unsigned int); - void RemoveLimits(unsigned int); - void SetValue(unsigned int, double); - void SetError(unsigned int, double); - void SetLimits(unsigned int, double, double); - void SetUpperLimit(unsigned int, double); - void SetLowerLimit(unsigned int, double); - void SetName(unsigned int, const std::string &); + //interaction via external number of Parameter + void Fix(unsigned int); + void Release(unsigned int); + void RemoveLimits(unsigned int); + void SetValue(unsigned int, double); + void SetError(unsigned int, double); + void SetLimits(unsigned int, double, double); + void SetUpperLimit(unsigned int, double); + void SetLowerLimit(unsigned int, double); + void SetName(unsigned int, const std::string &); - double Value(unsigned int) const; - double Error(unsigned int) const; + double Value(unsigned int) const; + double Error(unsigned int) const; - //interaction via Name of Parameter - void Fix(const std::string &); - void Release(const std::string &); - void SetValue(const std::string &, double); - void SetError(const std::string &, double); - void SetLimits(const std::string &, double, double); - void SetUpperLimit(const std::string &, double); - void SetLowerLimit(const std::string &, double); - void RemoveLimits(const std::string &); + //interaction via Name of Parameter + void Fix(const std::string &); + void Release(const std::string &); + void SetValue(const std::string &, double); + void SetError(const std::string &, double); + void SetLimits(const std::string &, double, double); + void SetUpperLimit(const std::string &, double); + void SetLowerLimit(const std::string &, double); + void RemoveLimits(const std::string &); - double Value(const std::string &) const; - double Error(const std::string &) const; + double Value(const std::string &) const; + double Error(const std::string &) const; - //convert Name into external number of Parameter (will assert if parameter is not found) - unsigned int Index(const std::string &) const; - // find parameter index given a name. If it is not found return a -1 - int FindIndex(const std::string & ) const; + //convert Name into external number of Parameter (will assert if parameter is not found) + unsigned int Index(const std::string &) const; + // find parameter index given a name. If it is not found return a -1 + int FindIndex(const std::string & ) const; - //convert external number into Name of Parameter (will assert if index is out of range) - const std::string & GetName(unsigned int) const; - // mantain interface with const char * for backward compatibility - const char* Name(unsigned int) const; + //convert external number into Name of Parameter (will assert if index is out of range) + const std::string & GetName(unsigned int) const; + // mantain interface with const char * for backward compatibility + const char* Name(unsigned int) const; -private: + private: - MnMachinePrecision fPrecision; + MnMachinePrecision fPrecision; - std::vector fParameters; - std::vector fExtOfInt; + std::vector fParameters; + std::vector fExtOfInt; - SinParameterTransformation fDoubleLimTrafo; - SqrtUpParameterTransformation fUpperLimTrafo; - SqrtLowParameterTransformation fLowerLimTrafo; + SinParameterTransformation fDoubleLimTrafo; + SqrtUpParameterTransformation fUpperLimTrafo; + SqrtLowParameterTransformation fLowerLimTrafo; - mutable std::vector fCache; + mutable std::vector fCache; -}; + }; - } // namespace Minuit2 + } // namespace Minuit2 } // namespace ROOT diff --git a/math/minuit2/inc/Minuit2/SinParameterTransformation.h b/math/minuit2/inc/Minuit2/SinParameterTransformation.h index a6924dfe1d79d..da6aadae560ed 100644 --- a/math/minuit2/inc/Minuit2/SinParameterTransformation.h +++ b/math/minuit2/inc/Minuit2/SinParameterTransformation.h @@ -1,9 +1,10 @@ // @(#)root/minuit2:$Id$ -// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005 +// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei, E.G.P. Bos 2003-2017 /********************************************************************** * * * Copyright (c) 2005 LCG ROOT Math team, CERN/PH-SFT * + * Copyright (c) 2017 Patrick Bos, Netherlands eScience Center * * * **********************************************************************/ @@ -12,32 +13,35 @@ namespace ROOT { - namespace Minuit2 { + namespace Minuit2 { -class MnMachinePrecision; + class MnMachinePrecision; /** class for the transformation for double-limited parameter Using a sin function one goes from a double-limited parameter range to an unlimited one */ -class SinParameterTransformation { + class SinParameterTransformation { -public: + public: - SinParameterTransformation() {} + SinParameterTransformation() {} - ~SinParameterTransformation() {} + ~SinParameterTransformation() {} - double Int2ext(double Value, double Upper, double Lower) const; - double Ext2int(double Value, double Upper, double Lower, - const MnMachinePrecision&) const; - double DInt2Ext(double Value, double Upper, double Lower) const; + double Int2ext(double Value, double Upper, double Lower) const; + double Ext2int(double Value, double Upper, double Lower, + const MnMachinePrecision&) const; + double DInt2Ext(double Value, double Upper, double Lower) const; -private: + double D2Int2Ext(double Value, double Upper, double Lower) const; + double GStepInt2Ext(double Value, double Upper, double Lower) const; -}; + private: + + }; } // namespace Minuit2 diff --git a/math/minuit2/inc/Minuit2/SqrtLowParameterTransformation.h b/math/minuit2/inc/Minuit2/SqrtLowParameterTransformation.h index 61571d4b3f13b..7ccbbe1af89ab 100644 --- a/math/minuit2/inc/Minuit2/SqrtLowParameterTransformation.h +++ b/math/minuit2/inc/Minuit2/SqrtLowParameterTransformation.h @@ -1,9 +1,10 @@ // @(#)root/minuit2:$Id$ -// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005 +// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei, E.G.P. Bos 2003-2017 /********************************************************************** * * * Copyright (c) 2005 LCG ROOT Math team, CERN/PH-SFT * + * Copyright (c) 2017 Patrick Bos, Netherlands eScience Center * * * **********************************************************************/ @@ -18,9 +19,9 @@ namespace ROOT { - namespace Minuit2 { + namespace Minuit2 { -class MnMachinePrecision; + class MnMachinePrecision; /** @@ -29,25 +30,28 @@ class MnMachinePrecision; * This transformation applies for the case of single side Lower Parameter limits */ -class SqrtLowParameterTransformation /* : public ParameterTransformation */ { + class SqrtLowParameterTransformation /* : public ParameterTransformation */ { -public: + public: - SqrtLowParameterTransformation() {} + SqrtLowParameterTransformation() {} - ~SqrtLowParameterTransformation() {} + ~SqrtLowParameterTransformation() {} - // transformation from internal to external - double Int2ext(double Value, double Lower) const; + // transformation from internal to external + double Int2ext(double Value, double Lower) const; - // transformation from external to internal - double Ext2int(double Value, double Lower, const MnMachinePrecision&) const; + // transformation from external to internal + double Ext2int(double Value, double Lower, const MnMachinePrecision&) const; - // derivative of transformation from internal to external - double DInt2Ext(double Value, double Lower) const; + // derivative of transformation from internal to external + double DInt2Ext(double Value, double Lower) const; -private: -}; + double D2Int2Ext(double Value, double Lower) const; + double GStepInt2Ext(double Value, double Lower) const; + + private: + }; } // namespace Minuit2 diff --git a/math/minuit2/inc/Minuit2/SqrtUpParameterTransformation.h b/math/minuit2/inc/Minuit2/SqrtUpParameterTransformation.h index ff127ffed5e4d..b82e193effaeb 100644 --- a/math/minuit2/inc/Minuit2/SqrtUpParameterTransformation.h +++ b/math/minuit2/inc/Minuit2/SqrtUpParameterTransformation.h @@ -1,9 +1,10 @@ // @(#)root/minuit2:$Id$ -// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005 +// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei, E.G.P. Bos 2003-2017 /********************************************************************** * * * Copyright (c) 2005 LCG ROOT Math team, CERN/PH-SFT * + * Copyright (c) 2017 Patrick Bos, Netherlands eScience Center * * * **********************************************************************/ @@ -19,10 +20,10 @@ namespace ROOT { - namespace Minuit2 { + namespace Minuit2 { -class MnMachinePrecision; + class MnMachinePrecision; /** @@ -31,27 +32,30 @@ class MnMachinePrecision; * This transformation applies for the case of single side Upper Parameter limits */ -class SqrtUpParameterTransformation /* : public ParameterTransformation */ { + class SqrtUpParameterTransformation /* : public ParameterTransformation */ { -public: + public: - // create with user defined precision - SqrtUpParameterTransformation() {} + // create with user defined precision + SqrtUpParameterTransformation() {} - ~SqrtUpParameterTransformation() {} + ~SqrtUpParameterTransformation() {} - // transformation from internal to external - double Int2ext(double Value, double Upper) const; + // transformation from internal to external + double Int2ext(double Value, double Upper) const; - // transformation from external to internal - double Ext2int(double Value, double Upper, const MnMachinePrecision&) const; + // transformation from external to internal + double Ext2int(double Value, double Upper, const MnMachinePrecision&) const; - // derivative of transformation from internal to external - double DInt2Ext(double Value, double Upper) const; + // derivative of transformation from internal to external + double DInt2Ext(double Value, double Upper) const; -private: + double D2Int2Ext(double Value, double Upper) const; + double GStepInt2Ext(double Value, double Upper) const; -}; + private: + + }; } // namespace Minuit2 diff --git a/math/minuit2/src/AnalyticalGradientCalculator.cxx b/math/minuit2/src/AnalyticalGradientCalculator.cxx index 2fa69b5d74589..186d03ad28a94 100644 --- a/math/minuit2/src/AnalyticalGradientCalculator.cxx +++ b/math/minuit2/src/AnalyticalGradientCalculator.cxx @@ -1,11 +1,10 @@ // @(#)root/minuit2:$Id$ -// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005 -// E.G.P. Bos 2017 +// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei, E.G.P. Bos 2003-2017 /********************************************************************** * * * Copyright (c) 2005 LCG ROOT Math team, CERN/PH-SFT * - * Copyright (c) 2017 Patrick Bos, NL eScience Center * + * Copyright (c) 2017 Patrick Bos, Netherlands eScience Center * * * **********************************************************************/ @@ -59,10 +58,10 @@ namespace ROOT { for(unsigned int i = 0; i < par.Vec().size(); i++) { unsigned int ext = fTransformation.ExtOfInt(i); if(fTransformation.Parameter(ext).HasLimits()) { - double jacobian_g2 = fTransformation.D2Int2Ext(i, par.Vec()(i)); - double jacobian_gstep = fTransformation.GStepInt2Ext(i, par.Vec()(i)); - vg2(i) = jacobian_g2 * g2[ext]; - vgstep(i) = jacobian_gstep * gstep[ext]; + double int2ext_g2 = fTransformation.D2Int2Ext(i, par.Vec()(i)); + double int2ext_gstep = fTransformation.GStepInt2Ext(i, par.Vec()(i)); + vg2(i) = int2ext_g2 * g2[ext]; + vgstep(i) = int2ext_gstep * gstep[ext]; } else { vg2(i) = g2[ext]; vgstep(i) = gstep[ext]; diff --git a/math/minuit2/src/InitialGradientCalculator.cxx b/math/minuit2/src/InitialGradientCalculator.cxx index 4a28dd9653f36..112130f5647b7 100644 --- a/math/minuit2/src/InitialGradientCalculator.cxx +++ b/math/minuit2/src/InitialGradientCalculator.cxx @@ -33,6 +33,8 @@ namespace ROOT { FunctionGradient InitialGradientCalculator::operator()(const MinimumParameters& par) const { // initial rough estimate of the gradient using the parameter step size + std::cout << "########### InitialGradientCalculator::operator()" < prec.Eps2() ? 1./dgrad.G2()(i) : 1.); } @@ -92,7 +99,12 @@ MinimumSeed MnSeedGenerator::operator()(const MnFcn& fcn, const GradientCalculat } NegativeG2LineSearch ng2ls; - if(ng2ls.HasNegativeG2(dgrad, prec)) { + std::cout << dgrad.Vec() << std::endl; + std::cout << prec << std::endl; + std::cout << ng2ls.HasNegativeG2(dgrad, prec) << std::endl; + std::cout << dgrad.Vec().size() << std::endl; + + if(ng2ls.HasNegativeG2(dgrad, prec)) { #ifdef DEBUG std::cout << "MnSeedGenerator: Negative G2 Found: " << std::endl; std::cout << x << std::endl; @@ -127,20 +139,33 @@ MinimumSeed MnSeedGenerator::operator()(const MnFcn& fcn, const GradientCalculat MinimumSeed MnSeedGenerator::operator()(const MnFcn& fcn, const AnalyticalGradientCalculator& gc, const MnUserParameterState& st, const MnStrategy& stra) const { - // find seed (initial point for minimization) using analytical gradient + + std::cout << "MnSeedGenerator::operator() for AnalyticalGradientCalculator" << std::endl; + + // find seed (initial point for minimization) using analytical gradient unsigned int n = st.VariableParameters(); const MnMachinePrecision& prec = st.Precision(); - // initial starting values + int printLevel = MnPrint::Level(); + + // initial starting values MnAlgebraicVector x(n); for(unsigned int i = 0; i < n; i++) x(i) = st.IntParameters()[i]; double fcnmin = fcn(x); MinimumParameters pa(x, fcnmin); - InitialGradientCalculator igc(fcn, st.Trafo(), stra); - FunctionGradient tmp = igc(pa); +// std::cout << "... creating igc ..." << std::endl; +// InitialGradientCalculator igc(fcn, st.Trafo(), stra); +// std::cout << "... doing igc(pa) ..." << std::endl; +// FunctionGradient tmp = igc(pa); +// std::cout << "tmp.G2: " << tmp.G2() << std::endl; FunctionGradient grd = gc(pa); - FunctionGradient dgrad(grd.Grad(), tmp.G2(), tmp.Gstep()); + std::cout << "grd.G2: " << grd.G2() << std::endl; + +// FunctionGradient dgrad(grd.Grad(), tmp.G2(), tmp.Gstep()); + FunctionGradient dgrad(grd.Grad(), grd.G2(), grd.Gstep()); + std::cout << "dgrad.Vec: " << dgrad.Vec() << std::endl; + std::cout << "dgrad.G2: " << dgrad.G2() << std::endl; if(gc.CheckGradient()) { bool good = true; @@ -170,10 +195,12 @@ MinimumSeed MnSeedGenerator::operator()(const MnFcn& fcn, const AnalyticalGradie MnAlgebraicSymMatrix mat(n); double dcovar = 1.; if(st.HasCovariance()) { + std::cout << "has covariance" << std::endl; for(unsigned int i = 0; i < n; i++) for(unsigned int j = i; j < n; j++) mat(i,j) = st.IntCovariance()(i,j); dcovar = 0.; } else { + std::cout << "has no covariance" << std::endl; for(unsigned int i = 0; i < n; i++) mat(i,i) = (fabs(dgrad.G2()(i)) > prec.Eps2() ? 1./dgrad.G2()(i) : 1.); } @@ -182,9 +209,20 @@ MinimumSeed MnSeedGenerator::operator()(const MnFcn& fcn, const AnalyticalGradie MinimumState state(pa, err, dgrad, edm, fcn.NumOfCalls()); NegativeG2LineSearch ng2ls; + std::cout << dgrad.Vec() << std::endl; + std::cout << prec << std::endl; + std::cout << ng2ls.HasNegativeG2(dgrad, prec) << std::endl; + std::cout << dgrad.Vec().size() << std::endl; + if(ng2ls.HasNegativeG2(dgrad, prec)) { - Numerical2PGradientCalculator ngc(fcn, st.Trafo(), stra); - state = ng2ls(fcn, state, ngc, prec); +// Numerical2PGradientCalculator ngc(fcn, st.Trafo(), stra); +// state = ng2ls(fcn, state, ngc, prec); + state = ng2ls(fcn, state, gc, prec); + + if (printLevel >1) { + MnPrint::PrintState(std::cout, state, "MnSeedGenerator: Negative G2 found - new state: "); + } + } if(stra.Strategy() == 2 && !st.HasCovariance()) { diff --git a/math/minuit2/src/MnUserTransformation.cxx b/math/minuit2/src/MnUserTransformation.cxx index d469cfde79035..82658c51a067f 100644 --- a/math/minuit2/src/MnUserTransformation.cxx +++ b/math/minuit2/src/MnUserTransformation.cxx @@ -1,9 +1,10 @@ // @(#)root/minuit2:$Id$ -// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005 +// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei, E.G.P. Bos 2003-2017 /********************************************************************** * * * Copyright (c) 2005 LCG ROOT Math team, CERN/PH-SFT * + * Copyright (c) 2017 Patrick Bos, Netherlands eScience Center * * * **********************************************************************/ @@ -18,63 +19,63 @@ namespace ROOT { - namespace Minuit2 { + namespace Minuit2 { -class MnParStr { + class MnParStr { -public: + public: - MnParStr(const std::string & name) : fName(name) {} + MnParStr(const std::string & name) : fName(name) {} - ~MnParStr() {} + ~MnParStr() {} - bool operator()(const MinuitParameter& par) const { + bool operator()(const MinuitParameter& par) const { // return (strcmp(par.Name(), fName) == 0); - return par.GetName() == fName; - } + return par.GetName() == fName; + } -private: - const std::string & fName; -}; + private: + const std::string & fName; + }; -MnUserTransformation::MnUserTransformation(const std::vector& par, const std::vector& err) : fPrecision(MnMachinePrecision()), fParameters(std::vector()), fExtOfInt(std::vector()), fDoubleLimTrafo(SinParameterTransformation()),fUpperLimTrafo(SqrtUpParameterTransformation()), fLowerLimTrafo(SqrtLowParameterTransformation()), fCache(std::vector()) { - // constructor from a vector of parameter values and a vector of errors (step sizes) - // class has as data member the transformation objects (all of the types), - // the std::vector of MinuitParameter objects and the vector with the index conversions from - // internal to external (fExtOfInt) + MnUserTransformation::MnUserTransformation(const std::vector& par, const std::vector& err) : fPrecision(MnMachinePrecision()), fParameters(std::vector()), fExtOfInt(std::vector()), fDoubleLimTrafo(SinParameterTransformation()),fUpperLimTrafo(SqrtUpParameterTransformation()), fLowerLimTrafo(SqrtLowParameterTransformation()), fCache(std::vector()) { + // constructor from a vector of parameter values and a vector of errors (step sizes) + // class has as data member the transformation objects (all of the types), + // the std::vector of MinuitParameter objects and the vector with the index conversions from + // internal to external (fExtOfInt) - fParameters.reserve(par.size()); - fExtOfInt.reserve(par.size()); - fCache.reserve(par.size()); + fParameters.reserve(par.size()); + fExtOfInt.reserve(par.size()); + fCache.reserve(par.size()); - std::string parName; - for(unsigned int i = 0; i < par.size(); i++) { - std::ostringstream buf; - buf << "p" << i; - parName = buf.str(); - Add(parName, par[i], err[i]); - } -} + std::string parName; + for(unsigned int i = 0; i < par.size(); i++) { + std::ostringstream buf; + buf << "p" << i; + parName = buf.str(); + Add(parName, par[i], err[i]); + } + } //#ifdef MINUIT2_THREAD_SAFE // this if a thread-safe implementation needed if want to share transformation object between the threads -std::vector MnUserTransformation::operator()(const MnAlgebraicVector& pstates ) const { - // transform an internal Minuit vector of internal values in a std::vector of external values - // fixed parameters will have their fixed values - unsigned int n = pstates.size(); - // need to initialize to the stored (initial values) parameter values for the fixed ones - std::vector pcache( fCache ); - for(unsigned int i = 0; i < n; i++) { - if(fParameters[fExtOfInt[i]].HasLimits()) { - pcache[fExtOfInt[i]] = Int2ext(i, pstates(i)); - } else { - pcache[fExtOfInt[i]] = pstates(i); + std::vector MnUserTransformation::operator()(const MnAlgebraicVector& pstates ) const { + // transform an internal Minuit vector of internal values in a std::vector of external values + // fixed parameters will have their fixed values + unsigned int n = pstates.size(); + // need to initialize to the stored (initial values) parameter values for the fixed ones + std::vector pcache( fCache ); + for(unsigned int i = 0; i < n; i++) { + if(fParameters[fExtOfInt[i]].HasLimits()) { + pcache[fExtOfInt[i]] = Int2ext(i, pstates(i)); + } else { + pcache[fExtOfInt[i]] = pstates(i); + } } - } - return pcache; -} + return pcache; + } // #else // const std::vector & MnUserTransformation::operator()(const MnAlgebraicVector& pstates) const { @@ -92,107 +93,144 @@ std::vector MnUserTransformation::operator()(const MnAlgebraicVector& ps // } // #endif -double MnUserTransformation::Int2ext(unsigned int i, double val) const { - // return external value from internal value for parameter i - if(fParameters[fExtOfInt[i]].HasLimits()) { - if(fParameters[fExtOfInt[i]].HasUpperLimit() && fParameters[fExtOfInt[i]].HasLowerLimit()) - return fDoubleLimTrafo.Int2ext(val, fParameters[fExtOfInt[i]].UpperLimit(), fParameters[fExtOfInt[i]].LowerLimit()); - else if(fParameters[fExtOfInt[i]].HasUpperLimit() && !fParameters[fExtOfInt[i]].HasLowerLimit()) - return fUpperLimTrafo.Int2ext(val, fParameters[fExtOfInt[i]].UpperLimit()); - else - return fLowerLimTrafo.Int2ext(val, fParameters[fExtOfInt[i]].LowerLimit()); - } - - return val; -} - -double MnUserTransformation::Int2extError(unsigned int i, double val, double err) const { - // return external error from internal error for parameter i - - //err = sigma Value == sqrt(cov(i,i)) - double dx = err; - - if(fParameters[fExtOfInt[i]].HasLimits()) { - double ui = Int2ext(i, val); - double du1 = Int2ext(i, val+dx) - ui; - double du2 = Int2ext(i, val-dx) - ui; - if(fParameters[fExtOfInt[i]].HasUpperLimit() && fParameters[fExtOfInt[i]].HasLowerLimit()) { - // double al = fParameters[fExtOfInt[i]].Lower(); - // double ba = fParameters[fExtOfInt[i]].Upper() - al; - // double du1 = al + 0.5*(sin(val + dx) + 1.)*ba - ui; - // double du2 = al + 0.5*(sin(val - dx) + 1.)*ba - ui; - // if(dx > 1.) du1 = ba; - if(dx > 1.) du1 = fParameters[fExtOfInt[i]].UpperLimit() - fParameters[fExtOfInt[i]].LowerLimit(); - dx = 0.5*(fabs(du1) + fabs(du2)); - } else { - dx = 0.5*(fabs(du1) + fabs(du2)); + double MnUserTransformation::Int2ext(unsigned int i, double val) const { + // return external value from internal value for parameter i + if(fParameters[fExtOfInt[i]].HasLimits()) { + if(fParameters[fExtOfInt[i]].HasUpperLimit() && fParameters[fExtOfInt[i]].HasLowerLimit()) + return fDoubleLimTrafo.Int2ext(val, fParameters[fExtOfInt[i]].UpperLimit(), fParameters[fExtOfInt[i]].LowerLimit()); + else if(fParameters[fExtOfInt[i]].HasUpperLimit() && !fParameters[fExtOfInt[i]].HasLowerLimit()) + return fUpperLimTrafo.Int2ext(val, fParameters[fExtOfInt[i]].UpperLimit()); + else + return fLowerLimTrafo.Int2ext(val, fParameters[fExtOfInt[i]].LowerLimit()); } - } - return dx; -} + return val; + } + + double MnUserTransformation::Int2extError(unsigned int i, double val, double err) const { + // return external error from internal error for parameter i -MnUserCovariance MnUserTransformation::Int2extCovariance(const MnAlgebraicVector& vec, const MnAlgebraicSymMatrix& cov) const { - // return the external covariance matrix from the internal error matrix and the internal parameter value - // the vector of internal parameter is needed for the derivatives (Jacobian of the transformation) - // Vext(i,j) = Vint(i,j) * dPext(i)/dPint(i) * dPext(j)/dPint(j) + //err = sigma Value == sqrt(cov(i,i)) + double dx = err; - MnUserCovariance result(cov.Nrow()); - for(unsigned int i = 0; i < vec.size(); i++) { - double dxdi = 1.; if(fParameters[fExtOfInt[i]].HasLimits()) { - // dxdi = 0.5*fabs((fParameters[fExtOfInt[i]].Upper() - fParameters[fExtOfInt[i]].Lower())*cos(vec(i))); - dxdi = DInt2Ext(i, vec(i)); + double ui = Int2ext(i, val); + double du1 = Int2ext(i, val+dx) - ui; + double du2 = Int2ext(i, val-dx) - ui; + if(fParameters[fExtOfInt[i]].HasUpperLimit() && fParameters[fExtOfInt[i]].HasLowerLimit()) { + // double al = fParameters[fExtOfInt[i]].Lower(); + // double ba = fParameters[fExtOfInt[i]].Upper() - al; + // double du1 = al + 0.5*(sin(val + dx) + 1.)*ba - ui; + // double du2 = al + 0.5*(sin(val - dx) + 1.)*ba - ui; + // if(dx > 1.) du1 = ba; + if(dx > 1.) du1 = fParameters[fExtOfInt[i]].UpperLimit() - fParameters[fExtOfInt[i]].LowerLimit(); + dx = 0.5*(fabs(du1) + fabs(du2)); + } else { + dx = 0.5*(fabs(du1) + fabs(du2)); + } } - for(unsigned int j = i; j < vec.size(); j++) { - double dxdj = 1.; - if(fParameters[fExtOfInt[j]].HasLimits()) { + + return dx; + } + + MnUserCovariance MnUserTransformation::Int2extCovariance(const MnAlgebraicVector& vec, const MnAlgebraicSymMatrix& cov) const { + // return the external covariance matrix from the internal error matrix and the internal parameter value + // the vector of internal parameter is needed for the derivatives (Jacobian of the transformation) + // Vext(i,j) = Vint(i,j) * dPext(i)/dPint(i) * dPext(j)/dPint(j) + + MnUserCovariance result(cov.Nrow()); + for(unsigned int i = 0; i < vec.size(); i++) { + double dxdi = 1.; + if(fParameters[fExtOfInt[i]].HasLimits()) { + // dxdi = 0.5*fabs((fParameters[fExtOfInt[i]].Upper() - fParameters[fExtOfInt[i]].Lower())*cos(vec(i))); + dxdi = DInt2Ext(i, vec(i)); + } + for(unsigned int j = i; j < vec.size(); j++) { + double dxdj = 1.; + if(fParameters[fExtOfInt[j]].HasLimits()) { // dxdj = 0.5*fabs((fParameters[fExtOfInt[j]].Upper() - fParameters[fExtOfInt[j]].Lower())*cos(vec(j))); dxdj = DInt2Ext(j, vec(j)); - } - result(i,j) = dxdi*cov(i,j)*dxdj; + } + result(i,j) = dxdi*cov(i,j)*dxdj; + } + // double diag = Int2extError(i, vec(i), sqrt(cov(i,i))); + // result(i,i) = diag*diag; + } + + return result; + } + + double MnUserTransformation::Ext2int(unsigned int i, double val) const { + // return the internal value for parameter i with external value val + + if(fParameters[i].HasLimits()) { + if(fParameters[i].HasUpperLimit() && fParameters[i].HasLowerLimit()) + return fDoubleLimTrafo.Ext2int(val, fParameters[i].UpperLimit(), fParameters[i].LowerLimit(), Precision()); + else if(fParameters[i].HasUpperLimit() && !fParameters[i].HasLowerLimit()) + return fUpperLimTrafo.Ext2int(val, fParameters[i].UpperLimit(), Precision()); + else + return fLowerLimTrafo.Ext2int(val, fParameters[i].LowerLimit(), Precision()); + } + + return val; + } + + double MnUserTransformation::DInt2Ext(unsigned int i, double val) const { + // return the derivative of the int->ext transformation: dPext(i) / dPint(i) + // for the parameter i with value val + + double dd = 1.; + if(fParameters[fExtOfInt[i]].HasLimits()) { + if(fParameters[fExtOfInt[i]].HasUpperLimit() && fParameters[fExtOfInt[i]].HasLowerLimit()) + // dd = 0.5*fabs((fParameters[fExtOfInt[i]].Upper() - fParameters[fExtOfInt[i]].Lower())*cos(vec(i))); + dd = fDoubleLimTrafo.DInt2Ext(val, fParameters[fExtOfInt[i]].UpperLimit(), fParameters[fExtOfInt[i]].LowerLimit()); + else if(fParameters[fExtOfInt[i]].HasUpperLimit() && !fParameters[fExtOfInt[i]].HasLowerLimit()) + dd = fUpperLimTrafo.DInt2Ext(val, fParameters[fExtOfInt[i]].UpperLimit()); + else + dd = fLowerLimTrafo.DInt2Ext(val, fParameters[fExtOfInt[i]].LowerLimit()); + } + + return dd; + } + + double MnUserTransformation::D2Int2Ext(unsigned int i, double val) const { + // return the 2nd derivative of the int->ext transformation: d^2{Pext(i)} / d{Pint(i)}^2 + // for the parameter i with value val + + double dd = 1.; + if(fParameters[fExtOfInt[i]].HasLimits()) { + if (fParameters[fExtOfInt[i]].HasUpperLimit() && fParameters[fExtOfInt[i]].HasLowerLimit()) { + dd = fDoubleLimTrafo.D2Int2Ext(val, fParameters[fExtOfInt[i]].UpperLimit(), + fParameters[fExtOfInt[i]].LowerLimit()); + } else if (fParameters[fExtOfInt[i]].HasUpperLimit() && !fParameters[fExtOfInt[i]].HasLowerLimit()) { + dd = fUpperLimTrafo.D2Int2Ext(val, fParameters[fExtOfInt[i]].UpperLimit()); + } else { + dd = fLowerLimTrafo.D2Int2Ext(val, fParameters[fExtOfInt[i]].LowerLimit()); + } + } + + return dd; + } + + double MnUserTransformation::GStepInt2Ext(unsigned int i, double val) const { + // return the conversion factor of the int->ext transformation for the step size + // for the parameter i with value val + + double dd = 1.; + if(fParameters[fExtOfInt[i]].HasLimits()) { + if(fParameters[fExtOfInt[i]].HasUpperLimit() && fParameters[fExtOfInt[i]].HasLowerLimit()) { + dd = fDoubleLimTrafo.GStepInt2Ext(val, fParameters[fExtOfInt[i]].UpperLimit(), fParameters[fExtOfInt[i]].LowerLimit()); + } else if(fParameters[fExtOfInt[i]].HasUpperLimit() && !fParameters[fExtOfInt[i]].HasLowerLimit()) { + dd = fUpperLimTrafo.GStepInt2Ext(val, fParameters[fExtOfInt[i]].UpperLimit()); + } else { + dd = fLowerLimTrafo.GStepInt2Ext(val, fParameters[fExtOfInt[i]].LowerLimit()); + } } - // double diag = Int2extError(i, vec(i), sqrt(cov(i,i))); - // result(i,i) = diag*diag; - } - - return result; -} - -double MnUserTransformation::Ext2int(unsigned int i, double val) const { - // return the internal value for parameter i with external value val - - if(fParameters[i].HasLimits()) { - if(fParameters[i].HasUpperLimit() && fParameters[i].HasLowerLimit()) - return fDoubleLimTrafo.Ext2int(val, fParameters[i].UpperLimit(), fParameters[i].LowerLimit(), Precision()); - else if(fParameters[i].HasUpperLimit() && !fParameters[i].HasLowerLimit()) - return fUpperLimTrafo.Ext2int(val, fParameters[i].UpperLimit(), Precision()); - else - return fLowerLimTrafo.Ext2int(val, fParameters[i].LowerLimit(), Precision()); - } - - return val; -} - -double MnUserTransformation::DInt2Ext(unsigned int i, double val) const { - // return the derivative of the int->ext transformation: dPext(i) / dPint(i) - // for the parameter i with value val - - double dd = 1.; - if(fParameters[fExtOfInt[i]].HasLimits()) { - if(fParameters[fExtOfInt[i]].HasUpperLimit() && fParameters[fExtOfInt[i]].HasLowerLimit()) - // dd = 0.5*fabs((fParameters[fExtOfInt[i]].Upper() - fParameters[fExtOfInt[i]].Lower())*cos(vec(i))); - dd = fDoubleLimTrafo.DInt2Ext(val, fParameters[fExtOfInt[i]].UpperLimit(), fParameters[fExtOfInt[i]].LowerLimit()); - else if(fParameters[fExtOfInt[i]].HasUpperLimit() && !fParameters[fExtOfInt[i]].HasLowerLimit()) - dd = fUpperLimTrafo.DInt2Ext(val, fParameters[fExtOfInt[i]].UpperLimit()); - else - dd = fLowerLimTrafo.DInt2Ext(val, fParameters[fExtOfInt[i]].LowerLimit()); - } - - return dd; -} - -/* + + return dd; + } + + /* double MnUserTransformation::dExt2Int(unsigned int, double) const { double dd = 1.; @@ -210,42 +248,42 @@ double MnUserTransformation::DInt2Ext(unsigned int i, double val) const { } */ -unsigned int MnUserTransformation::IntOfExt(unsigned int ext) const { - // return internal index given external one ext - assert(ext < fParameters.size()); - assert(!fParameters[ext].IsFixed()); - assert(!fParameters[ext].IsConst()); - std::vector::const_iterator iind = std::find(fExtOfInt.begin(), fExtOfInt.end(), ext); - assert(iind != fExtOfInt.end()); - - return (iind - fExtOfInt.begin()); -} - -std::vector MnUserTransformation::Params() const { - // return std::vector of double with parameter values - unsigned int n = fParameters.size(); - std::vector result(n); - for(unsigned int i = 0; i < n; ++i) - result[i] = fParameters[i].Value(); - - return result; -} - -std::vector MnUserTransformation::Errors() const { - // return std::vector of double with parameter errors - std::vector result; result.reserve(fParameters.size()); - for(std::vector::const_iterator ipar = Parameters().begin(); - ipar != Parameters().end(); ipar++) - result.push_back((*ipar).Error()); - - return result; -} - -const MinuitParameter& MnUserTransformation::Parameter(unsigned int n) const { - // return the MinuitParameter object for index n (external) - assert(n < fParameters.size()); - return fParameters[n]; -} + unsigned int MnUserTransformation::IntOfExt(unsigned int ext) const { + // return internal index given external one ext + assert(ext < fParameters.size()); + assert(!fParameters[ext].IsFixed()); + assert(!fParameters[ext].IsConst()); + std::vector::const_iterator iind = std::find(fExtOfInt.begin(), fExtOfInt.end(), ext); + assert(iind != fExtOfInt.end()); + + return (iind - fExtOfInt.begin()); + } + + std::vector MnUserTransformation::Params() const { + // return std::vector of double with parameter values + unsigned int n = fParameters.size(); + std::vector result(n); + for(unsigned int i = 0; i < n; ++i) + result[i] = fParameters[i].Value(); + + return result; + } + + std::vector MnUserTransformation::Errors() const { + // return std::vector of double with parameter errors + std::vector result; result.reserve(fParameters.size()); + for(std::vector::const_iterator ipar = Parameters().begin(); + ipar != Parameters().end(); ipar++) + result.push_back((*ipar).Error()); + + return result; + } + + const MinuitParameter& MnUserTransformation::Parameter(unsigned int n) const { + // return the MinuitParameter object for index n (external) + assert(n < fParameters.size()); + return fParameters[n]; + } // bool MnUserTransformation::Remove(const std::string & name) { // // remove parameter with name @@ -261,198 +299,198 @@ const MinuitParameter& MnUserTransformation::Parameter(unsigned int n) const { // if (iind != fExtOfInt.end()) fExtOfInt.erase(iind); // } -bool MnUserTransformation::Add(const std::string & name, double val, double err) { - // add a new unlimited parameter giving name, value and err (step size) - // return false if parameter already exists - if (std::find_if(fParameters.begin(), fParameters.end(), MnParStr(name)) != fParameters.end() ) - return false; - fExtOfInt.push_back(fParameters.size()); - fCache.push_back(val); - fParameters.push_back(MinuitParameter(fParameters.size(), name, val, err)); - return true; -} - -bool MnUserTransformation::Add(const std::string & name, double val, double err, double low, double up) { - // add a new limited parameter giving name, value, err (step size) and lower/upper limits - // return false if parameter already exists - if (std::find_if(fParameters.begin(), fParameters.end(), MnParStr(name)) != fParameters.end() ) - return false; - fExtOfInt.push_back(fParameters.size()); - fCache.push_back(val); - fParameters.push_back(MinuitParameter(fParameters.size(), name, val, err, low, up)); - return true; -} - -bool MnUserTransformation::Add(const std::string & name, double val) { - // add a new constant parameter giving name and value - // return false if parameter already exists - if (std::find_if(fParameters.begin(), fParameters.end(), MnParStr(name)) != fParameters.end() ) - return false; - fCache.push_back(val); - // costant parameter - do not add in list of internals (fExtOfInt) - fParameters.push_back(MinuitParameter(fParameters.size(), name, val)); - return true; -} - -void MnUserTransformation::Fix(unsigned int n) { - // fix parameter n (external index) - assert(n < fParameters.size()); - std::vector::iterator iind = std::find(fExtOfInt.begin(), fExtOfInt.end(), n); - if (iind != fExtOfInt.end()) - fExtOfInt.erase(iind, iind+1); - fParameters[n].Fix(); -} - -void MnUserTransformation::Release(unsigned int n) { - // release parameter n (external index) - assert(n < fParameters.size()); - std::vector::const_iterator iind = std::find(fExtOfInt.begin(), fExtOfInt.end(), n); - if (iind == fExtOfInt.end() ) { - fExtOfInt.push_back(n); - std::sort(fExtOfInt.begin(), fExtOfInt.end()); - } - fParameters[n].Release(); -} - -void MnUserTransformation::SetValue(unsigned int n, double val) { - // set value for parameter n (external index) - assert(n < fParameters.size()); - fParameters[n].SetValue(val); - fCache[n] = val; -} - -void MnUserTransformation::SetError(unsigned int n, double err) { - // set error for parameter n (external index) - assert(n < fParameters.size()); - fParameters[n].SetError(err); -} - -void MnUserTransformation::SetLimits(unsigned int n, double low, double up) { - // set limits (lower/upper) for parameter n (external index) - assert(n < fParameters.size()); - assert(low != up); - fParameters[n].SetLimits(low, up); -} - -void MnUserTransformation::SetUpperLimit(unsigned int n, double up) { - // set upper limit for parameter n (external index) - assert(n < fParameters.size()); - fParameters[n].SetUpperLimit(up); -} - -void MnUserTransformation::SetLowerLimit(unsigned int n, double lo) { - // set lower limit for parameter n (external index) - assert(n < fParameters.size()); - fParameters[n].SetLowerLimit(lo); -} - -void MnUserTransformation::RemoveLimits(unsigned int n) { - // remove limits for parameter n (external index) - assert(n < fParameters.size()); - fParameters[n].RemoveLimits(); -} - -void MnUserTransformation::SetName(unsigned int n, const std::string & name) { - // set name for parameter n (external index) - assert(n < fParameters.size()); - fParameters[n].SetName(name); -} - - -double MnUserTransformation::Value(unsigned int n) const { - // get value for parameter n (external index) - assert(n < fParameters.size()); - return fParameters[n].Value(); -} - -double MnUserTransformation::Error(unsigned int n) const { - // get error for parameter n (external index) - assert(n < fParameters.size()); - return fParameters[n].Error(); -} + bool MnUserTransformation::Add(const std::string & name, double val, double err) { + // add a new unlimited parameter giving name, value and err (step size) + // return false if parameter already exists + if (std::find_if(fParameters.begin(), fParameters.end(), MnParStr(name)) != fParameters.end() ) + return false; + fExtOfInt.push_back(fParameters.size()); + fCache.push_back(val); + fParameters.push_back(MinuitParameter(fParameters.size(), name, val, err)); + return true; + } + + bool MnUserTransformation::Add(const std::string & name, double val, double err, double low, double up) { + // add a new limited parameter giving name, value, err (step size) and lower/upper limits + // return false if parameter already exists + if (std::find_if(fParameters.begin(), fParameters.end(), MnParStr(name)) != fParameters.end() ) + return false; + fExtOfInt.push_back(fParameters.size()); + fCache.push_back(val); + fParameters.push_back(MinuitParameter(fParameters.size(), name, val, err, low, up)); + return true; + } + + bool MnUserTransformation::Add(const std::string & name, double val) { + // add a new constant parameter giving name and value + // return false if parameter already exists + if (std::find_if(fParameters.begin(), fParameters.end(), MnParStr(name)) != fParameters.end() ) + return false; + fCache.push_back(val); + // costant parameter - do not add in list of internals (fExtOfInt) + fParameters.push_back(MinuitParameter(fParameters.size(), name, val)); + return true; + } + + void MnUserTransformation::Fix(unsigned int n) { + // fix parameter n (external index) + assert(n < fParameters.size()); + std::vector::iterator iind = std::find(fExtOfInt.begin(), fExtOfInt.end(), n); + if (iind != fExtOfInt.end()) + fExtOfInt.erase(iind, iind+1); + fParameters[n].Fix(); + } + + void MnUserTransformation::Release(unsigned int n) { + // release parameter n (external index) + assert(n < fParameters.size()); + std::vector::const_iterator iind = std::find(fExtOfInt.begin(), fExtOfInt.end(), n); + if (iind == fExtOfInt.end() ) { + fExtOfInt.push_back(n); + std::sort(fExtOfInt.begin(), fExtOfInt.end()); + } + fParameters[n].Release(); + } + + void MnUserTransformation::SetValue(unsigned int n, double val) { + // set value for parameter n (external index) + assert(n < fParameters.size()); + fParameters[n].SetValue(val); + fCache[n] = val; + } + + void MnUserTransformation::SetError(unsigned int n, double err) { + // set error for parameter n (external index) + assert(n < fParameters.size()); + fParameters[n].SetError(err); + } + + void MnUserTransformation::SetLimits(unsigned int n, double low, double up) { + // set limits (lower/upper) for parameter n (external index) + assert(n < fParameters.size()); + assert(low != up); + fParameters[n].SetLimits(low, up); + } + + void MnUserTransformation::SetUpperLimit(unsigned int n, double up) { + // set upper limit for parameter n (external index) + assert(n < fParameters.size()); + fParameters[n].SetUpperLimit(up); + } + + void MnUserTransformation::SetLowerLimit(unsigned int n, double lo) { + // set lower limit for parameter n (external index) + assert(n < fParameters.size()); + fParameters[n].SetLowerLimit(lo); + } + + void MnUserTransformation::RemoveLimits(unsigned int n) { + // remove limits for parameter n (external index) + assert(n < fParameters.size()); + fParameters[n].RemoveLimits(); + } + + void MnUserTransformation::SetName(unsigned int n, const std::string & name) { + // set name for parameter n (external index) + assert(n < fParameters.size()); + fParameters[n].SetName(name); + } + + + double MnUserTransformation::Value(unsigned int n) const { + // get value for parameter n (external index) + assert(n < fParameters.size()); + return fParameters[n].Value(); + } + + double MnUserTransformation::Error(unsigned int n) const { + // get error for parameter n (external index) + assert(n < fParameters.size()); + return fParameters[n].Error(); + } // interface by parameter name -void MnUserTransformation::Fix(const std::string & name) { - // fix parameter - Fix(Index(name)); -} - -void MnUserTransformation::Release(const std::string & name) { - // release parameter - Release(Index(name)); -} - -void MnUserTransformation::SetValue(const std::string & name, double val) { - // set value for parameter - SetValue(Index(name), val); -} - -void MnUserTransformation::SetError(const std::string & name, double err) { - // set error - SetError(Index(name), err); -} - -void MnUserTransformation::SetLimits(const std::string & name, double low, double up) { - // set lower/upper limits - SetLimits(Index(name), low, up); -} - -void MnUserTransformation::SetUpperLimit(const std::string & name, double up) { - // set upper limit - SetUpperLimit(Index(name), up); -} - -void MnUserTransformation::SetLowerLimit(const std::string & name, double lo) { - // set lower limit - SetLowerLimit(Index(name), lo); -} - -void MnUserTransformation::RemoveLimits(const std::string & name) { - // remove limits - RemoveLimits(Index(name)); -} - -double MnUserTransformation::Value(const std::string & name) const { - // get parameter value - return Value(Index(name)); -} - -double MnUserTransformation::Error(const std::string & name) const { - // get parameter error - return Error(Index(name)); -} - -unsigned int MnUserTransformation::Index(const std::string & name) const { - // get index (external) corresponding to name - std::vector::const_iterator ipar = - std::find_if(fParameters.begin(), fParameters.end(), MnParStr(name)); - assert(ipar != fParameters.end()); - // return (ipar - fParameters.begin()); - return (*ipar).Number(); -} - -int MnUserTransformation::FindIndex(const std::string & name) const { - // find index (external) corresponding to name - return -1 if not found - std::vector::const_iterator ipar = - std::find_if(fParameters.begin(), fParameters.end(), MnParStr(name)); - if (ipar == fParameters.end() ) return -1; - return (*ipar).Number(); -} - - -const std::string & MnUserTransformation::GetName(unsigned int n) const { - // get name corresponding to index (external) - assert(n < fParameters.size()); - return fParameters[n].GetName(); -} - -const char* MnUserTransformation::Name(unsigned int n) const { - // get name corresponding to index (external) - return GetName(n).c_str(); -} - - - } // namespace Minuit2 + void MnUserTransformation::Fix(const std::string & name) { + // fix parameter + Fix(Index(name)); + } + + void MnUserTransformation::Release(const std::string & name) { + // release parameter + Release(Index(name)); + } + + void MnUserTransformation::SetValue(const std::string & name, double val) { + // set value for parameter + SetValue(Index(name), val); + } + + void MnUserTransformation::SetError(const std::string & name, double err) { + // set error + SetError(Index(name), err); + } + + void MnUserTransformation::SetLimits(const std::string & name, double low, double up) { + // set lower/upper limits + SetLimits(Index(name), low, up); + } + + void MnUserTransformation::SetUpperLimit(const std::string & name, double up) { + // set upper limit + SetUpperLimit(Index(name), up); + } + + void MnUserTransformation::SetLowerLimit(const std::string & name, double lo) { + // set lower limit + SetLowerLimit(Index(name), lo); + } + + void MnUserTransformation::RemoveLimits(const std::string & name) { + // remove limits + RemoveLimits(Index(name)); + } + + double MnUserTransformation::Value(const std::string & name) const { + // get parameter value + return Value(Index(name)); + } + + double MnUserTransformation::Error(const std::string & name) const { + // get parameter error + return Error(Index(name)); + } + + unsigned int MnUserTransformation::Index(const std::string & name) const { + // get index (external) corresponding to name + std::vector::const_iterator ipar = + std::find_if(fParameters.begin(), fParameters.end(), MnParStr(name)); + assert(ipar != fParameters.end()); + // return (ipar - fParameters.begin()); + return (*ipar).Number(); + } + + int MnUserTransformation::FindIndex(const std::string & name) const { + // find index (external) corresponding to name - return -1 if not found + std::vector::const_iterator ipar = + std::find_if(fParameters.begin(), fParameters.end(), MnParStr(name)); + if (ipar == fParameters.end() ) return -1; + return (*ipar).Number(); + } + + + const std::string & MnUserTransformation::GetName(unsigned int n) const { + // get name corresponding to index (external) + assert(n < fParameters.size()); + return fParameters[n].GetName(); + } + + const char* MnUserTransformation::Name(unsigned int n) const { + // get name corresponding to index (external) + return GetName(n).c_str(); + } + + + } // namespace Minuit2 } // namespace ROOT diff --git a/math/minuit2/src/NegativeG2LineSearch.cxx b/math/minuit2/src/NegativeG2LineSearch.cxx index dc36155bad66d..45c026b8d4e70 100644 --- a/math/minuit2/src/NegativeG2LineSearch.cxx +++ b/math/minuit2/src/NegativeG2LineSearch.cxx @@ -22,6 +22,9 @@ #include #endif + +#include + namespace ROOT { namespace Minuit2 { @@ -115,12 +118,15 @@ MinimumState NegativeG2LineSearch::operator()(const MnFcn& fcn, const MinimumSta bool NegativeG2LineSearch::HasNegativeG2(const FunctionGradient& grad, const MnMachinePrecision& /*prec */ ) const { // check if function gradient has any component which is neegative - for(unsigned int i = 0; i < grad.Vec().size(); i++) - if(grad.G2()(i) <= 0 ) { + for(unsigned int i = 0; i < grad.Vec().size(); i++) { + + std::cout << grad.G2()(i) << std::endl; + if (grad.G2()(i) <= 0) { - return true; - } + return true; + } + } return false; } diff --git a/math/minuit2/src/Numerical2PGradientCalculator.cxx b/math/minuit2/src/Numerical2PGradientCalculator.cxx index f3137d1575933..f97332d803570 100644 --- a/math/minuit2/src/Numerical2PGradientCalculator.cxx +++ b/math/minuit2/src/Numerical2PGradientCalculator.cxx @@ -81,6 +81,8 @@ FunctionGradient Numerical2PGradientCalculator::operator()(const MinimumParamete // std::cout << "position[" << i <<"] = " << par.Vec()(i) << std::endl << std::endl; // } + std::cout << "########### Numerical2PDerivative::operator()" < (1. - prec.Eps2())) { + if(yy < 0.) { + // Lower limit + // std::cout<<"SinParameterTransformation warning: is at its Lower allowed limit. "< (1. - prec.Eps2())) { - if(yy < 0.) { - // Lower limit - // std::cout<<"SinParameterTransformation warning: is at its Lower allowed limit. "< Differentiate(const double* x); - std::vector operator() (const double* x); - double GetFValue() const { - return fVal; - } - const double * GetG2() { - return &fG2[0]; - } - void SetStepTolerance(double value); - void SetGradTolerance(double value); - void SetNCycles(int value); - - void SetInitialValues(const double* g, const double* g2, const double* s); - - double Int2ext(const ROOT::Fit::ParameterSettings& parameter, double val) const; - double Ext2int(const ROOT::Fit::ParameterSettings& parameter, double val) const; - double DInt2Ext(const ROOT::Fit::ParameterSettings& parameter, double val) const; - - void SetInitialGradient(std::vector& parameters) const; - void SetParameterHasLimits(std::vector& parameters) const; - -private: - - // these are mutable because SetInitialGradient must be const because it's called - // from InitGradient which is const because DoDerivative must be const because the - // ROOT::Math::IMultiGradFunction interface requires this - mutable std::vector fGrd; - mutable std::vector fG2; - mutable std::vector fGstep; - // same story for SetParameterHasLimits - mutable std::vector _parameter_has_limits; - - const ROOT::Math::IBaseFunctionMultiDim* fFunction; - double fStepTolerance; - double fGradTolerance; - double fNCycles; - double fVal; - unsigned int fN; - double Up; - - // MODIFIED: Minuit2 determines machine precision itself in MnMachinePrecision.cxx, but - // mathcore isn't linked with minuit, so easier to pass in the correct eps from RooFit. - // This means precision is the caller's responsibility, beware! + class NumericalDerivatorMinuit2 { + public: + + NumericalDerivatorMinuit2(); + NumericalDerivatorMinuit2(const NumericalDerivatorMinuit2 &other); + NumericalDerivatorMinuit2& operator=(const NumericalDerivatorMinuit2 &other); + NumericalDerivatorMinuit2(const ROOT::Math::IBaseFunctionMultiDim &f, double step_tolerance, double grad_tolerance, unsigned int ncycles, double error_level);//, double precision); + // NumericalDerivatorMinuit2(const ROOT::Math::IBaseFunctionMultiDim &f, const ROOT::Fit::Fitter &fitter); + // NumericalDerivatorMinuit2(const ROOT::Math::IBaseFunctionMultiDim &f, const ROOT::Fit::Fitter &fitter, const ROOT::Minuit2::MnStrategy &strategy); + virtual ~NumericalDerivatorMinuit2(); + + ROOT::Minuit2::FunctionGradient Differentiate(const double* x, const std::vector& parameters); + ROOT::Minuit2::FunctionGradient operator() (const double* x, const std::vector& parameters); + double GetFValue() const { + return fVal; + } + const double * GetG2() { + return fG.G2().Data(); + } + void SetStepTolerance(double value); + void SetGradTolerance(double value); + void SetNCycles(int value); + + double Int2ext(const ROOT::Fit::ParameterSettings& parameter, double val) const; + double Ext2int(const ROOT::Fit::ParameterSettings& parameter, double val) const; + double DInt2Ext(const ROOT::Fit::ParameterSettings& parameter, double val) const; + double D2Int2Ext(const ROOT::Fit::ParameterSettings& parameter, double val) const; + double GStepInt2Ext(const ROOT::Fit::ParameterSettings& parameter, double val) const; + + void SetInitialGradient(std::vector& parameters) const; + void SetParameterHasLimits(std::vector& parameters) const; + + private: + + const ROOT::Math::IBaseFunctionMultiDim* fFunction; + + double fStepTolerance; + double fGradTolerance; + unsigned int fNCycles; + double Up; + double fVal; + unsigned int fN; + + // these are mutable because SetInitialGradient must be const because it's called + // from InitGradient which is const because DoDerivative must be const because the + // ROOT::Math::IMultiGradFunction interface requires this +// mutable std::vector fGrd; +// mutable std::vector fG2; +// mutable std::vector fGstep; + mutable ROOT::Minuit2::FunctionGradient fG; + // same story for SetParameterHasLimits + mutable std::vector _parameter_has_limits; + + + // MODIFIED: Minuit2 determines machine precision itself in MnMachinePrecision.cxx, but + // mathcore isn't linked with minuit, so easier to pass in the correct eps from RooFit. + // This means precision is the caller's responsibility, beware! // double eps; // double eps2; - // MODIFIED: Minuit2 determines machine precision in a slightly different way than - // std::numeric_limits::epsilon()). We go with the Minuit2 one. - ROOT::Minuit2::MnMachinePrecision precision; - - ROOT::Minuit2::SinParameterTransformation fDoubleLimTrafo; - ROOT::Minuit2::SqrtUpParameterTransformation fUpperLimTrafo; - ROOT::Minuit2::SqrtLowParameterTransformation fLowerLimTrafo; + // MODIFIED: Minuit2 determines machine precision in a slightly different way than + // std::numeric_limits::epsilon()). We go with the Minuit2 one. + ROOT::Minuit2::MnMachinePrecision precision; -}; + ROOT::Minuit2::SinParameterTransformation fDoubleLimTrafo; + ROOT::Minuit2::SqrtUpParameterTransformation fUpperLimTrafo; + ROOT::Minuit2::SqrtLowParameterTransformation fLowerLimTrafo; + }; } // namespace RooFit diff --git a/roofit/roofitcore/inc/RooGradMinimizerFcn.h b/roofit/roofitcore/inc/RooGradMinimizerFcn.h index c9afb1ea08f6a..3680207e6f9f9 100644 --- a/roofit/roofitcore/inc/RooGradMinimizerFcn.h +++ b/roofit/roofitcore/inc/RooGradMinimizerFcn.h @@ -31,65 +31,88 @@ #include #include +#include "Minuit2/FunctionGradient.h" + class RooGradMinimizer; class RooGradMinimizerFcn : public ROOT::Math::IMultiGradFunction { - public: +public: + + RooGradMinimizerFcn(RooAbsReal *funct, RooGradMinimizer *context, + bool verbose = false); + + RooGradMinimizerFcn(const RooGradMinimizerFcn &other); - RooGradMinimizerFcn(RooAbsReal *funct, RooGradMinimizer *context, - bool verbose = false); - RooGradMinimizerFcn(const RooGradMinimizerFcn& other); virtual ~RooGradMinimizerFcn(); - virtual ROOT::Math::IMultiGradFunction* Clone() const; - virtual unsigned int NDim() const { return _nDim; } + ROOT::Math::IMultiGradFunction *Clone() const override; + + unsigned int NDim() const override { return _nDim; } + + RooArgList *GetFloatParamList() { return _floatParamList; } + + RooArgList *GetConstParamList() { return _constParamList; } + + RooArgList *GetInitFloatParamList() { return _initFloatParamList; } - RooArgList* GetFloatParamList() { return _floatParamList; } - RooArgList* GetConstParamList() { return _constParamList; } - RooArgList* GetInitFloatParamList() { return _initFloatParamList; } - RooArgList* GetInitConstParamList() { return _initConstParamList; } + RooArgList *GetInitConstParamList() { return _initConstParamList; } - void SetEvalErrorWall(Bool_t flag) { _doEvalErrorWall = flag ; } - void SetPrintEvalErrors(Int_t numEvalErrors) { _printEvalErrors = numEvalErrors ; } - Bool_t SetLogFile(const char* inLogfile); - std::ofstream* GetLogFile() { return _logfile; } - void SetVerbose(Bool_t flag=kTRUE); + void SetEvalErrorWall(Bool_t flag) { _doEvalErrorWall = flag; } + + void SetPrintEvalErrors(Int_t numEvalErrors) { _printEvalErrors = numEvalErrors; } + + Bool_t SetLogFile(const char *inLogfile); + + std::ofstream *GetLogFile() { return _logfile; } + + void SetVerbose(Bool_t flag = kTRUE); + + Double_t &GetMaxFCN() { return _maxFCN; } - Double_t& GetMaxFCN() { return _maxFCN; } Int_t GetNumInvalidNLL() { return _numBadNLL; } - Bool_t Synchronize(std::vector& parameters, - Bool_t optConst, Bool_t verbose); - void BackProp(const ROOT::Fit::FitResult &results); - void ApplyCovarianceMatrix(TMatrixDSym& V); + Bool_t Synchronize(std::vector ¶meters, + Bool_t optConst, Bool_t verbose); + + void BackProp(const ROOT::Fit::FitResult &results); + + void ApplyCovarianceMatrix(TMatrixDSym &V); - Int_t evalCounter() const { return _evalCounter ; } - void zeroEvalCount() { _evalCounter = 0 ; } + Int_t evalCounter() const { return _evalCounter; } - void SynchronizeGradient(std::vector& parameters) const; + void zeroEvalCount() { _evalCounter = 0; } + + void SynchronizeGradient(std::vector ¶meters) const; + +private: - private: - Double_t GetPdfParamVal(Int_t index); + Double_t GetPdfParamErr(Int_t index); + void SetPdfParamErr(Int_t index, Double_t value); + void ClearPdfParamAsymErr(Int_t index); + void SetPdfParamErr(Int_t index, Double_t loVal, Double_t hiVal); inline Bool_t SetPdfParamVal(const Int_t &index, const Double_t &value) const; - virtual double DoEval(const double * x) const; - void updateFloatVec() ; + double DoEval(const double *x) const override; - virtual double DoDerivative(const double *x, unsigned int icoord) const; + void updateFloatVec(); + + double DoDerivative(const double *x, unsigned int icoord) const override; + + void run_derivator(const double *x) const; private: - mutable Int_t _evalCounter ; - + mutable Int_t _evalCounter; + RooAbsReal *_funct; RooGradMinimizer *_context; @@ -102,21 +125,27 @@ class RooGradMinimizerFcn : public ROOT::Math::IMultiGradFunction { std::ofstream *_logfile; bool _verbose; - RooArgList* _floatParamList; - std::vector _floatParamVec ; - RooArgList* _constParamList; - RooArgList* _initFloatParamList; - RooArgList* _initConstParamList; + RooArgList *_floatParamList; + std::vector _floatParamVec; + RooArgList *_constParamList; + RooArgList *_initFloatParamList; + RooArgList *_initConstParamList; // Before using any of the following members, call InitGradient! // this all needs to be mutable since ROOT::Math::IMultiGradFunction insists on DoDerivative being const mutable RooFit::NumericalDerivatorMinuit2 _gradf; - mutable std::vector _grad; +// mutable std::vector _grad; + mutable ROOT::Minuit2::FunctionGradient _grad; mutable std::vector _grad_params; mutable bool _grad_initialized; void InitGradient() const; -}; + bool hasG2ndDerivative() const override; + bool hasGStepSize() const override; + double DoSecondDerivative(const double *x, unsigned int icoord) const override; + double DoStepSize(const double *x, unsigned int icoord) const override; + +}; #endif #endif diff --git a/roofit/roofitcore/src/NumericalDerivatorMinuit2.cxx b/roofit/roofitcore/src/NumericalDerivatorMinuit2.cxx index 57e0d12b8d9ef..95f11ded1896b 100644 --- a/roofit/roofitcore/src/NumericalDerivatorMinuit2.cxx +++ b/roofit/roofitcore/src/NumericalDerivatorMinuit2.cxx @@ -1,10 +1,10 @@ // @(#)root/mathcore:$Id$ -// Authors: L. Moneta, J.T. Offermann 08/2013 -// E.G.P. Bos 09/2017 +// Authors: L. Moneta, J.T. Offermann, E.G.P. Bos 2013-2017 +// /********************************************************************** * * * Copyright (c) 2013 , LCG ROOT MathLib Team * - * Copyright (c) 2017 , Netherlands eScience Center * + * Copyright (c) 2017 Patrick Bos, Netherlands eScience Center * * * **********************************************************************/ /* @@ -37,84 +37,90 @@ namespace RooFit { -NumericalDerivatorMinuit2::NumericalDerivatorMinuit2() : - fFunction(0), - fStepTolerance(0.5), - fGradTolerance(0.1), - fNCycles(2), - fVal(0), - fN(0), - Up(1)//, + NumericalDerivatorMinuit2::NumericalDerivatorMinuit2() : + fFunction(0), + fStepTolerance(0.5), + fGradTolerance(0.1), + fNCycles(2), + Up(1), + fVal(0), + fN(0), + fG(0)//, // eps(std::numeric_limits::epsilon()), // eps2(2 * std::sqrt(eps)) -{} - - -NumericalDerivatorMinuit2::NumericalDerivatorMinuit2(const ROOT::Math::IBaseFunctionMultiDim &f, double step_tolerance, double grad_tolerance, unsigned int ncycles, double error_level)://, double precision): - fFunction(&f), - fStepTolerance(step_tolerance), - fGradTolerance(grad_tolerance), - fNCycles(ncycles), - Up(error_level)//, + {} + + + NumericalDerivatorMinuit2::NumericalDerivatorMinuit2(const ROOT::Math::IBaseFunctionMultiDim &f, double step_tolerance, double grad_tolerance, unsigned int ncycles, double error_level)://, double precision): + fFunction(&f), + fStepTolerance(step_tolerance), + fGradTolerance(grad_tolerance), + fNCycles(ncycles), + Up(error_level), + fVal(0), + fN(f.NDim()), + fG(f.NDim())//, // eps(precision), // eps2(2 * std::sqrt(eps)) -{ - // constructor with function, and tolerances (coordinates must be specified for differentiate function, not constructor) + { + // constructor with function, and tolerances (coordinates must be specified for differentiate function, not constructor) // fStepTolerance=step_tolerance; // fGradTolerance=grad_tolerance; // fFunction=&f; - - fN = f.NDim(); //number of dimensions, will look at vector size - fGrd.resize(fN); - fGstep.resize(fN); - fG2.resize(fN); + + ; //number of dimensions, will look at vector size +// fGrd.resize(fN); +// fGstep.resize(fN); +// fG2.resize(fN); _parameter_has_limits.resize(fN); - for (unsigned int i = 0; i NumericalDerivatorMinuit2::Differentiate(const double* cx) { + ROOT::Minuit2::FunctionGradient NumericalDerivatorMinuit2::Differentiate(const double* cx, const std::vector& parameters) { // std::cout << std::endl << "Start:" << std::endl; // for (unsigned int i = 0; i vx(fFunction->NDim()); assert (vx.size() > 0); + + TODO: x is anders in Numerical2PGradient (zie huidige test output), fix dat! + double *x = &vx[0]; std::copy (cx, cx+fFunction->NDim(), x); double step_tolerance = fStepTolerance; @@ -175,6 +179,10 @@ std::vector NumericalDerivatorMinuit2::Differentiate(const double* cx) { const ROOT::Math::IBaseFunctionMultiDim &f = *fFunction; fVal = f(x); //value of function at given points + ROOT::Minuit2::MnAlgebraicVector grad_vec(fG.Grad()), + gr2_vec(fG.G2()), + gstep_vec(fG.Gstep()); + // MODIFIED: Up // In Minuit2, this depends on the type of function to minimize, e.g. // chi-squared or negative log likelihood. It is set in the RooMinimizer @@ -185,158 +193,213 @@ std::vector NumericalDerivatorMinuit2::Differentiate(const double* cx) { double eps = precision.Eps(); double eps2 = precision.Eps2(); - // MODIFIED: two redundant double casts removed, for dfmin and for epspri + // MODIFIED: two redundant double casts removed, for dfmin and for epspri double dfmin = 8. * eps2 * (std::abs(fVal) + Up); double vrysml = 8.*eps*eps; unsigned int ncycle = fNCycles; for (int i = 0; i < int(fN); i++) { - double xtf = x[i]; - double epspri = eps2 + std::abs(fGrd[i] * eps2); - double step_old = 0.; - for (unsigned int j = 0; j < ncycle; ++ j) { - - double optstp = std::sqrt(dfmin/(std::abs(fG2[i])+epspri)); - double step = std::max(optstp, std::abs(0.1*fGstep[i])); - - // MODIFIED: in Minuit2 we have here the following condition: - // if(Trafo().Parameter(Trafo().ExtOfInt(i)).HasLimits()) { - // We replaced it by this: - if (_parameter_has_limits[i]) { - if(step > 0.5) step = 0.5; - } - // See the discussion above NumericalDerivatorMinuit2::SetInitialGradient - // below on how to pass parameter information to this derivator. - - double stpmax = 10.*std::abs(fGstep[i]); - if (step > stpmax) step = stpmax; - - double stpmin = std::max(vrysml, 8.*std::abs(eps2*x[i])); //8.*std::abs(double(eps2*x[i])) - if (step < stpmin) step = stpmin; - if (std::abs((step-step_old)/step) < step_tolerance) { - //answer = fGrd[i]; - break; - } - fGstep[i] = step; - step_old = step; - // std::cout << "step = " << step << std::endl; - x[i] = xtf + step; - //std::cout << "x[" << i << "] = " << x[i] <NDim()), + gr2_vec(fFunction->NDim()), + gstep_vec(fFunction->NDim()); + + unsigned ix = 0; + for (auto parameter = parameters.begin(); parameter != parameters.end(); ++parameter, ++ix) { // for (unsigned int i = 0; i < fN; ++i) { - // //double eps2 = TMath::Sqrt(fEpsilon); - // //double gsmin = 8.*eps2*(fabs(x[i])) + eps2; - // double dirin = s[i]; - // double g2 = 2.0 /(dirin*dirin); - // double gstep = 0.1*dirin; - // double grd = g2*dirin; - - // fGrd[i] = grd; - // fG2[i] = g2; - // fGstep[i] = gstep; - + // //double eps2 = TMath::Sqrt(fEpsilon); + // //double gsmin = 8.*eps2*(fabs(x[i])) + eps2; + // double dirin = s[i]; + // double g2 = 2.0 /(dirin*dirin); + // double gstep = 0.1*dirin; + // double grd = g2*dirin; + + // fGrd[i] = grd; + // fG2[i] = g2; + // fGstep[i] = gstep; + // unsigned int exOfIn = Trafo().ExtOfInt(i); // auto parameter = Trafo().Parameter(exOfIn); - // this should just be the parameter in the RooFit space ("external" in - // Minuit terms, since we're calculating the "external" gradient here) - // We get it from the loop. + // this should just be the parameter in the RooFit space ("external" in + // Minuit terms, since we're calculating the "external" gradient here) + // We get it from the loop. // double var = par.Vec()(i); - // I'm guessing par.Vec()(i) should give the value of variable i... + // I'm guessing par.Vec()(i) should give the value of variable i... // double var = parameter->Value(); - // Judging by the ParameterSettings.h constructor argument name "err", - // I'm guessing what MINUIT calls "Error" is stepsize on the ROOT side. + // Judging by the ParameterSettings.h constructor argument name "err", + // I'm guessing what MINUIT calls "Error" is stepsize on the ROOT side. // double werr = parameter->Error(); - double werr = parameter->StepSize(); + double werr = parameter->StepSize(); // double sav = Int2ext(*parameter, var); - // Actually, sav in Minuit2 is the external parameter value, so that is - // what we called var before and var is unnecessary here. - double sav = parameter->Value(); + // Actually, sav in Minuit2 is the external parameter value, so that is + // what we called var before and var is unnecessary here. + double sav = parameter->Value(); - // However, we do need var below, so let's calculate it using Ext2int: - double var = Ext2int(*parameter, sav); + // However, we do need var below, so let's calculate it using Ext2int: + double var = Ext2int(*parameter, sav); - // Int2Ext is not necessary, we're doing everything externally here - double sav2 = sav + werr; + // Int2Ext is not necessary, we're doing everything externally here + double sav2 = sav + werr; // double sav2 = var + werr; // if(parameter->HasLimits()) { // this if statement in MINUIT is superfluous - if(parameter->HasUpperLimit() && sav2 > parameter->UpperLimit()) { - sav2 = parameter->UpperLimit(); - } + if(parameter->HasUpperLimit() && sav2 > parameter->UpperLimit()) { + sav2 = parameter->UpperLimit(); + } - double var2 = Ext2int(*parameter, sav2); - // Ext2int is not necessary, we're doing everything externally here - double vplu = var2 - var; + double var2 = Ext2int(*parameter, sav2); + // Ext2int is not necessary, we're doing everything externally here + double vplu = var2 - var; // double vplu = sav2 - var; - sav2 = sav - werr; + sav2 = sav - werr; // sav2 = var - werr; // if(parameter->HasLimits()) { // this if statement in MINUIT is superfluous - if(parameter->HasLowerLimit() && sav2 < parameter->LowerLimit()) { - sav2 = parameter->LowerLimit(); - } + if(parameter->HasLowerLimit() && sav2 < parameter->LowerLimit()) { + sav2 = parameter->LowerLimit(); + } - var2 = Ext2int(*parameter, sav2); - // Ext2int is not necessary, we're doing everything externally here - double vmin = var2 - var; + var2 = Ext2int(*parameter, sav2); + // Ext2int is not necessary, we're doing everything externally here + double vmin = var2 - var; // double vmin = sav2 - var; - double gsmin = 8. * eps2 * (fabs(var) + eps2); - // protect against very small step sizes which can cause dirin to zero and then nan values in grd - double dirin = std::max(0.5*(fabs(vplu) + fabs(vmin)), gsmin ); + double gsmin = 8. * eps2 * (fabs(var) + eps2); + // protect against very small step sizes which can cause dirin to zero and then nan values in grd + double dirin = std::max(0.5*(fabs(vplu) + fabs(vmin)), gsmin ); // double g2 = 2.0*fFcn.ErrorDef()/(dirin*dirin); - // ErrorDef is the same as Up, which we already have in here - double g2 = 2.0*Up/(dirin*dirin); + // ErrorDef is the same as Up, which we already have in here + double g2 = 2.0*Up/(dirin*dirin); - double gstep = std::max(gsmin, 0.1*dirin); - double grd = g2*dirin; + double gstep = std::max(gsmin, 0.1*dirin); + double grd = g2*dirin; - if(parameter->IsBound()) { - if(gstep > 0.5) gstep = 0.5; - } + if(parameter->IsBound()) { + if(gstep > 0.5) gstep = 0.5; + } + + grad_vec(ix) = grd / DInt2Ext(*parameter, var); + gr2_vec(ix) = g2 / D2Int2Ext(*parameter, var); + gstep_vec(ix) = gstep / GStepInt2Ext(*parameter, var); - double inv_d_int_2_ext = 1 / DInt2Ext(*parameter, var); - fGrd[ix] = grd * inv_d_int_2_ext; - fG2[ix] = g2;// * inv_d_int_2_ext * inv_d_int_2_ext; - fGstep[ix] = gstep;// / inv_d_int_2_ext; + std::cout << "INTERNAL: fGrd[" << ix <<"] = " << grd << "\t"; + std::cout << "fG2[" << ix <<"] = " << g2 << "\t"; + std::cout << "fGstep[" << ix <<"] = " << gstep << std::endl; + std::cout << "EXTERNAL: fGrd[" << ix <<"] = " << grad_vec(ix) << "\t"; + std::cout << "fG2[" << ix <<"] = " << gr2_vec(ix) << "\t"; + std::cout << "fGstep[" << ix <<"] = " << gstep_vec(ix) << std::endl; - std::cout << "INTERNAL: fGrd[" << ix <<"] = " << grd << "\t"; - std::cout << "fG2[" << ix <<"] = " << g2 << "\t"; - std::cout << "fGstep[" << ix <<"] = " << gstep << std::endl; - std::cout << "EXTERNAL: fGrd[" << ix <<"] = " << fGrd[ix] << "\t"; - std::cout << "fG2[" << ix <<"] = " << fG2[ix] << "\t"; - std::cout << "fGstep[" << ix <<"] = " << fGstep[ix] << std::endl; + } + fG = ROOT::Minuit2::FunctionGradient(grad_vec, gr2_vec, gstep_vec); } -} } // namespace RooFit diff --git a/roofit/roofitcore/src/RooGradMinimizerFcn.cxx b/roofit/roofitcore/src/RooGradMinimizerFcn.cxx index 0eaf3b15f9dbc..c19a7fbf8ecf7 100644 --- a/roofit/roofitcore/src/RooGradMinimizerFcn.cxx +++ b/roofit/roofitcore/src/RooGradMinimizerFcn.cxx @@ -57,7 +57,8 @@ RooGradMinimizerFcn::RooGradMinimizerFcn(RooAbsReal *funct, RooGradMinimizer* co _maxFCN(-1e30), _numBadNLL(0), _printEvalErrors(10), _doEvalErrorWall(kTRUE), _nDim(0), _logfile(0), - _verbose(verbose), _grad_initialized(false) + _verbose(verbose), + _grad(0), _grad_initialized(false) { _evalCounter = 0 ; @@ -638,7 +639,7 @@ void RooGradMinimizerFcn::InitGradient() const { minimizer->ErrorDef());//, // precision.Eps()); _gradf = derivator; - _grad.resize(_nDim); +// _grad.resize(_nDim); _grad_params.resize(_nDim); SynchronizeGradient(fitter->Config().ParamsSettings()); @@ -647,31 +648,35 @@ void RooGradMinimizerFcn::InitGradient() const { } -double RooGradMinimizerFcn::DoDerivative(const double *x, unsigned int icoord) const { +void RooGradMinimizerFcn::run_derivator(const double *x) const { if (!_grad_initialized) { InitGradient(); } // check whether the derivative was already calculated for this set of parameters if (std::equal(_grad_params.begin(), _grad_params.end(), x)) { - std::cout << "grad value (cached) " << _grad[icoord] << std::endl; - return _grad[icoord]; - } - // if not, set the _grad_params to the current input parameters - std::vector new_grad_params(x, x + _nDim); - _grad_params = new_grad_params; + std::cout << "gradient already calculated for these parameters, use cached value" << std::endl; +// return _grad[icoord]; + } else { + // if not, set the _grad_params to the current input parameters + std::vector new_grad_params(x, x + _nDim); + _grad_params = new_grad_params; + + // Set the parameter values for this iteration + // TODO: this is already done in DoEval as well; find efficient way to do only once + for (int index = 0; index < _nDim; index++) { + if (_logfile) (*_logfile) << x[index] << " " ; + SetPdfParamVal(index,x[index]); + } - // Set the parameter values for this iteration - // TODO: this is already done in DoEval as well; find efficient way to do only once - for (int index = 0; index < _nDim; index++) { - if (_logfile) (*_logfile) << x[index] << " " ; - SetPdfParamVal(index,x[index]); + // Calculate the function for these parameters + _grad = _gradf(x, _context->fitter()->Config().ParamsSettings()); } +} - // Calculate the function for these parameters - _grad = _gradf(x); - +double RooGradMinimizerFcn::DoDerivative(const double *x, unsigned int icoord) const { + run_derivator(x); // std::cout << "grad value " << _grad[icoord] << std::endl; - return _grad[icoord]; + return _grad.Grad()(icoord); } @@ -680,5 +685,24 @@ void RooGradMinimizerFcn::SetVerbose(Bool_t flag) { } +bool RooGradMinimizerFcn::hasG2ndDerivative() const { + return true; +} + +bool RooGradMinimizerFcn::hasGStepSize() const { + return true; +} + +double RooGradMinimizerFcn::DoSecondDerivative(const double *x, unsigned int icoord) const { + run_derivator(x); + return _grad.G2()(icoord); +} + +double RooGradMinimizerFcn::DoStepSize(const double *x, unsigned int icoord) const { + run_derivator(x); + return _grad.Gstep()(icoord); +} + + #endif