Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
Numerical2PGradientCalculator.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#ifdef _OPENMP
23#include <omp.h>
24#endif
25
26#include <cmath>
27#include <cassert>
28#include <iomanip>
29
30namespace ROOT {
31
32namespace Minuit2 {
33
35{
36 // calculate gradient using Initial gradient calculator and from MinimumParameters object
37
39
40 return (*this)(par, gra);
41}
42
44operator()(const MinimumParameters &par, const FunctionGradient &Gradient) const
45{
46 // calculate numerical gradient from MinimumParameters object
47 // the algorithm takes correctly care when the gradient is approximately zero
48
50
51 // std::cout<<"########### Numerical2PDerivative"<<std::endl;
52 // std::cout<<"initial grd: "<<Gradient.Grad()<<std::endl;
53 // std::cout<<"position: "<<par.Vec()<<std::endl;
54 MnPrint print("Numerical2PGradientCalculator");
55
56 assert(par.IsValid());
57
58 double fcnmin = par.Fval();
59 // std::cout<<"fval: "<<fcnmin<<std::endl;
60
61 double eps2 = trafo.Precision().Eps2();
62 double eps = trafo.Precision().Eps();
63
64 //print.Trace("Assumed precision eps", eps, "eps2", eps2);
65
66 double dfmin = 8. * eps2 * (std::fabs(fcnmin) + fFcn.Up());
67 double vrysml = 8. * eps * eps;
68 // double vrysml = std::max(1.e-4, eps2);
69 // std::cout<<"dfmin= "<<dfmin<<std::endl;
70 // std::cout<<"vrysml= "<<vrysml<<std::endl;
71 // std::cout << " ncycle " << Ncycle() << std::endl;
72
73 unsigned int n = (par.Vec()).size();
74 unsigned int ncycle = fStrategy.GradientNCycles();
75 // MnAlgebraicVector vgrd(n), vgrd2(n), vgstp(n);
76 MnAlgebraicVector grd = Gradient.Grad();
77 MnAlgebraicVector g2 = Gradient.G2();
78 MnAlgebraicVector gstep = Gradient.Gstep();
79
80 print.Debug("Calculating gradient around function value", fcnmin, "\n\t at point", par.Vec());
81
82#ifndef _OPENMP
83
85
86 // for serial execution this can be outside the loop
87 MnAlgebraicVector x = par.Vec();
89
90 unsigned int startElementIndex = mpiproc.StartElementIndex();
91 unsigned int endElementIndex = mpiproc.EndElementIndex();
92
93 for (unsigned int i = startElementIndex; i < endElementIndex; i++) {
94
95#else
96
97 // parallelize this loop using OpenMP
98//#define N_PARALLEL_PAR 5
99#pragma omp parallel for if (fDoParallelOMP)
100 //#pragma omp for schedule (static, N_PARALLEL_PAR)
101
102 for (unsigned int i = 0; i < n; i++) {
103
104#endif
105
106#ifdef _OPENMP
107 // create in loop since each thread will use its own copy
108 MnAlgebraicVector x = par.Vec();
110#endif
111
112 double xtf = x(i);
113 double epspri = eps2 + std::fabs(grd(i) * eps2);
114 double stepb4 = 0.;
115 for (unsigned int j = 0; j < ncycle; j++) {
116 double optstp = std::sqrt(dfmin / (std::fabs(g2(i)) + epspri));
117 double step = std::max(optstp, std::fabs(0.1 * gstep(i)));
118 // std::cout<<"step: "<<step;
119 if (trafo.Parameter(trafo.ExtOfInt(i)).HasLimits()) {
120 if (step > 0.5)
121 step = 0.5;
122 }
123 double stpmax = 10. * std::fabs(gstep(i));
124 if (step > stpmax)
125 step = stpmax;
126 // std::cout<<" "<<step;
127 double stpmin = std::max(vrysml, 8. * std::fabs(eps2 * x(i)));
128 if (step < stpmin)
129 step = stpmin;
130 // std::cout<<" "<<step<<std::endl;
131 // std::cout<<"step: "<<step<<std::endl;
132 if (std::fabs((step - stepb4) / step) < fStrategy.GradientStepTolerance()) {
133 // std::cout<<"(step-stepb4)/step"<<std::endl;
134 // std::cout<<"j= "<<j<<std::endl;
135 // std::cout<<"step= "<<step<<std::endl;
136 break;
137 }
138 gstep(i) = step;
139 stepb4 = step;
140
141 x(i) = xtf + step;
142 double fs1 = mfcnCaller(x);
143 x(i) = xtf - step;
144 double fs2 = mfcnCaller(x);
145 x(i) = xtf;
146
147 double grdb4 = grd(i);
148 grd(i) = 0.5 * (fs1 - fs2) / step;
149 g2(i) = (fs1 + fs2 - 2. * fcnmin) / step / step;
150
151#ifdef _OPENMP
152#pragma omp critical
153#endif
154 {
155#ifdef _OPENMP
156 // must create thread-local MnPrint instances when printing inside threads
157 MnPrint printtl("Numerical2PGradientCalculator[OpenMP]");
158#endif
159 if (i == 0 && j == 0) {
160#ifdef _OPENMP
161 printtl.Trace([&](std::ostream &os) {
162#else
163 print.Trace([&](std::ostream &os) {
164#endif
165 os << std::setw(10) << "parameter" << std::setw(6) << "cycle" << std::setw(15) << "x" << std::setw(15)
166 << "step" << std::setw(15) << "f1" << std::setw(15) << "f2" << std::setw(15) << "grd"
167 << std::setw(15) << "g2" << std::endl;
168 });
169 }
170#ifdef _OPENMP
171 printtl.Trace([&](std::ostream &os) {
172#else
173 print.Trace([&](std::ostream &os) {
174#endif
175 const int pr = os.precision(13);
176 const int iext = trafo.ExtOfInt(i);
177 os << std::setw(10) << trafo.Name(iext) << std::setw(5) << j << " " << x(i) << " " << step << " "
178 << fs1 << " " << fs2 << " " << grd(i) << " " << g2(i) << std::endl;
179 os.precision(pr);
180 });
181 }
182
183 if (std::fabs(grdb4 - grd(i)) / (std::fabs(grd(i)) + dfmin / step) < fStrategy.GradientTolerance()) {
184 // std::cout<<"j= "<<j<<std::endl;
185 // std::cout<<"step= "<<step<<std::endl;
186 // std::cout<<"fs1, fs2: "<<fs1<<" "<<fs2<<std::endl;
187 // std::cout<<"fs1-fs2: "<<fs1-fs2<<std::endl;
188 break;
189 }
190 }
191
192 // vgrd(i) = grd;
193 // vgrd2(i) = g2;
194 // vgstp(i) = gstep;
195 }
196
197#ifndef _OPENMP
198 mpiproc.SyncVector(grd);
199 mpiproc.SyncVector(g2);
200 mpiproc.SyncVector(gstep);
201#endif
202
203 // print after parallel processing to avoid synchronization issues
204 print.Debug([&](std::ostream &os) {
205 const int pr = os.precision(13);
206 os << std::endl;
207 os << std::setw(14) << "Parameter" << std::setw(14) << "Gradient" << std::setw(14) << "g2 " << std::setw(14)
208 << "step" << std::endl;
209 for (int i = 0; i < int(n); i++) {
210 const int iext = trafo.ExtOfInt(i);
211 os << std::setw(14) << trafo.Name(iext) << " " << grd(i) << " " << g2(i) << " " << gstep(i) << std::endl;
212 }
213 os.precision(pr);
214 });
215
216 return FunctionGradient(grd, g2, gstep);
217}
218
219} // namespace Minuit2
220
221} // namespace ROOT
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
const MnAlgebraicVector & Vec() const
void Debug(const Ts &... args)
Definition MnPrint.h:135
void Trace(const Ts &... args)
Definition MnPrint.h:141
class dealing with the transformation between user specified parameters (external) and internal param...
FunctionGradient operator()(const MinimumParameters &) const override
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.
Namespace for new ROOT classes and functions.