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 
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 
248  // return gradient step tolerance (set in strategy object)
249  return Strategy().GradientStepTolerance();
250 }
251 
253  // return gradient tolerance (set in strategy object)
254  return Strategy().GradientTolerance();
255 }
256 
257 
258  } // namespace Minuit2
259 
260 } // namespace ROOT
ROOT::Minuit2::MnMachinePrecision::Eps2
double Eps2() const
eps2 returns 2*sqrt(eps)
Definition: MnMachinePrecision.h:56
n
const Int_t n
Definition: legend1.C:16
ROOT::Minuit2::Numerical2PGradientCalculator::Strategy
const MnStrategy & Strategy() const
Definition: Numerical2PGradientCalculator.h:68
ROOT::Minuit2::Numerical2PGradientCalculator::Trafo
const MnUserTransformation & Trafo() const
Definition: Numerical2PGradientCalculator.h:66
ROOT::Minuit2::MPIProcess
Definition: MPIProcess.h:58
MnMachinePrecision.h
ROOT::Minuit2::LAVector
Definition: LAVector.h:41
MnStrategy.h
ROOT::Minuit2::MnStrategy::GradientNCycles
unsigned int GradientNCycles() const
Definition: MnStrategy.h:53
ROOT::Minuit2::Numerical2PGradientCalculator::StepTolerance
double StepTolerance() const
Definition: Numerical2PGradientCalculator.cxx:255
MPIProcess.h
ROOT::Minuit2::MnUserTransformation::Precision
const MnMachinePrecision & Precision() const
forwarded interface
Definition: MnUserTransformation.h:131
MnUserTransformation.h
ROOT::Minuit2::Numerical2PGradientCalculator::Fcn
const MnFcn & Fcn() const
Definition: Numerical2PGradientCalculator.h:65
ROOT::Minuit2::Numerical2PGradientCalculator::fFcn
const MnFcn & fFcn
Definition: Numerical2PGradientCalculator.h:76
ROOT::Minuit2::MnStrategy::GradientTolerance
double GradientTolerance() const
Definition: MnStrategy.h:55
x
Double_t x[n]
Definition: legend1.C:17
ROOT::Minuit2::FunctionGradient
Definition: FunctionGradient.h:29
ROOT::Minuit2::MnStrategy::GradientStepTolerance
double GradientStepTolerance() const
Definition: MnStrategy.h:54
MnFcn.h
ROOT::Minuit2::MnUserTransformation::Name
const char * Name(unsigned int) const
Definition: MnUserTransformation.cxx:458
ROOT::Minuit2::MnMachinePrecision::Eps
double Eps() const
eps returns the smallest possible number so that 1.+eps > 1.
Definition: MnMachinePrecision.h:53
Numerical2PGradientCalculator.h
ROOT::Math::fabs
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
Definition: UnaryOperators.h:131
ROOT::Minuit2::MinimumParameters
Definition: MinimumParameters.h:29
FunctionGradient.h
ROOT::Minuit2::MnFcn::Up
double Up() const
Definition: MnFcn.cxx:43
ROOT::Minuit2::Numerical2PGradientCalculator::fStrategy
const MnStrategy & fStrategy
Definition: Numerical2PGradientCalculator.h:78
ROOT::Minuit2::Numerical2PGradientCalculator::operator()
virtual FunctionGradient operator()(const MinimumParameters &) const
Definition: Numerical2PGradientCalculator.cxx:49
Parameter
CPyCppyy::Parameter Parameter
Definition: clingwrapper.cxx:48
sqrt
double sqrt(double)
ROOT::Minuit2::Numerical2PGradientCalculator::Ncycle
unsigned int Ncycle() const
Definition: Numerical2PGradientCalculator.cxx:250
ROOT::Minuit2::Numerical2PGradientCalculator::Precision
const MnMachinePrecision & Precision() const
Definition: Numerical2PGradientCalculator.cxx:245
ROOT::Minuit2::Numerical2PGradientCalculator::GradTolerance
double GradTolerance() const
Definition: Numerical2PGradientCalculator.cxx:260
InitialGradientCalculator.h
MinimumParameters.h
ROOT::Minuit2::Numerical2PGradientCalculator::fTransformation
const MnUserTransformation & fTransformation
Definition: Numerical2PGradientCalculator.h:77
MnPrint.h
ROOT::Minuit2::MnUserTransformation::ExtOfInt
unsigned int ExtOfInt(unsigned int internal) const
Definition: MnUserTransformation.h:112
ROOT
VSD Structures.
Definition: StringConv.hxx:21
int