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
29namespace 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
52}
53
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
69std::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
#define d(i)
Definition: RSha256.hxx:102
double sqrt(double)
const MnAlgebraicVector & Gstep() const
const MnAlgebraicVector & Grad() const
const MnAlgebraicVector & G2() const
std::pair< FunctionGradient, MnAlgebraicVector > DeltaGradient(const MinimumParameters &, const FunctionGradient &) const
virtual FunctionGradient operator()(const MinimumParameters &) const
const MnMachinePrecision & Precision() const
Class to calculate an initial estimate of the gradient.
bool SyncVector(ROOT::Minuit2::MnAlgebraicVector &mnvector)
Definition: MPIProcess.cxx:160
unsigned int StartElementIndex() const
Definition: MPIProcess.h:57
unsigned int EndElementIndex() const
Definition: MPIProcess.h:61
const MnAlgebraicVector & Vec() const
double Up() const
Definition: MnFcn.cxx:35
determines the relative floating point arithmetic precision.
double Eps() const
eps returns the smallest possible number so that 1.+eps > 1.
double Eps2() const
eps2 returns 2*sqrt(eps)
unsigned int HessianGradientNCycles() const
Definition: MnStrategy.h:48
double GradientStepTolerance() const
Definition: MnStrategy.h:42
double GradientTolerance() const
Definition: MnStrategy.h:43
const MnMachinePrecision & Precision() const
forwarded interface
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
VSD Structures.
Definition: StringConv.hxx:21