Logo ROOT   6.08/07
Reference Guide
InitialGradientCalculator.cxx
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 
28 namespace 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 
89  return FunctionGradient(gr, gr2, gst);
90 }
91 
93  // Base class interface
94  return (*this)(par);
95 }
96 
98  // return precision (is set in trasformation class)
99  return fTransformation.Precision();
100 }
101 
102 unsigned int InitialGradientCalculator::Ncycle() const {
103  // return ncyles (from Strategy)
104  return Strategy().GradientNCycles();
105 }
106 
108  // return Gradient step tolerance (from Strategy)
109  return Strategy().GradientStepTolerance();
110 }
111 
113  // return Gradient tolerance
114  return Strategy().GradientTolerance();
115 }
116 
117 
118  } // namespace Minuit2
119 
120 } // namespace ROOT
double par[1]
Definition: unuranDistr.cxx:38
This namespace contains pre-defined functions to be used in conjuction with TExecutor::Map and TExecu...
Definition: StringConv.hxx:21
const MnAlgebraicVector & Vec() const
double GradientTolerance() const
Definition: MnStrategy.h:43
double Int2ext(unsigned int, double) const
unsigned int GradientNCycles() const
Definition: MnStrategy.h:41
double Ext2int(unsigned int, double) const
unsigned int ExtOfInt(unsigned int internal) const
determines the relative floating point arithmetic precision.
double GradientStepTolerance() const
Definition: MnStrategy.h:42
const MnMachinePrecision & Precision() const
forwarded interface
virtual FunctionGradient operator()(const MinimumParameters &) const
const char * Name(unsigned int) const
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
unsigned int size() const
Definition: LAVector.h:198
TGraphErrors * gr
Definition: legend1.C:25
double Eps2() const
eps2 returns 2*sqrt(eps)
const MnMachinePrecision & Precision() const
const MinuitParameter & Parameter(unsigned int) const
double ErrorDef() const
Definition: MnFcn.cxx:33
const MnUserTransformation & Trafo() const
const Int_t n
Definition: legend1.C:16