Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
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 
10 #include "Minuit2/Numerical2PGradientCalculator.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 
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 
36 namespace ROOT {
37 
38  namespace Minuit2 {
39 
40 
41 FunctionGradient Numerical2PGradientCalculator::operator()(const MinimumParameters& par) const {
42  // calculate gradient using Initial gradient calculator and from MinimumParameters object
43 
44  InitialGradientCalculator gc(fFcn, fTransformation, fStrategy);
45  FunctionGradient gra = gc(par);
46 
47  return (*this)(par, gra);
48 }
49 
50 
51 // comment it, because it was added
52 FunctionGradient 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 
72 FunctionGradient Numerical2PGradientCalculator::operator()(const MinimumParameters& par, const FunctionGradient& Gradient) const {
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 
237 const MnMachinePrecision& Numerical2PGradientCalculator::Precision() const {
238  // return global precision (set in transformation)
239  return fTransformation.Precision();
240 }
241 
242 unsigned int Numerical2PGradientCalculator::Ncycle() const {
243  // return number of cycles for gradient calculation (set in strategy object)
244  return Strategy().GradientNCycles();
245 }
246 
247 double Numerical2PGradientCalculator::StepTolerance() const {
248  // return gradient step tolerance (set in strategy object)
249  return Strategy().GradientStepTolerance();
250 }
251 
252 double Numerical2PGradientCalculator::GradTolerance() const {
253  // return gradient tolerance (set in strategy object)
254  return Strategy().GradientTolerance();
255 }
256 
257 
258  } // namespace Minuit2
259 
260 } // namespace ROOT