Logo 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
20#include "./MPIProcess.h"
21
22#include <cmath>
23#include <cassert>
24
25namespace ROOT {
26
27namespace Minuit2 {
28
30{
31 // use initial gradient as starting point
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
49 return fTransformation.Precision();
50}
51
53{
54 // return number of calculation cycles (defined in strategy)
55 return Strategy().HessianGradientNCycles();
56}
57
58std::pair<FunctionGradient, MnAlgebraicVector>
60{
61 // calculate gradient for Hessian
62 assert(par.IsValid());
63
64 MnPrint print("HessianGradientCalculator");
65
66 MnAlgebraicVector x = par.Vec();
67 MnAlgebraicVector grd = Gradient.Grad();
68 const MnAlgebraicVector &g2 = Gradient.G2();
69 // const MnAlgebraicVector& gstep = Gradient.Gstep();
70 // update also gradient step sizes
71 MnAlgebraicVector gstep = Gradient.Gstep();
72
73 double fcnmin = par.Fval();
74 // std::cout<<"fval: "<<fcnmin<<std::endl;
75
76 double dfmin = 4. * Precision().Eps2() * (fabs(fcnmin) + Fcn().Up());
77
78 unsigned int n = x.size();
80
82 // initial starting values
83 unsigned int startElementIndex = mpiproc.StartElementIndex();
84 unsigned int endElementIndex = mpiproc.EndElementIndex();
85
87
88 for (unsigned int i = startElementIndex; i < endElementIndex; i++) {
89 double xtf = x(i);
90 double dmin = 4. * Precision().Eps2() * (xtf + Precision().Eps2());
91 double epspri = Precision().Eps2() + fabs(grd(i) * Precision().Eps2());
92 double optstp = sqrt(dfmin / (fabs(g2(i)) + epspri));
93 double d = 0.2 * fabs(gstep(i));
94 if (d > optstp)
95 d = optstp;
96 if (d < dmin)
97 d = dmin;
98 double chgold = 10000.;
99 double dgmin = 0.;
100 double grdold = 0.;
101 double grdnew = 0.;
102 for (unsigned int j = 0; j < Ncycle(); j++) {
103 x(i) = xtf + d;
104 double fs1 = fcnCaller(x);
105 x(i) = xtf - d;
106 double fs2 = fcnCaller(x);
107 x(i) = xtf;
108 // double sag = 0.5*(fs1+fs2-2.*fcnmin);
109 // LM: should I calculate also here second derivatives ???
110
111 grdold = grd(i);
112 grdnew = (fs1 - fs2) / (2. * d);
113 dgmin = Precision().Eps() * (std::fabs(fs1) + std::fabs(fs2)) / d;
114 // if(fabs(grdnew) < Precision().Eps()) break;
115 if (grdnew == 0)
116 break;
117 double change = fabs((grdold - grdnew) / grdnew);
118 if (change > chgold && j > 1)
119 break;
120 chgold = change;
121 grd(i) = grdnew;
122 // LM : update also the step sizes
123 gstep(i) = d;
124
125 if (change < 0.05)
126 break;
127 if (fabs(grdold - grdnew) < dgmin)
128 break;
129 if (d < dmin)
130 break;
131 d *= 0.2;
132 }
133
134 dgrd(i) = std::max(dgmin, std::fabs(grdold - grdnew));
135
136 print.Debug("HGC Param :", i, "\t new g1 =", grd(i), "gstep =", d, "dgrd =", dgrd(i));
137 }
138
139 mpiproc.SyncVector(grd);
140 mpiproc.SyncVector(gstep);
141 mpiproc.SyncVector(dgrd);
142
143 return std::pair<FunctionGradient, MnAlgebraicVector>(FunctionGradient(grd, g2, gstep), dgrd);
144}
145
146} // namespace Minuit2
147
148} // namespace ROOT
#define d(i)
Definition RSha256.hxx:102
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
FunctionGradient operator()(const MinimumParameters &) const override
std::pair< FunctionGradient, MnAlgebraicVector > DeltaGradient(const MinimumParameters &, const FunctionGradient &) const
const MnAlgebraicVector & Vec() const
Sets the relative floating point (double) arithmetic precision.
void Debug(const Ts &... args)
Definition MnPrint.h:135
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
FunctionGradient calculateInitialGradient(const MinimumParameters &, const MnUserTransformation &, double errorDef)
Initial rough estimate of the gradient using the parameter step size.
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...