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
10 #include "Minuit2/HessianGradientCalculator.h"
11 #include "Minuit2/InitialGradientCalculator.h"
12 #include "Minuit2/MnFcn.h"
13 #include "Minuit2/MnUserTransformation.h"
14 #include "Minuit2/MnMachinePrecision.h"
15 #include "Minuit2/MinimumParameters.h"
16 #include "Minuit2/FunctionGradient.h"
17 #include "Minuit2/MnStrategy.h"
18 #include "Minuit2/MnPrint.h"
19 #include "Minuit2/MPIProcess.h"
20
21 #include <cmath>
22 #include <cassert>
23
24 namespace ROOT {
25
26 namespace Minuit2 {
27
operator ()(const MinimumParameters & par) const28 FunctionGradient HessianGradientCalculator::operator()(const MinimumParameters &par) const
29 {
30 // use initial gradient as starting point
31 InitialGradientCalculator gc(fFcn, fTransformation, fStrategy);
32 FunctionGradient gra = gc(par);
33
34 return (*this)(par, gra);
35 }
36
37 FunctionGradient HessianGradientCalculator::
operator ()(const MinimumParameters & par,const FunctionGradient & Gradient) const38 operator()(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
Precision() const46 const MnMachinePrecision &HessianGradientCalculator::Precision() const
47 {
48 // return the precision
49 return fTransformation.Precision();
50 }
51
Ncycle() const52 unsigned int HessianGradientCalculator::Ncycle() const
53 {
54 // return number of calculation cycles (defined in strategy)
55 return Strategy().HessianGradientNCycles();
56 }
57
StepTolerance() const58 double HessianGradientCalculator::StepTolerance() const
59 {
60 // return tolerance on step size (defined in strategy)
61 return Strategy().GradientStepTolerance();
62 }
63
GradTolerance() const64 double HessianGradientCalculator::GradTolerance() const
65 {
66 // return gradient tolerance (defines in strategy)
67 return Strategy().GradientTolerance();
68 }
69
70 std::pair<FunctionGradient, MnAlgebraicVector>
DeltaGradient(const MinimumParameters & par,const FunctionGradient & Gradient) const71 HessianGradientCalculator::DeltaGradient(const MinimumParameters &par, const FunctionGradient &Gradient) const
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
159