Logo ROOT  
Reference Guide
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#ifdef _OPENMP
21#include <omp.h>
22#endif
23
24#include <cmath>
25#include <cassert>
26#include <iomanip>
27
28#include "Minuit2/MPIProcess.h"
29
30namespace ROOT {
31
32namespace Minuit2 {
33
35{
36 // calculate gradient using Initial gradient calculator and from MinimumParameters object
37
39 FunctionGradient gra = gc(par);
40
41 return (*this)(par, gra);
42}
43
44// comment it, because it was added
45FunctionGradient Numerical2PGradientCalculator::operator()(const std::vector<double> &params) const
46{
47 // calculate gradient from an std;:vector of paramteters
48
49 int npar = params.size();
50
51 MnAlgebraicVector par(npar);
52 for (int i = 0; i < npar; ++i) {
53 par(i) = params[i];
54 }
55
56 double fval = Fcn()(par);
57
58 MinimumParameters minpars = MinimumParameters(par, fval);
59
60 return (*this)(minpars);
61}
62
64operator()(const MinimumParameters &par, const FunctionGradient &Gradient) const
65{
66 // calculate numerical gradient from MinimumParameters object
67 // the algorithm takes correctly care when the gradient is approximatly zero
68
69 // std::cout<<"########### Numerical2PDerivative"<<std::endl;
70 // std::cout<<"initial grd: "<<Gradient.Grad()<<std::endl;
71 // std::cout<<"position: "<<par.Vec()<<std::endl;
72 MnPrint print("Numerical2PGradientCalculator");
73
74 assert(par.IsValid());
75
76 double fcnmin = par.Fval();
77 // std::cout<<"fval: "<<fcnmin<<std::endl;
78
79 double eps2 = Precision().Eps2();
80 double eps = Precision().Eps();
81
82 print.Debug("Assumed precision eps", eps, "eps2", eps2);
83
84 double dfmin = 8. * eps2 * (std::fabs(fcnmin) + Fcn().Up());
85 double vrysml = 8. * eps * eps;
86 // double vrysml = std::max(1.e-4, eps2);
87 // std::cout<<"dfmin= "<<dfmin<<std::endl;
88 // std::cout<<"vrysml= "<<vrysml<<std::endl;
89 // std::cout << " ncycle " << Ncycle() << std::endl;
90
91 unsigned int n = (par.Vec()).size();
92 unsigned int ncycle = Ncycle();
93 // MnAlgebraicVector vgrd(n), vgrd2(n), vgstp(n);
94 MnAlgebraicVector grd = Gradient.Grad();
95 MnAlgebraicVector g2 = Gradient.G2();
96 MnAlgebraicVector gstep = Gradient.Gstep();
97
98 print.Debug("Calculating gradient around value", fcnmin, "at point", par.Vec());
99
100#ifndef _OPENMP
101
102 MPIProcess mpiproc(n, 0);
103
104 // for serial execution this can be outside the loop
105 MnAlgebraicVector x = par.Vec();
106
107 unsigned int startElementIndex = mpiproc.StartElementIndex();
108 unsigned int endElementIndex = mpiproc.EndElementIndex();
109
110 for (unsigned int i = startElementIndex; i < endElementIndex; i++) {
111
112#else
113
114 // parallelize this loop using OpenMP
115//#define N_PARALLEL_PAR 5
116#pragma omp parallel
117#pragma omp for
118 //#pragma omp for schedule (static, N_PARALLEL_PAR)
119
120 for (int i = 0; i < int(n); i++) {
121
122#endif
123
124#ifdef _OPENMP
125 // create in loop since each thread will use its own copy
126 MnAlgebraicVector x = par.Vec();
127#endif
128
129 double xtf = x(i);
130 double epspri = eps2 + std::fabs(grd(i) * eps2);
131 double stepb4 = 0.;
132 for (unsigned int j = 0; j < ncycle; j++) {
133 double optstp = std::sqrt(dfmin / (std::fabs(g2(i)) + epspri));
134 double step = std::max(optstp, std::fabs(0.1 * gstep(i)));
135 // std::cout<<"step: "<<step;
136 if (Trafo().Parameter(Trafo().ExtOfInt(i)).HasLimits()) {
137 if (step > 0.5)
138 step = 0.5;
139 }
140 double stpmax = 10. * std::fabs(gstep(i));
141 if (step > stpmax)
142 step = stpmax;
143 // std::cout<<" "<<step;
144 double stpmin = std::max(vrysml, 8. * std::fabs(eps2 * x(i)));
145 if (step < stpmin)
146 step = stpmin;
147 // std::cout<<" "<<step<<std::endl;
148 // std::cout<<"step: "<<step<<std::endl;
149 if (std::fabs((step - stepb4) / step) < StepTolerance()) {
150 // std::cout<<"(step-stepb4)/step"<<std::endl;
151 // std::cout<<"j= "<<j<<std::endl;
152 // std::cout<<"step= "<<step<<std::endl;
153 break;
154 }
155 gstep(i) = step;
156 stepb4 = step;
157 // MnAlgebraicVector pstep(n);
158 // pstep(i) = step;
159 // double fs1 = Fcn()(pstate + pstep);
160 // double fs2 = Fcn()(pstate - pstep);
161
162 x(i) = xtf + step;
163 double fs1 = Fcn()(x);
164 x(i) = xtf - step;
165 double fs2 = Fcn()(x);
166 x(i) = xtf;
167
168 double grdb4 = grd(i);
169 grd(i) = 0.5 * (fs1 - fs2) / step;
170 g2(i) = (fs1 + fs2 - 2. * fcnmin) / step / step;
171
172#ifdef _OPENMP
173#pragma omp critical
174#endif
175 {
176#ifdef _OPENMP
177 // must create thread-local MnPrint instances when printing inside threads
178 MnPrint printtl("Numerical2PGradientCalculator[OpenMP]");
179#endif
180 if (i == 0 && j == 0) {
181#ifdef _OPENMP
182 printtl.Debug([&](std::ostream &os) {
183#else
184 print.Debug([&](std::ostream &os) {
185#endif
186 os << std::setw(10) << "parameter" << std::setw(6) << "cycle" << std::setw(15) << "x" << std::setw(15)
187 << "step" << std::setw(15) << "f1" << std::setw(15) << "f2" << std::setw(15) << "grd"
188 << std::setw(15) << "g2" << std::endl;
189 });
190 }
191#ifdef _OPENMP
192 printtl.Debug([&](std::ostream &os) {
193#else
194 print.Debug([&](std::ostream &os) {
195#endif
196 const int pr = os.precision(13);
197 const int iext = Trafo().ExtOfInt(i);
198 os << std::setw(10) << Trafo().Name(iext) << std::setw(5) << j << " " << x(i) << " " << step << " "
199 << fs1 << " " << fs2 << " " << grd(i) << " " << g2(i) << std::endl;
200 os.precision(pr);
201 });
202 }
203
204 if (std::fabs(grdb4 - grd(i)) / (std::fabs(grd(i)) + dfmin / step) < GradTolerance()) {
205 // std::cout<<"j= "<<j<<std::endl;
206 // std::cout<<"step= "<<step<<std::endl;
207 // std::cout<<"fs1, fs2: "<<fs1<<" "<<fs2<<std::endl;
208 // std::cout<<"fs1-fs2: "<<fs1-fs2<<std::endl;
209 break;
210 }
211 }
212
213 // vgrd(i) = grd;
214 // vgrd2(i) = g2;
215 // vgstp(i) = gstep;
216 }
217
218#ifndef _OPENMP
219 mpiproc.SyncVector(grd);
220 mpiproc.SyncVector(g2);
221 mpiproc.SyncVector(gstep);
222#endif
223
224 // print after parallel processing to avoid synchronization issues
225 print.Debug([&](std::ostream &os) {
226 const int pr = os.precision(13);
227 os << std::endl;
228 os << std::setw(14) << "Parameter" << std::setw(14) << "Gradient" << std::setw(14) << "g2 " << std::setw(14)
229 << "step" << std::endl;
230 for (int i = 0; i < int(n); i++) {
231 const int iext = Trafo().ExtOfInt(i);
232 os << std::setw(14) << Trafo().Name(iext) << " " << grd(i) << " " << g2(i) << " " << gstep(i) << std::endl;
233 }
234 os.precision(pr);
235 });
236
237 return FunctionGradient(grd, g2, gstep);
238}
239
241{
242 // return global precision (set in transformation)
243 return fTransformation.Precision();
244}
245
247{
248 // return number of cycles for gradient calculation (set in strategy object)
249 return Strategy().GradientNCycles();
250}
251
253{
254 // return gradient step tolerance (set in strategy object)
256}
257
259{
260 // return gradient tolerance (set in strategy object)
261 return Strategy().GradientTolerance();
262}
263
264} // namespace Minuit2
265
266} // namespace ROOT
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void gc
const MnAlgebraicVector & Gstep() const
const MnAlgebraicVector & Grad() const
const MnAlgebraicVector & G2() const
Class to calculate an initial estimate of the gradient.
bool SyncVector(ROOT::Minuit2::MnAlgebraicVector &mnvector)
Definition: MPIProcess.cxx:153
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:147
double GradientStepTolerance() const
Definition: MnStrategy.h:41
double GradientTolerance() const
Definition: MnStrategy.h:42
unsigned int GradientNCycles() const
Definition: MnStrategy.h:40
unsigned int ExtOfInt(unsigned int internal) const
const char * Name(unsigned int) const
const MnMachinePrecision & Precision() const
forwarded interface
FunctionGradient operator()(const MinimumParameters &) const override
CPyCppyy::Parameter Parameter
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.