ROOT  6.06/09
Reference Guide
HessianGradientCalculator.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 
12 #include "Minuit2/MnFcn.h"
17 #include "Minuit2/MnStrategy.h"
18 
19 #include <math.h>
20 
21 //#define DEBUG
22 
23 #if defined(DEBUG) || defined(WARNINGMSG)
24 #include "Minuit2/MnPrint.h"
25 #endif
26 
27 #include "Minuit2/MPIProcess.h"
28 
29 namespace ROOT {
30 
31  namespace Minuit2 {
32 
33 
35  // use initial gradient as starting point
37  FunctionGradient gra = gc(par);
38 
39  return (*this)(par, gra);
40 }
41 
43  // interface of the base class. Use DeltaGradient for op.
44  std::pair<FunctionGradient, MnAlgebraicVector> mypair = DeltaGradient(par, Gradient);
45 
46  return mypair.first;
47 }
48 
50  // return the precision
51  return fTransformation.Precision();
52 }
53 
54 unsigned int HessianGradientCalculator::Ncycle() const {
55  // return number of calculation cycles (defined in strategy)
57 }
58 
60  // return tolerance on step size (defined in strategy)
62 }
63 
65  // return gradient tolerance (defines in strategy)
66  return Strategy().GradientTolerance();
67 }
68 
69 std::pair<FunctionGradient, MnAlgebraicVector> HessianGradientCalculator::DeltaGradient(const MinimumParameters& par, const FunctionGradient& Gradient) const {
70  // calculate gradient for Hessian
71  assert(par.IsValid());
72 
73  MnAlgebraicVector x = par.Vec();
74  MnAlgebraicVector grd = Gradient.Grad();
75  const MnAlgebraicVector& g2 = Gradient.G2();
76  //const MnAlgebraicVector& gstep = Gradient.Gstep();
77  // update also gradient step sizes
78  MnAlgebraicVector gstep = Gradient.Gstep();
79 
80  double fcnmin = par.Fval();
81  // std::cout<<"fval: "<<fcnmin<<std::endl;
82 
83  double dfmin = 4.*Precision().Eps2()*(fabs(fcnmin)+Fcn().Up());
84 
85  unsigned int n = x.size();
86  MnAlgebraicVector dgrd(n);
87 
88  MPIProcess mpiproc(n,0);
89  // initial starting values
90  unsigned int startElementIndex = mpiproc.StartElementIndex();
91  unsigned int endElementIndex = mpiproc.EndElementIndex();
92 
93  for(unsigned int i = startElementIndex; i < endElementIndex; i++) {
94  double xtf = x(i);
95  double dmin = 4.*Precision().Eps2()*(xtf + Precision().Eps2());
96  double epspri = Precision().Eps2() + fabs(grd(i)*Precision().Eps2());
97  double optstp = sqrt(dfmin/(fabs(g2(i))+epspri));
98  double d = 0.2*fabs(gstep(i));
99  if(d > optstp) d = optstp;
100  if(d < dmin) d = dmin;
101  double chgold = 10000.;
102  double dgmin = 0.;
103  double grdold = 0.;
104  double grdnew = 0.;
105  for(unsigned int j = 0; j < Ncycle(); j++) {
106  x(i) = xtf + d;
107  double fs1 = Fcn()(x);
108  x(i) = xtf - d;
109  double fs2 = Fcn()(x);
110  x(i) = xtf;
111  // double sag = 0.5*(fs1+fs2-2.*fcnmin);
112  //LM: should I calculate also here second derivatives ???
113 
114  grdold = grd(i);
115  grdnew = (fs1-fs2)/(2.*d);
116  dgmin = Precision().Eps()*(fabs(fs1) + fabs(fs2))/d;
117  //if(fabs(grdnew) < Precision().Eps()) break;
118  if (grdnew == 0) break;
119  double change = fabs((grdold-grdnew)/grdnew);
120  if(change > chgold && j > 1) break;
121  chgold = change;
122  grd(i) = grdnew;
123  //LM : update also the step sizes
124  gstep(i) = d;
125 
126  if(change < 0.05) break;
127  if(fabs(grdold-grdnew) < dgmin) break;
128  if(d < dmin) break;
129  d *= 0.2;
130  }
131 
132  dgrd(i) = std::max(dgmin, fabs(grdold-grdnew));
133 
134 #ifdef DEBUG
135  std::cout << "HGC Param : " << i << "\t new g1 = " << grd(i) << " gstep = " << d << " dgrd = " << dgrd(i) << std::endl;
136 #endif
137 
138  }
139 
140  mpiproc.SyncVector(grd);
141  mpiproc.SyncVector(gstep);
142  mpiproc.SyncVector(dgrd);
143 
144  return std::pair<FunctionGradient, MnAlgebraicVector>(FunctionGradient(grd, g2, gstep), dgrd);
145 }
146 
147  } // namespace Minuit2
148 
149 } // namespace ROOT
double par[1]
Definition: unuranDistr.cxx:38
Namespace for new ROOT classes and functions.
Definition: ROOT.py:1
#define assert(cond)
Definition: unittest.h:542
bool SyncVector(ROOT::Minuit2::MnAlgebraicVector &mnvector)
Definition: MPIProcess.cxx:160
double GradientStepTolerance() const
Definition: MnStrategy.h:42
determines the relative floating point arithmetic precision.
std::pair< FunctionGradient, MnAlgebraicVector > DeltaGradient(const MinimumParameters &, const FunctionGradient &) const
double sqrt(double)
Double_t x[n]
Definition: legend1.C:17
unsigned int StartElementIndex() const
Definition: MPIProcess.h:56
virtual FunctionGradient operator()(const MinimumParameters &) const
const MnAlgebraicVector & Vec() const
double GradientTolerance() const
Definition: MnStrategy.h:43
const MnMachinePrecision & Precision() const
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
double Eps2() const
eps2 returns 2*sqrt(eps)
Class to calculate an initial estimate of the gradient.
const MnAlgebraicVector & Gstep() const
double Eps() const
eps returns the smallest possible number so that 1.+eps > 1.
double Up() const
Definition: MnFcn.cxx:35
const MnAlgebraicVector & Grad() const
unsigned int EndElementIndex() const
Definition: MPIProcess.h:60
static Vc_ALWAYS_INLINE int_v max(const int_v &x, const int_v &y)
Definition: vector.h:440
const MnMachinePrecision & Precision() const
forwarded interface
unsigned int HessianGradientNCycles() const
Definition: MnStrategy.h:48
const MnAlgebraicVector & G2() const
const Int_t n
Definition: legend1.C:16