ROOT   Reference Guide
Go to the documentation of this file.
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
11#include "Minuit2/MnFcn.h"
16#include "Minuit2/MnStrategy.h"
17
18#include <math.h>
19
20//#define DEBUG
21
22#if defined(DEBUG) || defined(WARNINGMSG)
23#include "Minuit2/MnPrint.h"
24#endif
25
26
27
28namespace ROOT {
29
30 namespace Minuit2 {
31
32
34 // initial rough estimate of the gradient using the parameter step size
35
36 assert(par.IsValid());
37
38 unsigned int n = Trafo().VariableParameters();
39 assert(n == par.Vec().size());
40
41#ifdef DEBUG
42 std::cout << "Initial gradient calculator - params " << par.Vec() << std::endl;
43#endif
44
45 MnAlgebraicVector gr(n), gr2(n), gst(n);
46
47 for(unsigned int i = 0; i < n; i++) {
48 unsigned int exOfIn = Trafo().ExtOfInt(i);
49
50 double var = par.Vec()(i);
51 double werr = Trafo().Parameter(exOfIn).Error();
52 double sav = Trafo().Int2ext(i, var);
53 double sav2 = sav + werr;
54 if(Trafo().Parameter(exOfIn).HasLimits()) {
55 if(Trafo().Parameter(exOfIn).HasUpperLimit() &&
56 sav2 > Trafo().Parameter(exOfIn).UpperLimit())
57 sav2 = Trafo().Parameter(exOfIn).UpperLimit();
58 }
59 double var2 = Trafo().Ext2int(exOfIn, sav2);
60 double vplu = var2 - var;
61 sav2 = sav - werr;
62 if(Trafo().Parameter(exOfIn).HasLimits()) {
63 if(Trafo().Parameter(exOfIn).HasLowerLimit() &&
64 sav2 < Trafo().Parameter(exOfIn).LowerLimit())
65 sav2 = Trafo().Parameter(exOfIn).LowerLimit();
66 }
67 var2 = Trafo().Ext2int(exOfIn, sav2);
68 double vmin = var2 - var;
69 double gsmin = 8.*Precision().Eps2()*(fabs(var) + Precision().Eps2());
70 // protect against very small step sizes which can cause dirin to zero and then nan values in grd
71 double dirin = std::max(0.5*(fabs(vplu) + fabs(vmin)), gsmin );
72 double g2 = 2.0*fFcn.ErrorDef()/(dirin*dirin);
73 double gstep = std::max(gsmin, 0.1*dirin);
74 double grd = g2*dirin;
75 if(Trafo().Parameter(exOfIn).HasLimits()) {
76 if(gstep > 0.5) gstep = 0.5;
77 }
78 gr(i) = grd;
79 gr2(i) = g2;
80 gst(i) = gstep;
81
82#ifdef DEBUG
83 std::cout << "computing initial gradient for parameter " << Trafo().Name(exOfIn) << " value = " << var
84 << " [ " << vmin << " , " << vplu << " ] " << "dirin " << dirin << " grd " << grd << " g2 " << g2 << std::endl;
85#endif
86
87 }
88
90}
91
93 // Base class interface
94 return (*this)(par);
95}
96
98 // return precision (is set in trasformation class)
100}
101
103 // return ncyles (from Strategy)
105}
106
108 // return Gradient step tolerance (from Strategy)
110}
111
115}
116
117
118 } // namespace Minuit2
119
120} // namespace ROOT
const MnMachinePrecision & Precision() const
const MnUserTransformation & Trafo() const
virtual FunctionGradient operator()(const MinimumParameters &) const
unsigned int size() const
Definition: LAVector.h:198
const MnAlgebraicVector & Vec() const
double ErrorDef() const
Definition: MnFcn.cxx:33
determines the relative floating point arithmetic precision.
double Eps2() const
eps2 returns 2*sqrt(eps)
Definition: MnStrategy.h:42
Definition: MnStrategy.h:43
Definition: MnStrategy.h:41
double Ext2int(unsigned int, double) const
unsigned int ExtOfInt(unsigned int internal) const
const char * Name(unsigned int) const
double Int2ext(unsigned int, double) const
const MnMachinePrecision & Precision() const
forwarded interface
const MinuitParameter & Parameter(unsigned int) const
const Int_t n
Definition: legend1.C:16
TGraphErrors * gr
Definition: legend1.C:25
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
VSD Structures.
Definition: StringConv.hxx:21