ROOT   Reference Guide
Loading...
Searching...
No Matches
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#include "Minuit2/MnPrint.h"
19#include "Minuit2/MPIProcess.h"
20
21#include <cmath>
22#include <cassert>
23
24namespace ROOT {
25
26namespace Minuit2 {
27
29{
30 // use initial gradient as starting point
32 FunctionGradient gra = gc(par);
33
34 return (*this)(par, gra);
35}
36
38operator()(const MinimumParameters &par, const FunctionGradient &Gradient) const
39{
40 // interface of the base class. Use DeltaGradient for op.
41 std::pair<FunctionGradient, MnAlgebraicVector> mypair = DeltaGradient(par, Gradient);
42
43 return mypair.first;
44}
45
47{
48 // return the precision
50}
51
53{
54 // return number of calculation cycles (defined in strategy)
56}
57
59{
60 // return tolerance on step size (defined in strategy)
62}
63
65{
66 // return gradient tolerance (defines in strategy)
67 return Strategy().GradientTolerance();
68}
69
70std::pair<FunctionGradient, MnAlgebraicVector>
72{
73 // calculate gradient for Hessian
74 assert(par.IsValid());
75
76 MnPrint print("HessianGradientCalculator");
77
78 MnAlgebraicVector x = par.Vec();
79 MnAlgebraicVector grd = Gradient.Grad();
80 const MnAlgebraicVector &g2 = Gradient.G2();
81 // const MnAlgebraicVector& gstep = Gradient.Gstep();
82 // update also gradient step sizes
83 MnAlgebraicVector gstep = Gradient.Gstep();
84
85 double fcnmin = par.Fval();
86 // std::cout<<"fval: "<<fcnmin<<std::endl;
87
88 double dfmin = 4. * Precision().Eps2() * (fabs(fcnmin) + Fcn().Up());
89
90 unsigned int n = x.size();
91 MnAlgebraicVector dgrd(n);
92
93 MPIProcess mpiproc(n, 0);
94 // initial starting values
95 unsigned int startElementIndex = mpiproc.StartElementIndex();
96 unsigned int endElementIndex = mpiproc.EndElementIndex();
97
98 for (unsigned int i = startElementIndex; i < endElementIndex; i++) {
99 double xtf = x(i);
100 double dmin = 4. * Precision().Eps2() * (xtf + Precision().Eps2());
101 double epspri = Precision().Eps2() + fabs(grd(i) * Precision().Eps2());
102 double optstp = sqrt(dfmin / (fabs(g2(i)) + epspri));
103 double d = 0.2 * fabs(gstep(i));
104 if (d > optstp)
105 d = optstp;
106 if (d < dmin)
107 d = dmin;
108 double chgold = 10000.;
109 double dgmin = 0.;
110 double grdold = 0.;
111 double grdnew = 0.;
112 for (unsigned int j = 0; j < Ncycle(); j++) {
113 x(i) = xtf + d;
114 double fs1 = Fcn()(x);
115 x(i) = xtf - d;
116 double fs2 = Fcn()(x);
117 x(i) = xtf;
118 // double sag = 0.5*(fs1+fs2-2.*fcnmin);
119 // LM: should I calculate also here second derivatives ???
120
121 grdold = grd(i);
122 grdnew = (fs1 - fs2) / (2. * d);
123 dgmin = Precision().Eps() * (std::fabs(fs1) + std::fabs(fs2)) / d;
124 // if(fabs(grdnew) < Precision().Eps()) break;
125 if (grdnew == 0)
126 break;
127 double change = fabs((grdold - grdnew) / grdnew);
128 if (change > chgold && j > 1)
129 break;
130 chgold = change;
131 grd(i) = grdnew;
132 // LM : update also the step sizes
133 gstep(i) = d;
134
135 if (change < 0.05)
136 break;
137 if (fabs(grdold - grdnew) < dgmin)
138 break;
139 if (d < dmin)
140 break;
141 d *= 0.2;
142 }
143
144 dgrd(i) = std::max(dgmin, std::fabs(grdold - grdnew));
145
146 print.Debug("HGC Param :", i, "\t new g1 =", grd(i), "gstep =", d, "dgrd =", dgrd(i));
147 }
148
149 mpiproc.SyncVector(grd);
150 mpiproc.SyncVector(gstep);
151 mpiproc.SyncVector(dgrd);
152
153 return std::pair<FunctionGradient, MnAlgebraicVector>(FunctionGradient(grd, g2, gstep), dgrd);
154}
155
156} // namespace Minuit2
157
158} // 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
Class to calculate an initial estimate of the gradient.
bool SyncVector(ROOT::Minuit2::MnAlgebraicVector &mnvector)
unsigned int StartElementIndex() const
Definition MPIProcess.h:55
unsigned int EndElementIndex() const
Definition MPIProcess.h:61
const MnAlgebraicVector & Vec() const
double Up() const
Definition MnFcn.cxx:39
Sets the relative floating point (double) arithmetic precision.
double Eps() const
eps returns the smallest possible number so that 1.+eps > 1.
double Eps2() const
eps2 returns 2*sqrt(eps)
void Debug(const Ts &... args)
Definition MnPrint.h:138
unsigned int HessianGradientNCycles() const
Definition MnStrategy.h:47
double GradientStepTolerance() const
Definition MnStrategy.h:41
double GradientTolerance() const
Definition MnStrategy.h:42
const MnMachinePrecision & Precision() const
forwarded interface
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...