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..ee6203274f066 --- /dev/null +++ b/math/minuit2/inc/Minuit2/NumericalDerivator.h @@ -0,0 +1,110 @@ +// @(#)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 + +#include + +#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; } + 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; + double DInt2Ext(const ROOT::Fit::ParameterSettings ¶meter, double val) const; + + void SetInitialGradient(const ROOT::Math::IBaseFunctionMultiDim *function, + const std::vector ¶meters, + std::vector &gradient); + + inline bool AlwaysExactlyMimicMinuit2() const { return fAlwaysExactlyMimicMinuit2; } + inline void SetAlwaysExactlyMimicMinuit2(bool flag) { fAlwaysExactlyMimicMinuit2 = flag; } + +private: + double fStepTolerance = 0.5; + double fGradTolerance = 0.1; + double fUp = 1; + double fVal = 0; + + std::vector fVx; + std::vector fVxExternal; + std::vector fVxFValCache; + 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; + + unsigned int fNCycles = 2; + bool fAlwaysExactlyMimicMinuit2; + +}; + +std::ostream &operator<<(std::ostream &out, const DerivatorElement &value); + +} // namespace Minuit2 +} // namespace ROOT + +#endif // ROOT_Minuit2_NumericalDerivator diff --git a/math/minuit2/src/NumericalDerivator.cxx b/math/minuit2/src/NumericalDerivator.cxx new file mode 100644 index 0000000000000..6c2b3f2ac5ec0 --- /dev/null +++ b/math/minuit2/src/NumericalDerivator.cxx @@ -0,0 +1,293 @@ +// @(#)root/mathcore:$Id$ +// Authors: L. Moneta, J.T. Offermann, E.G.P. Bos 2013-2018 +// +/********************************************************************** + * * + * Copyright (c) 2013 , LCG ROOT MathLib Team * + * * + **********************************************************************/ +/** + * \class NumericalDerivator + * + * 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), fUp(error_level), fNCycles(ncycles), + fAlwaysExactlyMimicMinuit2(always_exactly_mimic_minuit2) +{ +} + +NumericalDerivator::NumericalDerivator(const NumericalDerivator &/*other*/) = default; + + +/// 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"); + + 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 + 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 * (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 * (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; + + if (parameter->IsBound()) { + gstep = std::min(gstep, 0.5); + } + + gradient[ix].derivative = grd; + gradient[ix].second_derivative = g2; + gradient[ix].step_size = gstep; + } +} + +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 afe21fe236a6f..371fad77eccea 100644 --- a/roofit/roofitcore/CMakeLists.txt +++ b/roofit/roofitcore/CMakeLists.txt @@ -454,6 +454,7 @@ ROOT_STANDARD_LIBRARY_PACKAGE(RooFitCore Matrix Tree Minuit + Minuit2 RIO MathCore Foam