Logo ROOT  
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 
34 FunctionGradient HessianGradientCalculator::operator()(const MinimumParameters& par) const {
35  // use initial gradient as starting point
36  InitialGradientCalculator gc(fFcn, fTransformation, fStrategy);
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
ROOT::Minuit2::MnMachinePrecision::Eps2
double Eps2() const
eps2 returns 2*sqrt(eps)
Definition: MnMachinePrecision.h:56
n
const Int_t n
Definition: legend1.C:16
ROOT::Minuit2::MPIProcess
Definition: MPIProcess.h:58
MnMachinePrecision.h
ROOT::Minuit2::LAVector
Definition: LAVector.h:41
MnStrategy.h
MPIProcess.h
ROOT::Minuit2::MnUserTransformation::Precision
const MnMachinePrecision & Precision() const
forwarded interface
Definition: MnUserTransformation.h:131
MnUserTransformation.h
ROOT::Minuit2::HessianGradientCalculator::Precision
const MnMachinePrecision & Precision() const
Definition: HessianGradientCalculator.cxx:57
ROOT::Minuit2::MinimumParameters::Vec
const MnAlgebraicVector & Vec() const
Definition: MinimumParameters.h:57
ROOT::Minuit2::MnStrategy::GradientTolerance
double GradientTolerance() const
Definition: MnStrategy.h:55
x
Double_t x[n]
Definition: legend1.C:17
ROOT::Minuit2::HessianGradientCalculator::Strategy
const MnStrategy & Strategy() const
Definition: HessianGradientCalculator.h:60
ROOT::Minuit2::FunctionGradient
Definition: FunctionGradient.h:29
ROOT::Minuit2::MnStrategy::GradientStepTolerance
double GradientStepTolerance() const
Definition: MnStrategy.h:54
ROOT::Minuit2::HessianGradientCalculator::fTransformation
const MnUserTransformation & fTransformation
Definition: HessianGradientCalculator.h:69
MnFcn.h
ROOT::Minuit2::FunctionGradient::G2
const MnAlgebraicVector & G2() const
Definition: FunctionGradient.h:63
ROOT::Minuit2::MinimumParameters::Fval
double Fval() const
Definition: MinimumParameters.h:59
ROOT::Minuit2::MnMachinePrecision::Eps
double Eps() const
eps returns the smallest possible number so that 1.+eps > 1.
Definition: MnMachinePrecision.h:53
ROOT::Minuit2::HessianGradientCalculator::operator()
virtual FunctionGradient operator()(const MinimumParameters &) const
Definition: HessianGradientCalculator.cxx:42
ROOT::Minuit2::HessianGradientCalculator::StepTolerance
double StepTolerance() const
Definition: HessianGradientCalculator.cxx:67
ROOT::Math::fabs
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
Definition: UnaryOperators.h:131
ROOT::Minuit2::MinimumParameters
Definition: MinimumParameters.h:29
FunctionGradient.h
ROOT::Minuit2::MnFcn::Up
double Up() const
Definition: MnFcn.cxx:43
ROOT::Minuit2::HessianGradientCalculator::Fcn
const MnFcn & Fcn() const
Definition: HessianGradientCalculator.h:57
ROOT::Minuit2::MnMachinePrecision
Sets the relative floating point (double) arithmetic precision.
Definition: MnMachinePrecision.h:41
ROOT::Minuit2::HessianGradientCalculator::fStrategy
const MnStrategy & fStrategy
Definition: HessianGradientCalculator.h:70
sqrt
double sqrt(double)
ROOT::Minuit2::HessianGradientCalculator::fFcn
const MnFcn & fFcn
Definition: HessianGradientCalculator.h:68
ROOT::Minuit2::FunctionGradient::Grad
const MnAlgebraicVector & Grad() const
Definition: FunctionGradient.h:58
ROOT::Minuit2::MinimumParameters::IsValid
bool IsValid() const
Definition: MinimumParameters.h:60
ROOT::Minuit2::FunctionGradient::Gstep
const MnAlgebraicVector & Gstep() const
Definition: FunctionGradient.h:64
HessianGradientCalculator.h
InitialGradientCalculator.h
ROOT::Minuit2::HessianGradientCalculator::Ncycle
unsigned int Ncycle() const
Definition: HessianGradientCalculator.cxx:62
ROOT::Minuit2::MnStrategy::HessianGradientNCycles
unsigned int HessianGradientNCycles() const
Definition: MnStrategy.h:60
d
#define d(i)
Definition: RSha256.hxx:120
ROOT::Minuit2::HessianGradientCalculator::DeltaGradient
std::pair< FunctionGradient, MnAlgebraicVector > DeltaGradient(const MinimumParameters &, const FunctionGradient &) const
Definition: HessianGradientCalculator.cxx:77
MinimumParameters.h
ROOT::Minuit2::HessianGradientCalculator::GradTolerance
double GradTolerance() const
Definition: HessianGradientCalculator.cxx:72
MnPrint.h
ROOT
VSD Structures.
Definition: StringConv.hxx:21