From f10ed0768c6458d06a50219d6a6b203983143b20 Mon Sep 17 00:00:00 2001 From: "E. G. Patrick Bos" Date: Tue, 29 Jun 2021 14:09:57 +0200 Subject: [PATCH 1/7] add Minuit2 numerical derivator class to RooFit This class, based on original code by @lmoneta (https://github.com/lmoneta/root/blob/lvmini/math/mathcore/src/NumericalDerivator.cxx), replicates the Minuit2 numerical derivator calculations. Because it is a stand alone class, it can be used in RooFit to manually calculate gradients outside of Minuit2, e.g. in a parallelized way. --- roofit/roofitcore/CMakeLists.txt | 2 + .../inc/NumericalDerivatorMinuit2.h | 117 ++++++ .../src/NumericalDerivatorMinuit2.cxx | 381 ++++++++++++++++++ 3 files changed, 500 insertions(+) create mode 100644 roofit/roofitcore/inc/NumericalDerivatorMinuit2.h create mode 100644 roofit/roofitcore/src/NumericalDerivatorMinuit2.cxx diff --git a/roofit/roofitcore/CMakeLists.txt b/roofit/roofitcore/CMakeLists.txt index afe21fe236a6f..9e82d412d442a 100644 --- a/roofit/roofitcore/CMakeLists.txt +++ b/roofit/roofitcore/CMakeLists.txt @@ -12,6 +12,7 @@ ROOT_STANDARD_LIBRARY_PACKAGE(RooFitCore HEADERS Floats.h + NumericalDerivatorMinuit2.h Roo1DTable.h RooAbsAnaConvPdf.h RooAbsArg.h @@ -233,6 +234,7 @@ ROOT_STANDARD_LIBRARY_PACKAGE(RooFitCore SOURCES src/BidirMMapPipe.cxx src/BidirMMapPipe.h + src/NumericalDerivatorMinuit2.cxx src/Roo1DTable.cxx src/RooAbsAnaConvPdf.cxx src/RooAbsArg.cxx diff --git a/roofit/roofitcore/inc/NumericalDerivatorMinuit2.h b/roofit/roofitcore/inc/NumericalDerivatorMinuit2.h new file mode 100644 index 0000000000000..e23be6b7c2204 --- /dev/null +++ b/roofit/roofitcore/inc/NumericalDerivatorMinuit2.h @@ -0,0 +1,117 @@ +// @(#)root/mathcore:$Id$ +// Authors: L. Moneta, J.T. Offermann, E.G.P. Bos 2013-2018 +// +/********************************************************************** + * * + * Copyright (c) 2013 , LCG ROOT MathLib Team * + * * + **********************************************************************/ +/* + * NumericalDerivatorMinuit2.h + * + * Original version (NumericalDerivator) created on: Aug 14, 2013 + * Authors: L. Moneta, J. T. Offermann + * Modified version (NumericalDerivatorMinuit2) created on: Sep 27, 2017 + * Author: E. G. P. Bos + */ + +#ifndef RooFit_NumericalDerivatorMinuit2 +#define RooFit_NumericalDerivatorMinuit2 + +#ifndef ROOT_Math_IFunctionfwd +#include +#endif + +#include +#include "Fit/ParameterSettings.h" +#include "Minuit2/SinParameterTransformation.h" +#include "Minuit2/SqrtUpParameterTransformation.h" +#include "Minuit2/SqrtLowParameterTransformation.h" +#include "Minuit2/MnMachinePrecision.h" + +namespace RooFit { + +// Holds all necessary derivatives and associated numbers (per parameter) used in the NumericalDerivatorMinuit2 class. +struct MinuitDerivatorElement { + double derivative; + double second_derivative; + double step_size; +}; + +class NumericalDerivatorMinuit2 { +public: + explicit NumericalDerivatorMinuit2(bool always_exactly_mimic_minuit2 = true); + NumericalDerivatorMinuit2(const NumericalDerivatorMinuit2 &other); + NumericalDerivatorMinuit2(double step_tolerance, double grad_tolerance, unsigned int ncycles, double error_level, + bool always_exactly_mimic_minuit2 = true); + virtual ~NumericalDerivatorMinuit2(); + + void setup_differentiate(const ROOT::Math::IBaseFunctionMultiDim *function, const double *cx, + const std::vector ¶meters); + std::vector + Differentiate(const ROOT::Math::IBaseFunctionMultiDim *function, const double *x, + const std::vector ¶meters, + const std::vector& previous_gradient); + + MinuitDerivatorElement + partial_derivative(const ROOT::Math::IBaseFunctionMultiDim *function, const double *x, + const std::vector ¶meters, unsigned int i_component, MinuitDerivatorElement previous); + MinuitDerivatorElement fast_partial_derivative(const ROOT::Math::IBaseFunctionMultiDim *function, + const std::vector ¶meters, unsigned int i_component, const MinuitDerivatorElement& previous); + MinuitDerivatorElement operator()(const ROOT::Math::IBaseFunctionMultiDim *function, const double *x, + const std::vector ¶meters, + unsigned int i_component, const MinuitDerivatorElement& previous); + + double GetFValue() const { return fVal; } + void SetStepTolerance(double value); + void SetGradTolerance(double value); + void SetNCycles(int value); + + double Int2ext(const ROOT::Fit::ParameterSettings ¶meter, double val) const; + double Ext2int(const ROOT::Fit::ParameterSettings ¶meter, double val) const; + double DInt2Ext(const ROOT::Fit::ParameterSettings ¶meter, double val) const; + + void SetInitialGradient(const ROOT::Math::IBaseFunctionMultiDim *function, + const std::vector ¶meters, + std::vector& gradient); + + void set_step_tolerance(double step_tolerance); + void set_grad_tolerance(double grad_tolerance); + void set_ncycles(unsigned int ncycles); + void set_error_level(double error_level); + +private: + double fStepTolerance = 0.5; + double fGradTolerance = 0.1; + unsigned int fNCycles = 2; + double Up = 1; + double fVal = 0; + + std::vector vx, vx_external; + double dfmin; + double vrysml; + + // 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; + +private: + bool _always_exactly_mimic_minuit2; + +public: + bool always_exactly_mimic_minuit2() const; + void set_always_exactly_mimic_minuit2(bool flag); + +private: + std::vector vx_fVal_cache; +}; + +std::ostream& operator<<(std::ostream& out, const MinuitDerivatorElement value); + +} // namespace RooFit + +#endif /* NumericalDerivatorMinuit2_H_ */ \ No newline at end of file diff --git a/roofit/roofitcore/src/NumericalDerivatorMinuit2.cxx b/roofit/roofitcore/src/NumericalDerivatorMinuit2.cxx new file mode 100644 index 0000000000000..de9168c11d239 --- /dev/null +++ b/roofit/roofitcore/src/NumericalDerivatorMinuit2.cxx @@ -0,0 +1,381 @@ +// @(#)root/mathcore:$Id$ +// Authors: L. Moneta, J.T. Offermann, E.G.P. Bos 2013-2018 +// +/********************************************************************** + * * + * Copyright (c) 2013 , LCG ROOT MathLib Team * + * * + **********************************************************************/ +/* + * NumericalDerivatorMinuit2.cxx + * + * Original version (NumericalDerivator) created on: Aug 14, 2013 + * Authors: L. Moneta, J. T. Offermann + * Modified version (NumericalDerivatorMinuit2) created on: Sep 27, 2017 + * Author: E. G. P. Bos + * + * NumericalDerivator was essentially a slightly modified copy of code + * written by M. Winkler, F. James, L. Moneta, and A. Zsenei for Minuit2, + * Copyright (c) 2005 LCG ROOT Math team, CERN/PH-SFT. Original version: + * https://github.com/lmoneta/root/blob/lvmini/math/mathcore/src/NumericalDerivator.cxx + * + * This class attempts to more closely follow the Minuit2 implementation. + * Modified things (w.r.t. NumericalDerivator) are indicated by MODIFIED. + */ + +#include "NumericalDerivatorMinuit2.h" +#include +#include +#include +#include +#include +#include +#include "Fit/ParameterSettings.h" + +#include // needed here because in Fitter is only a forward declaration + +#include + +namespace RooFit { + +NumericalDerivatorMinuit2::NumericalDerivatorMinuit2(bool always_exactly_mimic_minuit2) + : _always_exactly_mimic_minuit2(always_exactly_mimic_minuit2) +{ +} + +NumericalDerivatorMinuit2::NumericalDerivatorMinuit2(double step_tolerance, double grad_tolerance, unsigned int ncycles, + double error_level, bool always_exactly_mimic_minuit2) + : fStepTolerance(step_tolerance), fGradTolerance(grad_tolerance), fNCycles(ncycles), Up(error_level), + _always_exactly_mimic_minuit2(always_exactly_mimic_minuit2) +{ +} + +// deep copy constructor +NumericalDerivatorMinuit2::NumericalDerivatorMinuit2(const RooFit::NumericalDerivatorMinuit2 &other) + : fStepTolerance(other.fStepTolerance), fGradTolerance(other.fGradTolerance), fNCycles(other.fNCycles), Up(other.Up), + fVal(other.fVal), vx(other.vx), vx_external(other.vx_external), dfmin(other.dfmin), vrysml(other.vrysml), + precision(other.precision), _always_exactly_mimic_minuit2(other._always_exactly_mimic_minuit2), + vx_fVal_cache(other.vx_fVal_cache) +{ +} + +void NumericalDerivatorMinuit2::SetStepTolerance(double value) +{ + fStepTolerance = value; +} + +void NumericalDerivatorMinuit2::SetGradTolerance(double value) +{ + fGradTolerance = value; +} + +void NumericalDerivatorMinuit2::SetNCycles(int value) +{ + fNCycles = value; +} + +NumericalDerivatorMinuit2::~NumericalDerivatorMinuit2() +{ + // TODO Auto-generated destructor stub +} + +// This function sets internal state based on input parameters. This state +// setup is used in the actual (partial) derivative calculations. +void NumericalDerivatorMinuit2::setup_differentiate(const ROOT::Math::IBaseFunctionMultiDim *function, const double *cx, + const std::vector ¶meters) +{ + + assert(function != nullptr && "function is a nullptr"); + + auto get_time = []() { + return std::chrono::duration_cast( + std::chrono::high_resolution_clock::now().time_since_epoch()) + .count(); + }; + decltype(get_time()) t1, t2, t3, t4, t5, t6, t7, t8; + + t1 = get_time(); + if (vx.size() != function->NDim()) { + vx.resize(function->NDim()); + } + t2 = get_time(); + if (vx_external.size() != function->NDim()) { + vx_external.resize(function->NDim()); + } + t3 = get_time(); + if (vx_fVal_cache.size() != function->NDim()) { + vx_fVal_cache.resize(function->NDim()); + } + t4 = get_time(); + std::copy(cx, cx + function->NDim(), vx.data()); + t5 = get_time(); + + // convert to Minuit external parameters + for (unsigned i = 0; i < function->NDim(); i++) { + vx_external[i] = Int2ext(parameters[i], vx[i]); + } + + t6 = get_time(); + + if (vx != vx_fVal_cache) { + vx_fVal_cache = vx; + fVal = (*function)(vx_external.data()); // value of function at given points + } + t7 = get_time(); + + dfmin = 8. * precision.Eps2() * (std::abs(fVal) + Up); + vrysml = 8. * precision.Eps() * precision.Eps(); + + t8 = get_time(); +} + +MinuitDerivatorElement +NumericalDerivatorMinuit2::partial_derivative(const ROOT::Math::IBaseFunctionMultiDim *function, const double *x, + const std::vector ¶meters, + unsigned int i_component, MinuitDerivatorElement previous) +{ + setup_differentiate(function, x, parameters); + return fast_partial_derivative(function, parameters, i_component, previous); +} + +// leaves the parameter setup to the caller +MinuitDerivatorElement NumericalDerivatorMinuit2::fast_partial_derivative(const ROOT::Math::IBaseFunctionMultiDim *function, + const std::vector ¶meters, + unsigned int ix, const MinuitDerivatorElement& previous) +{ + MinuitDerivatorElement deriv {previous}; + + double xtf = vx[ix]; + double epspri = precision.Eps2() + std::abs(deriv.derivative * precision.Eps2()); + double step_old = 0.; + + for (unsigned int j = 0; j < fNCycles; ++j) { + double optstp = std::sqrt(dfmin / (std::abs(deriv.second_derivative) + epspri)); + double step = std::max(optstp, std::abs(0.1 * deriv.step_size)); + + if (parameters[ix].IsBound()) { + if (step > 0.5) + step = 0.5; + } + + double stpmax = 10. * std::abs(deriv.step_size); + if (step > stpmax) + step = stpmax; + + double stpmin = std::max(vrysml, 8. * std::abs(precision.Eps2() * vx[ix])); + if (step < stpmin) + step = stpmin; + if (std::abs((step - step_old) / step) < fStepTolerance) { + break; + } + deriv.step_size = step; + step_old = step; + vx[ix] = xtf + step; + vx_external[ix] = Int2ext(parameters[ix], vx[ix]); + double fs1 = (*function)(vx_external.data()); + + vx[ix] = xtf - step; + vx_external[ix] = Int2ext(parameters[ix], vx[ix]); + double fs2 = (*function)(vx_external.data()); + + vx[ix] = xtf; + vx_external[ix] = Int2ext(parameters[ix], vx[ix]); + + double fGrd_old = deriv.derivative; + deriv.derivative = 0.5 * (fs1 - fs2) / step; + + deriv.second_derivative = (fs1 + fs2 - 2. * fVal) / step / step; + + if (std::abs(fGrd_old - deriv.derivative) / (std::abs(deriv.derivative) + dfmin / step) < fGradTolerance) { + break; + } + } + return deriv; +} + +MinuitDerivatorElement +NumericalDerivatorMinuit2::operator()(const ROOT::Math::IBaseFunctionMultiDim *function, const double *x, + const std::vector ¶meters, + unsigned int i_component, const MinuitDerivatorElement& previous) +{ + return partial_derivative(function, x, parameters, i_component, previous); +} + +std::vector +NumericalDerivatorMinuit2::Differentiate(const ROOT::Math::IBaseFunctionMultiDim *function, const double *cx, + const std::vector ¶meters, + const std::vector& previous_gradient) +{ + setup_differentiate(function, cx, parameters); + + std::vector gradient; + gradient.reserve(function->NDim()); + + for (unsigned int ix = 0; ix < function->NDim(); ++ix) { + gradient.emplace_back(fast_partial_derivative(function, parameters, ix, previous_gradient[ix])); + } + + return gradient; +} + + +double NumericalDerivatorMinuit2::Int2ext(const ROOT::Fit::ParameterSettings ¶meter, double val) const +{ + // return external value from internal value for parameter i + if (parameter.IsBound()) { + if (parameter.IsDoubleBound()) { + return fDoubleLimTrafo.Int2ext(val, parameter.UpperLimit(), parameter.LowerLimit()); + } else if (parameter.HasUpperLimit() && !parameter.HasLowerLimit()) { + return fUpperLimTrafo.Int2ext(val, parameter.UpperLimit()); + } else { + return fLowerLimTrafo.Int2ext(val, parameter.LowerLimit()); + } + } + + return val; +} + +double NumericalDerivatorMinuit2::Ext2int(const ROOT::Fit::ParameterSettings ¶meter, double val) const +{ + // return the internal value for parameter i with external value val + + if (parameter.IsBound()) { + if (parameter.IsDoubleBound()) { + return fDoubleLimTrafo.Ext2int(val, parameter.UpperLimit(), parameter.LowerLimit(), precision); + } else if (parameter.HasUpperLimit() && !parameter.HasLowerLimit()) { + return fUpperLimTrafo.Ext2int(val, parameter.UpperLimit(), precision); + } else { + return fLowerLimTrafo.Ext2int(val, parameter.LowerLimit(), precision); + } + } + + return val; +} + +double NumericalDerivatorMinuit2::DInt2Ext(const ROOT::Fit::ParameterSettings ¶meter, 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 (parameter.IsBound()) { + if (parameter.IsDoubleBound()) { + dd = fDoubleLimTrafo.DInt2Ext(val, parameter.UpperLimit(), parameter.LowerLimit()); + } else if (parameter.HasUpperLimit() && !parameter.HasLowerLimit()) { + dd = fUpperLimTrafo.DInt2Ext(val, parameter.UpperLimit()); + } else { + dd = fLowerLimTrafo.DInt2Ext(val, parameter.LowerLimit()); + } + } + + return dd; +} + +// MODIFIED: +// This function was not implemented as in Minuit2. Now it copies the behavior +// of InitialGradientCalculator. See https://github.com/roofit-dev/root/issues/10 +void NumericalDerivatorMinuit2::SetInitialGradient(const ROOT::Math::IBaseFunctionMultiDim *function, + const std::vector ¶meters, + std::vector& gradient) +{ + // set an initial gradient using some given steps + // (used in the first iteration) + auto get_time = []() { + return std::chrono::duration_cast( + std::chrono::high_resolution_clock::now().time_since_epoch()) + .count(); + }; + decltype(get_time()) t1, t2; + + t1 = get_time(); + + assert(function != nullptr && "function is a nullptr"); + + double eps2 = precision.Eps2(); + + unsigned ix = 0; + for (auto parameter = parameters.begin(); parameter != parameters.end(); ++parameter, ++ix) { + // What Minuit2 calls "Error" is stepsize on the ROOT side. + double werr = parameter->StepSize(); + + // 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); + + if (_always_exactly_mimic_minuit2) { + // this transformation can lose a few bits, but Minuit2 does it too + sav = Int2ext(*parameter, var); + } + + double sav2 = sav + werr; + + if (parameter->HasUpperLimit() && sav2 > parameter->UpperLimit()) { + sav2 = parameter->UpperLimit(); + } + + double var2 = Ext2int(*parameter, sav2); + double vplu = var2 - var; + + sav2 = sav - werr; + + if (parameter->HasLowerLimit() && sav2 < parameter->LowerLimit()) { + sav2 = parameter->LowerLimit(); + } + + var2 = Ext2int(*parameter, sav2); + double vmin = var2 - 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 g2 = 2.0 * Up / (dirin * dirin); + double gstep = std::max(gsmin, 0.1 * dirin); + double grd = g2 * dirin; + + if (parameter->IsBound()) { + if (gstep > 0.5) + gstep = 0.5; + } + + gradient[ix].derivative = grd; + gradient[ix].second_derivative = g2; + gradient[ix].step_size = gstep; + } + + t2 = get_time(); +} + +bool NumericalDerivatorMinuit2::always_exactly_mimic_minuit2() const +{ + return _always_exactly_mimic_minuit2; +}; + +void NumericalDerivatorMinuit2::set_always_exactly_mimic_minuit2(bool flag) +{ + _always_exactly_mimic_minuit2 = flag; +} + +void NumericalDerivatorMinuit2::set_step_tolerance(double step_tolerance) +{ + fStepTolerance = step_tolerance; +} +void NumericalDerivatorMinuit2::set_grad_tolerance(double grad_tolerance) +{ + fGradTolerance = grad_tolerance; +} +void NumericalDerivatorMinuit2::set_ncycles(unsigned int ncycles) +{ + fNCycles = ncycles; +} +void NumericalDerivatorMinuit2::set_error_level(double error_level) +{ + Up = error_level; +} + +std::ostream& operator<<(std::ostream& out, const MinuitDerivatorElement value){ + return out << "(derivative: " << value.derivative << ", second_derivative: " << value.second_derivative << ", step_size: " << value.step_size << ")"; +} + +} // namespace RooFit From cb46448214b91413a25ff10626e7cf068a48d172 Mon Sep 17 00:00:00 2001 From: "E. G. Patrick Bos" Date: Fri, 2 Jul 2021 08:59:11 +0200 Subject: [PATCH 2/7] Add Minuit2 cmake dependency --- roofit/roofitcore/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/roofit/roofitcore/CMakeLists.txt b/roofit/roofitcore/CMakeLists.txt index 9e82d412d442a..53c2e17fd0c30 100644 --- a/roofit/roofitcore/CMakeLists.txt +++ b/roofit/roofitcore/CMakeLists.txt @@ -456,6 +456,7 @@ ROOT_STANDARD_LIBRARY_PACKAGE(RooFitCore Matrix Tree Minuit + Minuit2 RIO MathCore Foam From e3e056d75bc737212ec44a5897d75464187dd071 Mon Sep 17 00:00:00 2001 From: "E. G. Patrick Bos" Date: Fri, 2 Jul 2021 09:39:21 +0200 Subject: [PATCH 3/7] Remove timing code --- .../src/NumericalDerivatorMinuit2.cxx | 28 ------------------- 1 file changed, 28 deletions(-) diff --git a/roofit/roofitcore/src/NumericalDerivatorMinuit2.cxx b/roofit/roofitcore/src/NumericalDerivatorMinuit2.cxx index de9168c11d239..da531ee698e4e 100644 --- a/roofit/roofitcore/src/NumericalDerivatorMinuit2.cxx +++ b/roofit/roofitcore/src/NumericalDerivatorMinuit2.cxx @@ -84,49 +84,31 @@ NumericalDerivatorMinuit2::~NumericalDerivatorMinuit2() void NumericalDerivatorMinuit2::setup_differentiate(const ROOT::Math::IBaseFunctionMultiDim *function, const double *cx, const std::vector ¶meters) { - assert(function != nullptr && "function is a nullptr"); - auto get_time = []() { - return std::chrono::duration_cast( - std::chrono::high_resolution_clock::now().time_since_epoch()) - .count(); - }; - decltype(get_time()) t1, t2, t3, t4, t5, t6, t7, t8; - - t1 = get_time(); if (vx.size() != function->NDim()) { vx.resize(function->NDim()); } - t2 = get_time(); if (vx_external.size() != function->NDim()) { vx_external.resize(function->NDim()); } - t3 = get_time(); if (vx_fVal_cache.size() != function->NDim()) { vx_fVal_cache.resize(function->NDim()); } - t4 = get_time(); std::copy(cx, cx + function->NDim(), vx.data()); - t5 = get_time(); // convert to Minuit external parameters for (unsigned i = 0; i < function->NDim(); i++) { vx_external[i] = Int2ext(parameters[i], vx[i]); } - t6 = get_time(); - if (vx != vx_fVal_cache) { vx_fVal_cache = vx; fVal = (*function)(vx_external.data()); // value of function at given points } - t7 = get_time(); dfmin = 8. * precision.Eps2() * (std::abs(fVal) + Up); vrysml = 8. * precision.Eps() * precision.Eps(); - - t8 = get_time(); } MinuitDerivatorElement @@ -280,14 +262,6 @@ void NumericalDerivatorMinuit2::SetInitialGradient(const ROOT::Math::IBaseFuncti { // set an initial gradient using some given steps // (used in the first iteration) - auto get_time = []() { - return std::chrono::duration_cast( - std::chrono::high_resolution_clock::now().time_since_epoch()) - .count(); - }; - decltype(get_time()) t1, t2; - - t1 = get_time(); assert(function != nullptr && "function is a nullptr"); @@ -343,8 +317,6 @@ void NumericalDerivatorMinuit2::SetInitialGradient(const ROOT::Math::IBaseFuncti gradient[ix].second_derivative = g2; gradient[ix].step_size = gstep; } - - t2 = get_time(); } bool NumericalDerivatorMinuit2::always_exactly_mimic_minuit2() const From a530b6305bf161e6671c9d4902673251d1e8c43b Mon Sep 17 00:00:00 2001 From: "E. G. Patrick Bos" Date: Wed, 7 Jul 2021 22:34:23 +0200 Subject: [PATCH 4/7] move NumericalDerivator from roofitcore to Minuit2 --- math/minuit2/CMakeLists.txt | 2 + math/minuit2/inc/Minuit2/NumericalDerivator.h | 114 ++++++ math/minuit2/src/NumericalDerivator.cxx | 336 +++++++++++++++++ roofit/roofitcore/CMakeLists.txt | 2 - .../inc/NumericalDerivatorMinuit2.h | 117 ------ .../src/NumericalDerivatorMinuit2.cxx | 353 ------------------ 6 files changed, 452 insertions(+), 472 deletions(-) create mode 100644 math/minuit2/inc/Minuit2/NumericalDerivator.h create mode 100644 math/minuit2/src/NumericalDerivator.cxx delete mode 100644 roofit/roofitcore/inc/NumericalDerivatorMinuit2.h delete mode 100644 roofit/roofitcore/src/NumericalDerivatorMinuit2.cxx diff --git a/math/minuit2/CMakeLists.txt b/math/minuit2/CMakeLists.txt index fac5ca399789e..6badb5fa4fd01 100644 --- a/math/minuit2/CMakeLists.txt +++ b/math/minuit2/CMakeLists.txt @@ -106,6 +106,7 @@ if(CMAKE_PROJECT_NAME STREQUAL ROOT) Minuit2/ModularFunctionMinimizer.h Minuit2/NegativeG2LineSearch.h Minuit2/Numerical2PGradientCalculator.h + Minuit2/NumericalDerivator.h Minuit2/ParametricFunction.h Minuit2/ScanBuilder.h Minuit2/ScanMinimizer.h @@ -176,6 +177,7 @@ if(CMAKE_PROJECT_NAME STREQUAL ROOT) src/ModularFunctionMinimizer.cxx src/NegativeG2LineSearch.cxx src/Numerical2PGradientCalculator.cxx + src/NumericalDerivator.cxx src/ParametricFunction.cxx src/ScanBuilder.cxx src/SimplexBuilder.cxx diff --git a/math/minuit2/inc/Minuit2/NumericalDerivator.h b/math/minuit2/inc/Minuit2/NumericalDerivator.h new file mode 100644 index 0000000000000..a7343d9bd9f10 --- /dev/null +++ b/math/minuit2/inc/Minuit2/NumericalDerivator.h @@ -0,0 +1,114 @@ +// @(#)root/mathcore:$Id$ +// Authors: L. Moneta, J.T. Offermann, E.G.P. Bos 2013-2018 +// +/********************************************************************** + * * + * Copyright (c) 2013 , LCG ROOT MathLib Team * + * * + **********************************************************************/ +/* + * NumericalDerivator.h + * + * Original version created on: Aug 14, 2013 + * Authors: L. Moneta, J. T. Offermann + * Modified version created on: Sep 27, 2017 + * Author: E. G. P. Bos + */ + +#ifndef ROOT_Minuit2_NumericalDerivator +#define ROOT_Minuit2_NumericalDerivator + +#ifndef ROOT_Math_IFunctionfwd +#include +#endif + +#include +#include "Fit/ParameterSettings.h" +#include "Minuit2/SinParameterTransformation.h" +#include "Minuit2/SqrtUpParameterTransformation.h" +#include "Minuit2/SqrtLowParameterTransformation.h" +#include "Minuit2/MnMachinePrecision.h" + +namespace ROOT { +namespace Minuit2 { + +// Holds all necessary derivatives and associated numbers (per parameter) used in the NumericalDerivator class. +struct DerivatorElement { + double derivative; + double second_derivative; + double step_size; +}; + +class NumericalDerivator { +public: + explicit NumericalDerivator(bool always_exactly_mimic_minuit2 = true); + NumericalDerivator(const NumericalDerivator &other); + NumericalDerivator(double step_tolerance, double grad_tolerance, unsigned int ncycles, double error_level, + bool always_exactly_mimic_minuit2 = true); + + void SetupDifferentiate(const ROOT::Math::IBaseFunctionMultiDim *function, const double *cx, + const std::vector ¶meters); + std::vector Differentiate(const ROOT::Math::IBaseFunctionMultiDim *function, const double *x, + const std::vector ¶meters, + const std::vector &previous_gradient); + + DerivatorElement PartialDerivative(const ROOT::Math::IBaseFunctionMultiDim *function, const double *x, + const std::vector ¶meters, + unsigned int i_component, DerivatorElement previous); + DerivatorElement FastPartialDerivative(const ROOT::Math::IBaseFunctionMultiDim *function, + const std::vector ¶meters, + unsigned int i_component, const DerivatorElement &previous); + DerivatorElement operator()(const ROOT::Math::IBaseFunctionMultiDim *function, const double *x, + const std::vector ¶meters, unsigned int i_component, + const DerivatorElement &previous); + + double GetValue() const { return fVal; } + void SetStepTolerance(double value); + void SetGradTolerance(double value); + void SetNCycles(unsigned int value); + void SetErrorLevel(double value); + + double Int2ext(const ROOT::Fit::ParameterSettings ¶meter, double val) const; + double Ext2int(const ROOT::Fit::ParameterSettings ¶meter, double val) const; + double DInt2Ext(const ROOT::Fit::ParameterSettings ¶meter, double val) const; + + void SetInitialGradient(const ROOT::Math::IBaseFunctionMultiDim *function, + const std::vector ¶meters, + std::vector &gradient); + +private: + double fStepTolerance = 0.5; + double fGradTolerance = 0.1; + unsigned int fNCycles = 2; + double fUp = 1; + double fVal = 0; + + std::vector fVx, fVxExternal; + double fDfmin; + double fVrysml; + + // MODIFIED: Minuit2 determines machine precision in a slightly different way than + // std::numeric_limits::epsilon()). We go with the Minuit2 one. + ROOT::Minuit2::MnMachinePrecision fPrecision; + + ROOT::Minuit2::SinParameterTransformation fDoubleLimTrafo; + ROOT::Minuit2::SqrtUpParameterTransformation fUpperLimTrafo; + ROOT::Minuit2::SqrtLowParameterTransformation fLowerLimTrafo; + +private: + bool fAlwaysExactlyMimicMinuit2; + +public: + bool AlwaysExactlyMimicMinuit2() const; + void SetAlwaysExactlyMimicMinuit2(bool flag); + +private: + std::vector fVxFValCache; +}; + +std::ostream &operator<<(std::ostream &out, const DerivatorElement &value); + +} // namespace Minuit2 +} // namespace ROOT + +#endif // ROOT_Minuit2_NumericalDerivator \ No newline at end of file diff --git a/math/minuit2/src/NumericalDerivator.cxx b/math/minuit2/src/NumericalDerivator.cxx new file mode 100644 index 0000000000000..e099560ae7a6d --- /dev/null +++ b/math/minuit2/src/NumericalDerivator.cxx @@ -0,0 +1,336 @@ +// @(#)root/mathcore:$Id$ +// Authors: L. Moneta, J.T. Offermann, E.G.P. Bos 2013-2018 +// +/********************************************************************** + * * + * Copyright (c) 2013 , LCG ROOT MathLib Team * + * * + **********************************************************************/ +/* + * NumericalDerivator.cxx + * + * Original version created on: Aug 14, 2013 + * Authors: L. Moneta, J. T. Offermann + * Modified version created on: Sep 27, 2017 + * Author: E. G. P. Bos + * + * NumericalDerivator was essentially a slightly modified copy of code + * written by M. Winkler, F. James, L. Moneta, and A. Zsenei for Minuit2, + * Copyright (c) 2005 LCG ROOT Math team, CERN/PH-SFT. Original version: + * https://github.com/lmoneta/root/blob/lvmini/math/mathcore/src/NumericalDerivator.cxx + * + * This class attempts to more closely follow the Minuit2 implementation. + * Modified things (w.r.t. original) are indicated by MODIFIED. + */ + +#include "Minuit2/NumericalDerivator.h" +#include +#include +#include +#include +#include +#include +#include "Fit/ParameterSettings.h" + +#include // needed here because in Fitter is only a forward declaration + +namespace ROOT { +namespace Minuit2 { + +NumericalDerivator::NumericalDerivator(bool always_exactly_mimic_minuit2) + : fAlwaysExactlyMimicMinuit2(always_exactly_mimic_minuit2) +{ +} + +NumericalDerivator::NumericalDerivator(double step_tolerance, double grad_tolerance, unsigned int ncycles, + double error_level, bool always_exactly_mimic_minuit2) + : fStepTolerance(step_tolerance), fGradTolerance(grad_tolerance), fNCycles(ncycles), fUp(error_level), + fAlwaysExactlyMimicMinuit2(always_exactly_mimic_minuit2) +{ +} + +// deep copy constructor +NumericalDerivator::NumericalDerivator(const NumericalDerivator &other) + : fStepTolerance(other.fStepTolerance), fGradTolerance(other.fGradTolerance), fNCycles(other.fNCycles), + fUp(other.fUp), fVal(other.fVal), fVx(other.fVx), fVxExternal(other.fVxExternal), fDfmin(other.fDfmin), + fVrysml(other.fVrysml), fPrecision(other.fPrecision), fAlwaysExactlyMimicMinuit2(other.fAlwaysExactlyMimicMinuit2), + fVxFValCache(other.fVxFValCache) +{ +} + +void NumericalDerivator::SetStepTolerance(double value) +{ + fStepTolerance = value; +} + +void NumericalDerivator::SetGradTolerance(double value) +{ + fGradTolerance = value; +} + +void NumericalDerivator::SetNCycles(unsigned int value) +{ + fNCycles = value; +} + +void NumericalDerivator::SetErrorLevel(double value) +{ + fUp = value; +} + +// This function sets internal state based on input parameters. This state +// setup is used in the actual (partial) derivative calculations. +void NumericalDerivator::SetupDifferentiate(const ROOT::Math::IBaseFunctionMultiDim *function, const double *cx, + const std::vector ¶meters) +{ + assert(function != nullptr && "function is a nullptr"); + + if (fVx.size() != function->NDim()) { + fVx.resize(function->NDim()); + } + if (fVxExternal.size() != function->NDim()) { + fVxExternal.resize(function->NDim()); + } + if (fVxFValCache.size() != function->NDim()) { + fVxFValCache.resize(function->NDim()); + } + std::copy(cx, cx + function->NDim(), fVx.data()); + + // convert to Minuit external parameters + for (unsigned i = 0; i < function->NDim(); i++) { + fVxExternal[i] = Int2ext(parameters[i], fVx[i]); + } + + if (fVx != fVxFValCache) { + fVxFValCache = fVx; + fVal = (*function)(fVxExternal.data()); // value of function at given points + } + + fDfmin = 8. * fPrecision.Eps2() * (std::abs(fVal) + fUp); + fVrysml = 8. * fPrecision.Eps() * fPrecision.Eps(); +} + +DerivatorElement NumericalDerivator::PartialDerivative(const ROOT::Math::IBaseFunctionMultiDim *function, + const double *x, + const std::vector ¶meters, + unsigned int i_component, DerivatorElement previous) +{ + SetupDifferentiate(function, x, parameters); + return FastPartialDerivative(function, parameters, i_component, previous); +} + +// leaves the parameter setup to the caller +DerivatorElement NumericalDerivator::FastPartialDerivative(const ROOT::Math::IBaseFunctionMultiDim *function, + const std::vector ¶meters, + unsigned int i_component, const DerivatorElement &previous) +{ + DerivatorElement deriv{previous}; + + double xtf = fVx[i_component]; + double epspri = fPrecision.Eps2() + std::abs(deriv.derivative * fPrecision.Eps2()); + double step_old = 0.; + + for (unsigned int j = 0; j < fNCycles; ++j) { + double optstp = std::sqrt(fDfmin / (std::abs(deriv.second_derivative) + epspri)); + double step = std::max(optstp, std::abs(0.1 * deriv.step_size)); + + if (parameters[i_component].IsBound()) { + if (step > 0.5) + step = 0.5; + } + + double stpmax = 10. * std::abs(deriv.step_size); + if (step > stpmax) + step = stpmax; + + double stpmin = std::max(fVrysml, 8. * std::abs(fPrecision.Eps2() * fVx[i_component])); + if (step < stpmin) + step = stpmin; + if (std::abs((step - step_old) / step) < fStepTolerance) { + break; + } + deriv.step_size = step; + step_old = step; + fVx[i_component] = xtf + step; + fVxExternal[i_component] = Int2ext(parameters[i_component], fVx[i_component]); + double fs1 = (*function)(fVxExternal.data()); + + fVx[i_component] = xtf - step; + fVxExternal[i_component] = Int2ext(parameters[i_component], fVx[i_component]); + double fs2 = (*function)(fVxExternal.data()); + + fVx[i_component] = xtf; + fVxExternal[i_component] = Int2ext(parameters[i_component], fVx[i_component]); + + double fGrd_old = deriv.derivative; + deriv.derivative = 0.5 * (fs1 - fs2) / step; + + deriv.second_derivative = (fs1 + fs2 - 2. * fVal) / step / step; + + if (std::abs(fGrd_old - deriv.derivative) / (std::abs(deriv.derivative) + fDfmin / step) < fGradTolerance) { + break; + } + } + return deriv; +} + +DerivatorElement NumericalDerivator::operator()(const ROOT::Math::IBaseFunctionMultiDim *function, const double *x, + const std::vector ¶meters, + unsigned int i_component, const DerivatorElement &previous) +{ + return PartialDerivative(function, x, parameters, i_component, previous); +} + +std::vector +NumericalDerivator::Differentiate(const ROOT::Math::IBaseFunctionMultiDim *function, const double *cx, + const std::vector ¶meters, + const std::vector &previous_gradient) +{ + SetupDifferentiate(function, cx, parameters); + + std::vector gradient; + gradient.reserve(function->NDim()); + + for (unsigned int ix = 0; ix < function->NDim(); ++ix) { + gradient.emplace_back(FastPartialDerivative(function, parameters, ix, previous_gradient[ix])); + } + + return gradient; +} + +double NumericalDerivator::Int2ext(const ROOT::Fit::ParameterSettings ¶meter, double val) const +{ + // return external value from internal value for parameter i + if (parameter.IsBound()) { + if (parameter.IsDoubleBound()) { + return fDoubleLimTrafo.Int2ext(val, parameter.UpperLimit(), parameter.LowerLimit()); + } else if (parameter.HasUpperLimit() && !parameter.HasLowerLimit()) { + return fUpperLimTrafo.Int2ext(val, parameter.UpperLimit()); + } else { + return fLowerLimTrafo.Int2ext(val, parameter.LowerLimit()); + } + } + + return val; +} + +double NumericalDerivator::Ext2int(const ROOT::Fit::ParameterSettings ¶meter, double val) const +{ + // return the internal value for parameter i with external value val + + if (parameter.IsBound()) { + if (parameter.IsDoubleBound()) { + return fDoubleLimTrafo.Ext2int(val, parameter.UpperLimit(), parameter.LowerLimit(), fPrecision); + } else if (parameter.HasUpperLimit() && !parameter.HasLowerLimit()) { + return fUpperLimTrafo.Ext2int(val, parameter.UpperLimit(), fPrecision); + } else { + return fLowerLimTrafo.Ext2int(val, parameter.LowerLimit(), fPrecision); + } + } + + return val; +} + +double NumericalDerivator::DInt2Ext(const ROOT::Fit::ParameterSettings ¶meter, 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 (parameter.IsBound()) { + if (parameter.IsDoubleBound()) { + dd = fDoubleLimTrafo.DInt2Ext(val, parameter.UpperLimit(), parameter.LowerLimit()); + } else if (parameter.HasUpperLimit() && !parameter.HasLowerLimit()) { + dd = fUpperLimTrafo.DInt2Ext(val, parameter.UpperLimit()); + } else { + dd = fLowerLimTrafo.DInt2Ext(val, parameter.LowerLimit()); + } + } + + return dd; +} + +// MODIFIED: +// This function was not implemented as in Minuit2. Now it copies the behavior +// of InitialGradientCalculator. See https://github.com/roofit-dev/root/issues/10 +void NumericalDerivator::SetInitialGradient(const ROOT::Math::IBaseFunctionMultiDim *function, + const std::vector ¶meters, + std::vector &gradient) +{ + // set an initial gradient using some given steps + // (used in the first iteration) + + assert(function != nullptr && "function is a nullptr"); + + double eps2 = fPrecision.Eps2(); + + unsigned ix = 0; + for (auto parameter = parameters.begin(); parameter != parameters.end(); ++parameter, ++ix) { + // What Minuit2 calls "Error" is stepsize on the ROOT side. + double werr = parameter->StepSize(); + + // 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); + + if (fAlwaysExactlyMimicMinuit2) { + // this transformation can lose a few bits, but Minuit2 does it too + sav = Int2ext(*parameter, var); + } + + double sav2 = sav + werr; + + if (parameter->HasUpperLimit() && sav2 > parameter->UpperLimit()) { + sav2 = parameter->UpperLimit(); + } + + double var2 = Ext2int(*parameter, sav2); + double vplu = var2 - var; + + sav2 = sav - werr; + + if (parameter->HasLowerLimit() && sav2 < parameter->LowerLimit()) { + sav2 = parameter->LowerLimit(); + } + + var2 = Ext2int(*parameter, sav2); + double vmin = var2 - 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 g2 = 2.0 * fUp / (dirin * dirin); + double gstep = std::max(gsmin, 0.1 * dirin); + double grd = g2 * dirin; + + if (parameter->IsBound()) { + if (gstep > 0.5) + gstep = 0.5; + } + + gradient[ix].derivative = grd; + gradient[ix].second_derivative = g2; + gradient[ix].step_size = gstep; + } +} + +bool NumericalDerivator::AlwaysExactlyMimicMinuit2() const +{ + return fAlwaysExactlyMimicMinuit2; +}; + +void NumericalDerivator::SetAlwaysExactlyMimicMinuit2(bool flag) +{ + fAlwaysExactlyMimicMinuit2 = flag; +} + +std::ostream &operator<<(std::ostream &out, const DerivatorElement &value) +{ + return out << "(derivative: " << value.derivative << ", second_derivative: " << value.second_derivative + << ", step_size: " << value.step_size << ")"; +} + +} // namespace Minuit2 +} // namespace ROOT diff --git a/roofit/roofitcore/CMakeLists.txt b/roofit/roofitcore/CMakeLists.txt index 53c2e17fd0c30..371fad77eccea 100644 --- a/roofit/roofitcore/CMakeLists.txt +++ b/roofit/roofitcore/CMakeLists.txt @@ -12,7 +12,6 @@ ROOT_STANDARD_LIBRARY_PACKAGE(RooFitCore HEADERS Floats.h - NumericalDerivatorMinuit2.h Roo1DTable.h RooAbsAnaConvPdf.h RooAbsArg.h @@ -234,7 +233,6 @@ ROOT_STANDARD_LIBRARY_PACKAGE(RooFitCore SOURCES src/BidirMMapPipe.cxx src/BidirMMapPipe.h - src/NumericalDerivatorMinuit2.cxx src/Roo1DTable.cxx src/RooAbsAnaConvPdf.cxx src/RooAbsArg.cxx diff --git a/roofit/roofitcore/inc/NumericalDerivatorMinuit2.h b/roofit/roofitcore/inc/NumericalDerivatorMinuit2.h deleted file mode 100644 index e23be6b7c2204..0000000000000 --- a/roofit/roofitcore/inc/NumericalDerivatorMinuit2.h +++ /dev/null @@ -1,117 +0,0 @@ -// @(#)root/mathcore:$Id$ -// Authors: L. Moneta, J.T. Offermann, E.G.P. Bos 2013-2018 -// -/********************************************************************** - * * - * Copyright (c) 2013 , LCG ROOT MathLib Team * - * * - **********************************************************************/ -/* - * NumericalDerivatorMinuit2.h - * - * Original version (NumericalDerivator) created on: Aug 14, 2013 - * Authors: L. Moneta, J. T. Offermann - * Modified version (NumericalDerivatorMinuit2) created on: Sep 27, 2017 - * Author: E. G. P. Bos - */ - -#ifndef RooFit_NumericalDerivatorMinuit2 -#define RooFit_NumericalDerivatorMinuit2 - -#ifndef ROOT_Math_IFunctionfwd -#include -#endif - -#include -#include "Fit/ParameterSettings.h" -#include "Minuit2/SinParameterTransformation.h" -#include "Minuit2/SqrtUpParameterTransformation.h" -#include "Minuit2/SqrtLowParameterTransformation.h" -#include "Minuit2/MnMachinePrecision.h" - -namespace RooFit { - -// Holds all necessary derivatives and associated numbers (per parameter) used in the NumericalDerivatorMinuit2 class. -struct MinuitDerivatorElement { - double derivative; - double second_derivative; - double step_size; -}; - -class NumericalDerivatorMinuit2 { -public: - explicit NumericalDerivatorMinuit2(bool always_exactly_mimic_minuit2 = true); - NumericalDerivatorMinuit2(const NumericalDerivatorMinuit2 &other); - NumericalDerivatorMinuit2(double step_tolerance, double grad_tolerance, unsigned int ncycles, double error_level, - bool always_exactly_mimic_minuit2 = true); - virtual ~NumericalDerivatorMinuit2(); - - void setup_differentiate(const ROOT::Math::IBaseFunctionMultiDim *function, const double *cx, - const std::vector ¶meters); - std::vector - Differentiate(const ROOT::Math::IBaseFunctionMultiDim *function, const double *x, - const std::vector ¶meters, - const std::vector& previous_gradient); - - MinuitDerivatorElement - partial_derivative(const ROOT::Math::IBaseFunctionMultiDim *function, const double *x, - const std::vector ¶meters, unsigned int i_component, MinuitDerivatorElement previous); - MinuitDerivatorElement fast_partial_derivative(const ROOT::Math::IBaseFunctionMultiDim *function, - const std::vector ¶meters, unsigned int i_component, const MinuitDerivatorElement& previous); - MinuitDerivatorElement operator()(const ROOT::Math::IBaseFunctionMultiDim *function, const double *x, - const std::vector ¶meters, - unsigned int i_component, const MinuitDerivatorElement& previous); - - double GetFValue() const { return fVal; } - void SetStepTolerance(double value); - void SetGradTolerance(double value); - void SetNCycles(int value); - - double Int2ext(const ROOT::Fit::ParameterSettings ¶meter, double val) const; - double Ext2int(const ROOT::Fit::ParameterSettings ¶meter, double val) const; - double DInt2Ext(const ROOT::Fit::ParameterSettings ¶meter, double val) const; - - void SetInitialGradient(const ROOT::Math::IBaseFunctionMultiDim *function, - const std::vector ¶meters, - std::vector& gradient); - - void set_step_tolerance(double step_tolerance); - void set_grad_tolerance(double grad_tolerance); - void set_ncycles(unsigned int ncycles); - void set_error_level(double error_level); - -private: - double fStepTolerance = 0.5; - double fGradTolerance = 0.1; - unsigned int fNCycles = 2; - double Up = 1; - double fVal = 0; - - std::vector vx, vx_external; - double dfmin; - double vrysml; - - // 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; - -private: - bool _always_exactly_mimic_minuit2; - -public: - bool always_exactly_mimic_minuit2() const; - void set_always_exactly_mimic_minuit2(bool flag); - -private: - std::vector vx_fVal_cache; -}; - -std::ostream& operator<<(std::ostream& out, const MinuitDerivatorElement value); - -} // namespace RooFit - -#endif /* NumericalDerivatorMinuit2_H_ */ \ No newline at end of file diff --git a/roofit/roofitcore/src/NumericalDerivatorMinuit2.cxx b/roofit/roofitcore/src/NumericalDerivatorMinuit2.cxx deleted file mode 100644 index da531ee698e4e..0000000000000 --- a/roofit/roofitcore/src/NumericalDerivatorMinuit2.cxx +++ /dev/null @@ -1,353 +0,0 @@ -// @(#)root/mathcore:$Id$ -// Authors: L. Moneta, J.T. Offermann, E.G.P. Bos 2013-2018 -// -/********************************************************************** - * * - * Copyright (c) 2013 , LCG ROOT MathLib Team * - * * - **********************************************************************/ -/* - * NumericalDerivatorMinuit2.cxx - * - * Original version (NumericalDerivator) created on: Aug 14, 2013 - * Authors: L. Moneta, J. T. Offermann - * Modified version (NumericalDerivatorMinuit2) created on: Sep 27, 2017 - * Author: E. G. P. Bos - * - * NumericalDerivator was essentially a slightly modified copy of code - * written by M. Winkler, F. James, L. Moneta, and A. Zsenei for Minuit2, - * Copyright (c) 2005 LCG ROOT Math team, CERN/PH-SFT. Original version: - * https://github.com/lmoneta/root/blob/lvmini/math/mathcore/src/NumericalDerivator.cxx - * - * This class attempts to more closely follow the Minuit2 implementation. - * Modified things (w.r.t. NumericalDerivator) are indicated by MODIFIED. - */ - -#include "NumericalDerivatorMinuit2.h" -#include -#include -#include -#include -#include -#include -#include "Fit/ParameterSettings.h" - -#include // needed here because in Fitter is only a forward declaration - -#include - -namespace RooFit { - -NumericalDerivatorMinuit2::NumericalDerivatorMinuit2(bool always_exactly_mimic_minuit2) - : _always_exactly_mimic_minuit2(always_exactly_mimic_minuit2) -{ -} - -NumericalDerivatorMinuit2::NumericalDerivatorMinuit2(double step_tolerance, double grad_tolerance, unsigned int ncycles, - double error_level, bool always_exactly_mimic_minuit2) - : fStepTolerance(step_tolerance), fGradTolerance(grad_tolerance), fNCycles(ncycles), Up(error_level), - _always_exactly_mimic_minuit2(always_exactly_mimic_minuit2) -{ -} - -// deep copy constructor -NumericalDerivatorMinuit2::NumericalDerivatorMinuit2(const RooFit::NumericalDerivatorMinuit2 &other) - : fStepTolerance(other.fStepTolerance), fGradTolerance(other.fGradTolerance), fNCycles(other.fNCycles), Up(other.Up), - fVal(other.fVal), vx(other.vx), vx_external(other.vx_external), dfmin(other.dfmin), vrysml(other.vrysml), - precision(other.precision), _always_exactly_mimic_minuit2(other._always_exactly_mimic_minuit2), - vx_fVal_cache(other.vx_fVal_cache) -{ -} - -void NumericalDerivatorMinuit2::SetStepTolerance(double value) -{ - fStepTolerance = value; -} - -void NumericalDerivatorMinuit2::SetGradTolerance(double value) -{ - fGradTolerance = value; -} - -void NumericalDerivatorMinuit2::SetNCycles(int value) -{ - fNCycles = value; -} - -NumericalDerivatorMinuit2::~NumericalDerivatorMinuit2() -{ - // TODO Auto-generated destructor stub -} - -// This function sets internal state based on input parameters. This state -// setup is used in the actual (partial) derivative calculations. -void NumericalDerivatorMinuit2::setup_differentiate(const ROOT::Math::IBaseFunctionMultiDim *function, const double *cx, - const std::vector ¶meters) -{ - assert(function != nullptr && "function is a nullptr"); - - if (vx.size() != function->NDim()) { - vx.resize(function->NDim()); - } - if (vx_external.size() != function->NDim()) { - vx_external.resize(function->NDim()); - } - if (vx_fVal_cache.size() != function->NDim()) { - vx_fVal_cache.resize(function->NDim()); - } - std::copy(cx, cx + function->NDim(), vx.data()); - - // convert to Minuit external parameters - for (unsigned i = 0; i < function->NDim(); i++) { - vx_external[i] = Int2ext(parameters[i], vx[i]); - } - - if (vx != vx_fVal_cache) { - vx_fVal_cache = vx; - fVal = (*function)(vx_external.data()); // value of function at given points - } - - dfmin = 8. * precision.Eps2() * (std::abs(fVal) + Up); - vrysml = 8. * precision.Eps() * precision.Eps(); -} - -MinuitDerivatorElement -NumericalDerivatorMinuit2::partial_derivative(const ROOT::Math::IBaseFunctionMultiDim *function, const double *x, - const std::vector ¶meters, - unsigned int i_component, MinuitDerivatorElement previous) -{ - setup_differentiate(function, x, parameters); - return fast_partial_derivative(function, parameters, i_component, previous); -} - -// leaves the parameter setup to the caller -MinuitDerivatorElement NumericalDerivatorMinuit2::fast_partial_derivative(const ROOT::Math::IBaseFunctionMultiDim *function, - const std::vector ¶meters, - unsigned int ix, const MinuitDerivatorElement& previous) -{ - MinuitDerivatorElement deriv {previous}; - - double xtf = vx[ix]; - double epspri = precision.Eps2() + std::abs(deriv.derivative * precision.Eps2()); - double step_old = 0.; - - for (unsigned int j = 0; j < fNCycles; ++j) { - double optstp = std::sqrt(dfmin / (std::abs(deriv.second_derivative) + epspri)); - double step = std::max(optstp, std::abs(0.1 * deriv.step_size)); - - if (parameters[ix].IsBound()) { - if (step > 0.5) - step = 0.5; - } - - double stpmax = 10. * std::abs(deriv.step_size); - if (step > stpmax) - step = stpmax; - - double stpmin = std::max(vrysml, 8. * std::abs(precision.Eps2() * vx[ix])); - if (step < stpmin) - step = stpmin; - if (std::abs((step - step_old) / step) < fStepTolerance) { - break; - } - deriv.step_size = step; - step_old = step; - vx[ix] = xtf + step; - vx_external[ix] = Int2ext(parameters[ix], vx[ix]); - double fs1 = (*function)(vx_external.data()); - - vx[ix] = xtf - step; - vx_external[ix] = Int2ext(parameters[ix], vx[ix]); - double fs2 = (*function)(vx_external.data()); - - vx[ix] = xtf; - vx_external[ix] = Int2ext(parameters[ix], vx[ix]); - - double fGrd_old = deriv.derivative; - deriv.derivative = 0.5 * (fs1 - fs2) / step; - - deriv.second_derivative = (fs1 + fs2 - 2. * fVal) / step / step; - - if (std::abs(fGrd_old - deriv.derivative) / (std::abs(deriv.derivative) + dfmin / step) < fGradTolerance) { - break; - } - } - return deriv; -} - -MinuitDerivatorElement -NumericalDerivatorMinuit2::operator()(const ROOT::Math::IBaseFunctionMultiDim *function, const double *x, - const std::vector ¶meters, - unsigned int i_component, const MinuitDerivatorElement& previous) -{ - return partial_derivative(function, x, parameters, i_component, previous); -} - -std::vector -NumericalDerivatorMinuit2::Differentiate(const ROOT::Math::IBaseFunctionMultiDim *function, const double *cx, - const std::vector ¶meters, - const std::vector& previous_gradient) -{ - setup_differentiate(function, cx, parameters); - - std::vector gradient; - gradient.reserve(function->NDim()); - - for (unsigned int ix = 0; ix < function->NDim(); ++ix) { - gradient.emplace_back(fast_partial_derivative(function, parameters, ix, previous_gradient[ix])); - } - - return gradient; -} - - -double NumericalDerivatorMinuit2::Int2ext(const ROOT::Fit::ParameterSettings ¶meter, double val) const -{ - // return external value from internal value for parameter i - if (parameter.IsBound()) { - if (parameter.IsDoubleBound()) { - return fDoubleLimTrafo.Int2ext(val, parameter.UpperLimit(), parameter.LowerLimit()); - } else if (parameter.HasUpperLimit() && !parameter.HasLowerLimit()) { - return fUpperLimTrafo.Int2ext(val, parameter.UpperLimit()); - } else { - return fLowerLimTrafo.Int2ext(val, parameter.LowerLimit()); - } - } - - return val; -} - -double NumericalDerivatorMinuit2::Ext2int(const ROOT::Fit::ParameterSettings ¶meter, double val) const -{ - // return the internal value for parameter i with external value val - - if (parameter.IsBound()) { - if (parameter.IsDoubleBound()) { - return fDoubleLimTrafo.Ext2int(val, parameter.UpperLimit(), parameter.LowerLimit(), precision); - } else if (parameter.HasUpperLimit() && !parameter.HasLowerLimit()) { - return fUpperLimTrafo.Ext2int(val, parameter.UpperLimit(), precision); - } else { - return fLowerLimTrafo.Ext2int(val, parameter.LowerLimit(), precision); - } - } - - return val; -} - -double NumericalDerivatorMinuit2::DInt2Ext(const ROOT::Fit::ParameterSettings ¶meter, 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 (parameter.IsBound()) { - if (parameter.IsDoubleBound()) { - dd = fDoubleLimTrafo.DInt2Ext(val, parameter.UpperLimit(), parameter.LowerLimit()); - } else if (parameter.HasUpperLimit() && !parameter.HasLowerLimit()) { - dd = fUpperLimTrafo.DInt2Ext(val, parameter.UpperLimit()); - } else { - dd = fLowerLimTrafo.DInt2Ext(val, parameter.LowerLimit()); - } - } - - return dd; -} - -// MODIFIED: -// This function was not implemented as in Minuit2. Now it copies the behavior -// of InitialGradientCalculator. See https://github.com/roofit-dev/root/issues/10 -void NumericalDerivatorMinuit2::SetInitialGradient(const ROOT::Math::IBaseFunctionMultiDim *function, - const std::vector ¶meters, - std::vector& gradient) -{ - // set an initial gradient using some given steps - // (used in the first iteration) - - assert(function != nullptr && "function is a nullptr"); - - double eps2 = precision.Eps2(); - - unsigned ix = 0; - for (auto parameter = parameters.begin(); parameter != parameters.end(); ++parameter, ++ix) { - // What Minuit2 calls "Error" is stepsize on the ROOT side. - double werr = parameter->StepSize(); - - // 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); - - if (_always_exactly_mimic_minuit2) { - // this transformation can lose a few bits, but Minuit2 does it too - sav = Int2ext(*parameter, var); - } - - double sav2 = sav + werr; - - if (parameter->HasUpperLimit() && sav2 > parameter->UpperLimit()) { - sav2 = parameter->UpperLimit(); - } - - double var2 = Ext2int(*parameter, sav2); - double vplu = var2 - var; - - sav2 = sav - werr; - - if (parameter->HasLowerLimit() && sav2 < parameter->LowerLimit()) { - sav2 = parameter->LowerLimit(); - } - - var2 = Ext2int(*parameter, sav2); - double vmin = var2 - 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 g2 = 2.0 * Up / (dirin * dirin); - double gstep = std::max(gsmin, 0.1 * dirin); - double grd = g2 * dirin; - - if (parameter->IsBound()) { - if (gstep > 0.5) - gstep = 0.5; - } - - gradient[ix].derivative = grd; - gradient[ix].second_derivative = g2; - gradient[ix].step_size = gstep; - } -} - -bool NumericalDerivatorMinuit2::always_exactly_mimic_minuit2() const -{ - return _always_exactly_mimic_minuit2; -}; - -void NumericalDerivatorMinuit2::set_always_exactly_mimic_minuit2(bool flag) -{ - _always_exactly_mimic_minuit2 = flag; -} - -void NumericalDerivatorMinuit2::set_step_tolerance(double step_tolerance) -{ - fStepTolerance = step_tolerance; -} -void NumericalDerivatorMinuit2::set_grad_tolerance(double grad_tolerance) -{ - fGradTolerance = grad_tolerance; -} -void NumericalDerivatorMinuit2::set_ncycles(unsigned int ncycles) -{ - fNCycles = ncycles; -} -void NumericalDerivatorMinuit2::set_error_level(double error_level) -{ - Up = error_level; -} - -std::ostream& operator<<(std::ostream& out, const MinuitDerivatorElement value){ - return out << "(derivative: " << value.derivative << ", second_derivative: " << value.second_derivative << ", step_size: " << value.step_size << ")"; -} - -} // namespace RooFit From f92404d83ce613aa1615f8a92938ab1aa4a1c927 Mon Sep 17 00:00:00 2001 From: Patrick Bos Date: Thu, 15 Jul 2021 14:56:55 +0200 Subject: [PATCH 5/7] memory efficiently reorder NumericalDerivator members Co-authored-by: Stephan Hageboeck --- math/minuit2/inc/Minuit2/NumericalDerivator.h | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/math/minuit2/inc/Minuit2/NumericalDerivator.h b/math/minuit2/inc/Minuit2/NumericalDerivator.h index a7343d9bd9f10..871537d13fd37 100644 --- a/math/minuit2/inc/Minuit2/NumericalDerivator.h +++ b/math/minuit2/inc/Minuit2/NumericalDerivator.h @@ -76,14 +76,18 @@ class NumericalDerivator { const std::vector ¶meters, std::vector &gradient); + bool AlwaysExactlyMimicMinuit2() const; + void SetAlwaysExactlyMimicMinuit2(bool flag); + private: double fStepTolerance = 0.5; double fGradTolerance = 0.1; - unsigned int fNCycles = 2; double fUp = 1; double fVal = 0; - std::vector fVx, fVxExternal; + std::vector fVx; + std::vector fVxExternal; + std::vector fVxFValCache; double fDfmin; double fVrysml; @@ -95,15 +99,10 @@ class NumericalDerivator { ROOT::Minuit2::SqrtUpParameterTransformation fUpperLimTrafo; ROOT::Minuit2::SqrtLowParameterTransformation fLowerLimTrafo; -private: + unsigned int fNCycles = 2; bool fAlwaysExactlyMimicMinuit2; -public: - bool AlwaysExactlyMimicMinuit2() const; - void SetAlwaysExactlyMimicMinuit2(bool flag); -private: - std::vector fVxFValCache; }; std::ostream &operator<<(std::ostream &out, const DerivatorElement &value); @@ -111,4 +110,4 @@ std::ostream &operator<<(std::ostream &out, const DerivatorElement &value); } // namespace Minuit2 } // namespace ROOT -#endif // ROOT_Minuit2_NumericalDerivator \ No newline at end of file +#endif // ROOT_Minuit2_NumericalDerivator From 46ff6046a51c1981c49f589026b016e777610dda Mon Sep 17 00:00:00 2001 From: Patrick Bos Date: Thu, 15 Jul 2021 15:28:06 +0200 Subject: [PATCH 6/7] Apply suggestions from code review Co-authored-by: Jonas Rembser Co-authored-by: Stephan Hageboeck --- math/minuit2/inc/Minuit2/NumericalDerivator.h | 2 -- math/minuit2/src/NumericalDerivator.cxx | 27 +++++++------------ 2 files changed, 10 insertions(+), 19 deletions(-) diff --git a/math/minuit2/inc/Minuit2/NumericalDerivator.h b/math/minuit2/inc/Minuit2/NumericalDerivator.h index 871537d13fd37..889ad9279e25f 100644 --- a/math/minuit2/inc/Minuit2/NumericalDerivator.h +++ b/math/minuit2/inc/Minuit2/NumericalDerivator.h @@ -18,9 +18,7 @@ #ifndef ROOT_Minuit2_NumericalDerivator #define ROOT_Minuit2_NumericalDerivator -#ifndef ROOT_Math_IFunctionfwd #include -#endif #include #include "Fit/ParameterSettings.h" diff --git a/math/minuit2/src/NumericalDerivator.cxx b/math/minuit2/src/NumericalDerivator.cxx index e099560ae7a6d..03939723573a8 100644 --- a/math/minuit2/src/NumericalDerivator.cxx +++ b/math/minuit2/src/NumericalDerivator.cxx @@ -6,8 +6,8 @@ * Copyright (c) 2013 , LCG ROOT MathLib Team * * * **********************************************************************/ -/* - * NumericalDerivator.cxx +/** + * \class NumericalDerivator * * Original version created on: Aug 14, 2013 * Authors: L. Moneta, J. T. Offermann @@ -78,22 +78,16 @@ void NumericalDerivator::SetErrorLevel(double value) fUp = value; } -// This function sets internal state based on input parameters. This state -// setup is used in the actual (partial) derivative calculations. +/// This function sets internal state based on input parameters. This state +/// setup is used in the actual (partial) derivative calculations. void NumericalDerivator::SetupDifferentiate(const ROOT::Math::IBaseFunctionMultiDim *function, const double *cx, const std::vector ¶meters) { assert(function != nullptr && "function is a nullptr"); - if (fVx.size() != function->NDim()) { - fVx.resize(function->NDim()); - } - if (fVxExternal.size() != function->NDim()) { - fVxExternal.resize(function->NDim()); - } - if (fVxFValCache.size() != function->NDim()) { - fVxFValCache.resize(function->NDim()); - } + fVx.resize(function->NDim()); + fVxExternal.resize(function->NDim()); + fVxFValCache.resize(function->NDim()); std::copy(cx, cx + function->NDim(), fVx.data()); // convert to Minuit external parameters @@ -251,8 +245,8 @@ double NumericalDerivator::DInt2Ext(const ROOT::Fit::ParameterSettings ¶mete } // MODIFIED: -// This function was not implemented as in Minuit2. Now it copies the behavior -// of InitialGradientCalculator. See https://github.com/roofit-dev/root/issues/10 +/// This function was not implemented as in Minuit2. Now it copies the behavior +/// of InitialGradientCalculator. See https://github.com/roofit-dev/root/issues/10 void NumericalDerivator::SetInitialGradient(const ROOT::Math::IBaseFunctionMultiDim *function, const std::vector ¶meters, std::vector &gradient) @@ -306,8 +300,7 @@ void NumericalDerivator::SetInitialGradient(const ROOT::Math::IBaseFunctionMulti double grd = g2 * dirin; if (parameter->IsBound()) { - if (gstep > 0.5) - gstep = 0.5; + gstep = std::min(gstep, 0.5); } gradient[ix].derivative = grd; From 435fe8fbd44318e2a89de45281054b15820119f1 Mon Sep 17 00:00:00 2001 From: "E. G. Patrick Bos" Date: Thu, 15 Jul 2021 15:30:06 +0200 Subject: [PATCH 7/7] apply more review comments From comments by @guitargeek and @hageboeck on PR #8567 (https://github.com/root-project/root/pull/8567): - Simplify copy constructor. - Inline simple getter/setters. - Replace fabs with std::abs. - Fix warning due to changed member ordering. --- math/minuit2/inc/Minuit2/NumericalDerivator.h | 15 +++---- math/minuit2/src/NumericalDerivator.cxx | 44 ++----------------- 2 files changed, 11 insertions(+), 48 deletions(-) diff --git a/math/minuit2/inc/Minuit2/NumericalDerivator.h b/math/minuit2/inc/Minuit2/NumericalDerivator.h index 889ad9279e25f..ee6203274f066 100644 --- a/math/minuit2/inc/Minuit2/NumericalDerivator.h +++ b/math/minuit2/inc/Minuit2/NumericalDerivator.h @@ -61,10 +61,10 @@ class NumericalDerivator { const DerivatorElement &previous); double GetValue() const { return fVal; } - void SetStepTolerance(double value); - void SetGradTolerance(double value); - void SetNCycles(unsigned int value); - void SetErrorLevel(double value); + inline void SetStepTolerance(double value) { fStepTolerance = value; } + inline void SetGradTolerance(double value) { fGradTolerance = value; } + inline void SetNCycles(unsigned int value) { fNCycles = value; } + inline void SetErrorLevel(double value) { fUp = value; } double Int2ext(const ROOT::Fit::ParameterSettings ¶meter, double val) const; double Ext2int(const ROOT::Fit::ParameterSettings ¶meter, double val) const; @@ -74,9 +74,9 @@ class NumericalDerivator { const std::vector ¶meters, std::vector &gradient); - bool AlwaysExactlyMimicMinuit2() const; - void SetAlwaysExactlyMimicMinuit2(bool flag); - + inline bool AlwaysExactlyMimicMinuit2() const { return fAlwaysExactlyMimicMinuit2; } + inline void SetAlwaysExactlyMimicMinuit2(bool flag) { fAlwaysExactlyMimicMinuit2 = flag; } + private: double fStepTolerance = 0.5; double fGradTolerance = 0.1; @@ -100,7 +100,6 @@ class NumericalDerivator { unsigned int fNCycles = 2; bool fAlwaysExactlyMimicMinuit2; - }; std::ostream &operator<<(std::ostream &out, const DerivatorElement &value); diff --git a/math/minuit2/src/NumericalDerivator.cxx b/math/minuit2/src/NumericalDerivator.cxx index 03939723573a8..6c2b3f2ac5ec0 100644 --- a/math/minuit2/src/NumericalDerivator.cxx +++ b/math/minuit2/src/NumericalDerivator.cxx @@ -44,39 +44,13 @@ NumericalDerivator::NumericalDerivator(bool always_exactly_mimic_minuit2) NumericalDerivator::NumericalDerivator(double step_tolerance, double grad_tolerance, unsigned int ncycles, double error_level, bool always_exactly_mimic_minuit2) - : fStepTolerance(step_tolerance), fGradTolerance(grad_tolerance), fNCycles(ncycles), fUp(error_level), + : fStepTolerance(step_tolerance), fGradTolerance(grad_tolerance), fUp(error_level), fNCycles(ncycles), fAlwaysExactlyMimicMinuit2(always_exactly_mimic_minuit2) { } -// deep copy constructor -NumericalDerivator::NumericalDerivator(const NumericalDerivator &other) - : fStepTolerance(other.fStepTolerance), fGradTolerance(other.fGradTolerance), fNCycles(other.fNCycles), - fUp(other.fUp), fVal(other.fVal), fVx(other.fVx), fVxExternal(other.fVxExternal), fDfmin(other.fDfmin), - fVrysml(other.fVrysml), fPrecision(other.fPrecision), fAlwaysExactlyMimicMinuit2(other.fAlwaysExactlyMimicMinuit2), - fVxFValCache(other.fVxFValCache) -{ -} - -void NumericalDerivator::SetStepTolerance(double value) -{ - fStepTolerance = value; -} - -void NumericalDerivator::SetGradTolerance(double value) -{ - fGradTolerance = value; -} - -void NumericalDerivator::SetNCycles(unsigned int value) -{ - fNCycles = value; -} +NumericalDerivator::NumericalDerivator(const NumericalDerivator &/*other*/) = default; -void NumericalDerivator::SetErrorLevel(double value) -{ - fUp = value; -} /// This function sets internal state based on input parameters. This state /// setup is used in the actual (partial) derivative calculations. @@ -292,9 +266,9 @@ void NumericalDerivator::SetInitialGradient(const ROOT::Math::IBaseFunctionMulti var2 = Ext2int(*parameter, sav2); double vmin = var2 - var; - double gsmin = 8. * eps2 * (fabs(var) + eps2); + double gsmin = 8. * eps2 * (std::abs(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 dirin = std::max(0.5 * (std::abs(vplu) + std::abs(vmin)), gsmin); double g2 = 2.0 * fUp / (dirin * dirin); double gstep = std::max(gsmin, 0.1 * dirin); double grd = g2 * dirin; @@ -309,16 +283,6 @@ void NumericalDerivator::SetInitialGradient(const ROOT::Math::IBaseFunctionMulti } } -bool NumericalDerivator::AlwaysExactlyMimicMinuit2() const -{ - return fAlwaysExactlyMimicMinuit2; -}; - -void NumericalDerivator::SetAlwaysExactlyMimicMinuit2(bool flag) -{ - fAlwaysExactlyMimicMinuit2 = flag; -} - std::ostream &operator<<(std::ostream &out, const DerivatorElement &value) { return out << "(derivative: " << value.derivative << ", second_derivative: " << value.second_derivative