Logo ROOT  
Reference Guide
NumericalDerivator.cxx
Go to the documentation of this file.
1// @(#)root/mathcore:$Id$
2// Authors: L. Moneta, J.T. Offermann, E.G.P. Bos 2013-2018
3//
4/**********************************************************************
5 * *
6 * Copyright (c) 2013 , LCG ROOT MathLib Team *
7 * *
8 **********************************************************************/
9/**
10 * \class NumericalDerivator
11 *
12 * Original version created on: Aug 14, 2013
13 * Authors: L. Moneta, J. T. Offermann
14 * Modified version created on: Sep 27, 2017
15 * Author: E. G. P. Bos
16 *
17 * NumericalDerivator was essentially a slightly modified copy of code
18 * written by M. Winkler, F. James, L. Moneta, and A. Zsenei for Minuit2,
19 * Copyright (c) 2005 LCG ROOT Math team, CERN/PH-SFT. Original version:
20 * https://github.com/lmoneta/root/blob/lvmini/math/mathcore/src/NumericalDerivator.cxx
21 *
22 * This class attempts to more closely follow the Minuit2 implementation.
23 * Modified things (w.r.t. original) are indicated by MODIFIED.
24 */
25
27#include <cmath>
28#include <algorithm>
29#include <Math/IFunction.h>
30#include <iostream>
31#include <TMath.h>
32#include <cassert>
34
35#include <Math/Minimizer.h> // needed here because in Fitter is only a forward declaration
36
37namespace ROOT {
38namespace Minuit2 {
39
40NumericalDerivator::NumericalDerivator(bool always_exactly_mimic_minuit2)
41 : fAlwaysExactlyMimicMinuit2(always_exactly_mimic_minuit2)
42{
43}
44
45NumericalDerivator::NumericalDerivator(double step_tolerance, double grad_tolerance, unsigned int ncycles,
46 double error_level, bool always_exactly_mimic_minuit2)
47 : fStepTolerance(step_tolerance), fGradTolerance(grad_tolerance), fUp(error_level), fNCycles(ncycles),
48 fAlwaysExactlyMimicMinuit2(always_exactly_mimic_minuit2)
49{
50}
51
53
54
55/// This function sets internal state based on input parameters. This state
56/// setup is used in the actual (partial) derivative calculations.
58 const std::vector<ROOT::Fit::ParameterSettings> &parameters)
59{
60 assert(function != nullptr && "function is a nullptr");
61
62 fVx.resize(function->NDim());
63 fVxExternal.resize(function->NDim());
64 fVxFValCache.resize(function->NDim());
65 std::copy(cx, cx + function->NDim(), fVx.data());
66
67 // convert to Minuit external parameters
68 for (unsigned i = 0; i < function->NDim(); i++) {
69 fVxExternal[i] = Int2ext(parameters[i], fVx[i]);
70 }
71
72 if (fVx != fVxFValCache) {
74 fVal = (*function)(fVxExternal.data()); // value of function at given points
75 }
76
77 fDfmin = 8. * fPrecision.Eps2() * (std::abs(fVal) + fUp);
78 fVrysml = 8. * fPrecision.Eps() * fPrecision.Eps();
79}
80
82 const double *x,
83 const std::vector<ROOT::Fit::ParameterSettings> &parameters,
84 unsigned int i_component, DerivatorElement previous)
85{
86 SetupDifferentiate(function, x, parameters);
87 return FastPartialDerivative(function, parameters, i_component, previous);
88}
89
90// leaves the parameter setup to the caller
92 const std::vector<ROOT::Fit::ParameterSettings> &parameters,
93 unsigned int i_component, const DerivatorElement &previous)
94{
95 DerivatorElement deriv{previous};
96
97 double xtf = fVx[i_component];
98 double epspri = fPrecision.Eps2() + std::abs(deriv.derivative * fPrecision.Eps2());
99 double step_old = 0.;
100
101 for (unsigned int j = 0; j < fNCycles; ++j) {
102 double optstp = std::sqrt(fDfmin / (std::abs(deriv.second_derivative) + epspri));
103 double step = std::max(optstp, std::abs(0.1 * deriv.step_size));
104
105 if (parameters[i_component].IsBound()) {
106 if (step > 0.5)
107 step = 0.5;
108 }
109
110 double stpmax = 10. * std::abs(deriv.step_size);
111 if (step > stpmax)
112 step = stpmax;
113
114 double stpmin = std::max(fVrysml, 8. * std::abs(fPrecision.Eps2() * fVx[i_component]));
115 if (step < stpmin)
116 step = stpmin;
117 if (std::abs((step - step_old) / step) < fStepTolerance) {
118 break;
119 }
120 deriv.step_size = step;
121 step_old = step;
122 fVx[i_component] = xtf + step;
123 fVxExternal[i_component] = Int2ext(parameters[i_component], fVx[i_component]);
124 double fs1 = (*function)(fVxExternal.data());
125
126 fVx[i_component] = xtf - step;
127 fVxExternal[i_component] = Int2ext(parameters[i_component], fVx[i_component]);
128 double fs2 = (*function)(fVxExternal.data());
129
130 fVx[i_component] = xtf;
131 fVxExternal[i_component] = Int2ext(parameters[i_component], fVx[i_component]);
132
133 double fGrd_old = deriv.derivative;
134 deriv.derivative = 0.5 * (fs1 - fs2) / step;
135
136 deriv.second_derivative = (fs1 + fs2 - 2. * fVal) / step / step;
137
138 if (std::abs(fGrd_old - deriv.derivative) / (std::abs(deriv.derivative) + fDfmin / step) < fGradTolerance) {
139 break;
140 }
141 }
142 return deriv;
143}
144
146 const std::vector<ROOT::Fit::ParameterSettings> &parameters,
147 unsigned int i_component, const DerivatorElement &previous)
148{
149 return PartialDerivative(function, x, parameters, i_component, previous);
150}
151
152std::vector<DerivatorElement>
154 const std::vector<ROOT::Fit::ParameterSettings> &parameters,
155 const std::vector<DerivatorElement> &previous_gradient)
156{
157 SetupDifferentiate(function, cx, parameters);
158
159 std::vector<DerivatorElement> gradient;
160 gradient.reserve(function->NDim());
161
162 for (unsigned int ix = 0; ix < function->NDim(); ++ix) {
163 gradient.emplace_back(FastPartialDerivative(function, parameters, ix, previous_gradient[ix]));
164 }
165
166 return gradient;
167}
168
169double NumericalDerivator::Int2ext(const ROOT::Fit::ParameterSettings &parameter, double val) const
170{
171 // return external value from internal value for parameter i
172 if (parameter.IsBound()) {
173 if (parameter.IsDoubleBound()) {
174 return fDoubleLimTrafo.Int2ext(val, parameter.UpperLimit(), parameter.LowerLimit());
175 } else if (parameter.HasUpperLimit() && !parameter.HasLowerLimit()) {
176 return fUpperLimTrafo.Int2ext(val, parameter.UpperLimit());
177 } else {
178 return fLowerLimTrafo.Int2ext(val, parameter.LowerLimit());
179 }
180 }
181
182 return val;
183}
184
185double NumericalDerivator::Ext2int(const ROOT::Fit::ParameterSettings &parameter, double val) const
186{
187 // return the internal value for parameter i with external value val
188
189 if (parameter.IsBound()) {
190 if (parameter.IsDoubleBound()) {
191 return fDoubleLimTrafo.Ext2int(val, parameter.UpperLimit(), parameter.LowerLimit(), fPrecision);
192 } else if (parameter.HasUpperLimit() && !parameter.HasLowerLimit()) {
193 return fUpperLimTrafo.Ext2int(val, parameter.UpperLimit(), fPrecision);
194 } else {
195 return fLowerLimTrafo.Ext2int(val, parameter.LowerLimit(), fPrecision);
196 }
197 }
198
199 return val;
200}
201
202double NumericalDerivator::DInt2Ext(const ROOT::Fit::ParameterSettings &parameter, double val) const
203{
204 // return the derivative of the int->ext transformation: dPext(i) / dPint(i)
205 // for the parameter i with value val
206
207 double dd = 1.;
208 if (parameter.IsBound()) {
209 if (parameter.IsDoubleBound()) {
210 dd = fDoubleLimTrafo.DInt2Ext(val, parameter.UpperLimit(), parameter.LowerLimit());
211 } else if (parameter.HasUpperLimit() && !parameter.HasLowerLimit()) {
212 dd = fUpperLimTrafo.DInt2Ext(val, parameter.UpperLimit());
213 } else {
214 dd = fLowerLimTrafo.DInt2Ext(val, parameter.LowerLimit());
215 }
216 }
217
218 return dd;
219}
220
221// MODIFIED:
222/// This function was not implemented as in Minuit2. Now it copies the behavior
223/// of InitialGradientCalculator. See https://github.com/roofit-dev/root/issues/10
225 const std::vector<ROOT::Fit::ParameterSettings> &parameters,
226 std::vector<DerivatorElement> &gradient)
227{
228 // set an initial gradient using some given steps
229 // (used in the first iteration)
230
231 double eps2 = fPrecision.Eps2();
232
233 unsigned ix = 0;
234 for (auto parameter = parameters.begin(); parameter != parameters.end(); ++parameter, ++ix) {
235 // What Minuit2 calls "Error" is stepsize on the ROOT side.
236 double werr = parameter->StepSize();
237
238 // Actually, sav in Minuit2 is the external parameter value, so that is
239 // what we called var before and var is unnecessary here.
240 double sav = parameter->Value();
241
242 // However, we do need var below, so let's calculate it using Ext2int:
243 double var = Ext2int(*parameter, sav);
244
246 // this transformation can lose a few bits, but Minuit2 does it too
247 sav = Int2ext(*parameter, var);
248 }
249
250 double sav2 = sav + werr;
251
252 if (parameter->HasUpperLimit() && sav2 > parameter->UpperLimit()) {
253 sav2 = parameter->UpperLimit();
254 }
255
256 double var2 = Ext2int(*parameter, sav2);
257 double vplu = var2 - var;
258
259 sav2 = sav - werr;
260
261 if (parameter->HasLowerLimit() && sav2 < parameter->LowerLimit()) {
262 sav2 = parameter->LowerLimit();
263 }
264
265 var2 = Ext2int(*parameter, sav2);
266 double vmin = var2 - var;
267 double gsmin = 8. * eps2 * (std::abs(var) + eps2);
268 // protect against very small step sizes which can cause dirin to zero and then nan values in grd
269 double dirin = std::max(0.5 * (std::abs(vplu) + std::abs(vmin)), gsmin);
270 double g2 = 2.0 * fUp / (dirin * dirin);
271 double gstep = std::max(gsmin, 0.1 * dirin);
272 double grd = g2 * dirin;
273
274 if (parameter->IsBound()) {
275 gstep = std::min(gstep, 0.5);
276 }
277
278 gradient[ix].derivative = grd;
279 gradient[ix].second_derivative = g2;
280 gradient[ix].step_size = gstep;
281 }
282}
283
284std::ostream &operator<<(std::ostream &out, const DerivatorElement &value)
285{
286 return out << "(derivative: " << value.derivative << ", second_derivative: " << value.second_derivative
287 << ", step_size: " << value.step_size << ")";
288}
289
290} // namespace Minuit2
291} // namespace ROOT
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void value
Class, describing value, limits and step size of the parameters Provides functionality also to set/re...
bool HasUpperLimit() const
check if parameter has upper limit
double LowerLimit() const
return lower limit value
bool HasLowerLimit() const
check if parameter has lower limit
double UpperLimit() const
return upper limit value
bool IsDoubleBound() const
check if is double bound (upper AND lower limit)
bool IsBound() const
check if is bound
double Eps() const
eps returns the smallest possible number so that 1.+eps > 1.
double Eps2() const
eps2 returns 2*sqrt(eps)
DerivatorElement operator()(const ROOT::Math::IBaseFunctionMultiDim *function, const double *x, const std::vector< ROOT::Fit::ParameterSettings > &parameters, unsigned int i_component, const DerivatorElement &previous)
double Int2ext(const ROOT::Fit::ParameterSettings &parameter, double val) const
ROOT::Minuit2::SqrtUpParameterTransformation fUpperLimTrafo
ROOT::Minuit2::MnMachinePrecision fPrecision
DerivatorElement PartialDerivative(const ROOT::Math::IBaseFunctionMultiDim *function, const double *x, const std::vector< ROOT::Fit::ParameterSettings > &parameters, unsigned int i_component, DerivatorElement previous)
void SetupDifferentiate(const ROOT::Math::IBaseFunctionMultiDim *function, const double *cx, const std::vector< ROOT::Fit::ParameterSettings > &parameters)
This function sets internal state based on input parameters.
void SetInitialGradient(const ROOT::Math::IBaseFunctionMultiDim *function, const std::vector< ROOT::Fit::ParameterSettings > &parameters, std::vector< DerivatorElement > &gradient)
This function was not implemented as in Minuit2.
NumericalDerivator(bool always_exactly_mimic_minuit2=true)
double DInt2Ext(const ROOT::Fit::ParameterSettings &parameter, double val) const
double Ext2int(const ROOT::Fit::ParameterSettings &parameter, double val) const
ROOT::Minuit2::SinParameterTransformation fDoubleLimTrafo
DerivatorElement FastPartialDerivative(const ROOT::Math::IBaseFunctionMultiDim *function, const std::vector< ROOT::Fit::ParameterSettings > &parameters, unsigned int i_component, const DerivatorElement &previous)
ROOT::Minuit2::SqrtLowParameterTransformation fLowerLimTrafo
std::vector< DerivatorElement > Differentiate(const ROOT::Math::IBaseFunctionMultiDim *function, const double *x, const std::vector< ROOT::Fit::ParameterSettings > &parameters, const std::vector< DerivatorElement > &previous_gradient)
long double DInt2Ext(long double Value, long double Upper, long double Lower) const
long double Int2ext(long double Value, long double Upper, long double Lower) const
long double Ext2int(long double Value, long double Upper, long double Lower, const MnMachinePrecision &) const
long double Ext2int(long double Value, long double Lower, const MnMachinePrecision &) const
long double DInt2Ext(long double Value, long double Lower) const
long double Int2ext(long double Value, long double Lower) const
long double Int2ext(long double Value, long double Upper) const
long double DInt2Ext(long double Value, long double Upper) const
long double Ext2int(long double Value, long double Upper, const MnMachinePrecision &) const
RVec< PromoteType< T > > abs(const RVec< T > &v)
Definition: RVec.hxx:1756
Double_t x[n]
Definition: legend1.C:17
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
std::ostream & operator<<(std::ostream &, const FunctionMinimum &)
Definition: MnPrint.cxx:335
void function(const Char_t *name_, T fun, const Char_t *docstring=0)
Definition: RExports.h:167
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.