1 // @(#)root/minuit2:$Id$
2 // Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005
3
4 /**********************************************************************
5 * *
6 * Copyright (c) 2005 LCG ROOT Math team, CERN/PH-SFT *
7 * *
8 **********************************************************************/
9
10 #include "Minuit2/InitialGradientCalculator.h"
11 #include "Minuit2/MnFcn.h"
12 #include "Minuit2/MnUserTransformation.h"
13 #include "Minuit2/MnMachinePrecision.h"
14 #include "Minuit2/MinimumParameters.h"
15 #include "Minuit2/FunctionGradient.h"
16 #include "Minuit2/MnStrategy.h"
17 #include "Minuit2/MnPrint.h"
18
19 #include <cmath>
20
21 namespace ROOT {
22
23 namespace Minuit2 {
24
operator ()(const MinimumParameters & par) const25 FunctionGradient InitialGradientCalculator::operator()(const MinimumParameters &par) const
26 {
27 // initial rough estimate of the gradient using the parameter step size
28
29 assert(par.IsValid());
30
31 unsigned int n = Trafo().VariableParameters();
32 assert(n == par.Vec().size());
33
34 MnPrint print("InitialGradientCalculator");
35
36 print.Debug("Calculating initial gradient at point", par.Vec());
37
38 MnAlgebraicVector gr(n), gr2(n), gst(n);
39
40 for (unsigned int i = 0; i < n; i++) {
41 unsigned int exOfIn = Trafo().ExtOfInt(i);
42
43 double var = par.Vec()(i);
44 double werr = Trafo().Parameter(exOfIn).Error();
45 double sav = Trafo().Int2ext(i, var);
46 double sav2 = sav + werr;
47 if (Trafo().Parameter(exOfIn).HasLimits()) {
48 if (Trafo().Parameter(exOfIn).HasUpperLimit() && sav2 > Trafo().Parameter(exOfIn).UpperLimit())
49 sav2 = Trafo().Parameter(exOfIn).UpperLimit();
50 }
51 double var2 = Trafo().Ext2int(exOfIn, sav2);
52 double vplu = var2 - var;
53 sav2 = sav - werr;
54 if (Trafo().Parameter(exOfIn).HasLimits()) {
55 if (Trafo().Parameter(exOfIn).HasLowerLimit() && sav2 < Trafo().Parameter(exOfIn).LowerLimit())
56 sav2 = Trafo().Parameter(exOfIn).LowerLimit();
57 }
58 var2 = Trafo().Ext2int(exOfIn, sav2);
59 double vmin = var2 - var;
60 double gsmin = 8. * Precision().Eps2() * (std::fabs(var) + Precision().Eps2());
61 // protect against very small step sizes which can cause dirin to zero and then nan values in grd
62 double dirin = std::max(0.5 * (std::fabs(vplu) + std::fabs(vmin)), gsmin);
63 double g2 = 2.0 * fFcn.ErrorDef() / (dirin * dirin);
64 double gstep = std::max(gsmin, 0.1 * dirin);
65 double grd = g2 * dirin;
66 if (Trafo().Parameter(exOfIn).HasLimits()) {
67 if (gstep > 0.5)
68 gstep = 0.5;
69 }
70 gr(i) = grd;
71 gr2(i) = g2;
72 gst(i) = gstep;
73
74 print.Debug("Computed initial gradient for parameter", Trafo().Name(exOfIn), "value", var, "[", vmin, ",", vplu,
75 "]", "dirin", dirin, "grd", grd, "g2", g2);
76 }
77
78 return FunctionGradient(gr, gr2, gst);
79 }
80
operator ()(const MinimumParameters & par,const FunctionGradient &) const81 FunctionGradient InitialGradientCalculator::operator()(const MinimumParameters &par, const FunctionGradient &) const
82 {
83 // Base class interface
84 return (*this)(par);
85 }
86
Precision() const87 const MnMachinePrecision &InitialGradientCalculator::Precision() const
88 {
89 // return precision (is set in trasformation class)
90 return fTransformation.Precision();
91 }
92
Ncycle() const93 unsigned int InitialGradientCalculator::Ncycle() const
94 {
95 // return ncyles (from Strategy)
96 return Strategy().GradientNCycles();
97 }
98
StepTolerance() const99 double InitialGradientCalculator::StepTolerance() const
100 {
101 // return Gradient step tolerance (from Strategy)
102 return Strategy().GradientStepTolerance();
103 }
104
GradTolerance() const105 double InitialGradientCalculator::GradTolerance() const
106 {
107 // return Gradient tolerance
108 return Strategy().GradientTolerance();
109 }
110
111 } // namespace Minuit2
112
113 } // namespace ROOT
114