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