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 
30 namespace ROOT {
31 
32 namespace 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
45 FunctionGradient 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 
64 operator()(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 print("Numerical2PGradientCalculator[OpenMP]");
179 #endif
180  if (i == 0 && j == 0) {
181  print.Debug([&](std::ostream &os) {
182  os << std::setw(10) << "parameter" << std::setw(6) << "cycle" << std::setw(15) << "x" << std::setw(15)
183  << "step" << std::setw(15) << "f1" << std::setw(15) << "f2" << std::setw(15) << "grd"
184  << std::setw(15) << "g2" << std::endl;
185  });
186  }
187  print.Debug([&](std::ostream &os) {
188  const int pr = os.precision(13);
189  const int iext = Trafo().ExtOfInt(i);
190  os << std::setw(10) << Trafo().Name(iext) << std::setw(5) << j << " " << x(i) << " " << step << " "
191  << fs1 << " " << fs2 << " " << grd(i) << " " << g2(i) << std::endl;
192  os.precision(pr);
193  });
194  }
195 
196  if (std::fabs(grdb4 - grd(i)) / (std::fabs(grd(i)) + dfmin / step) < GradTolerance()) {
197  // std::cout<<"j= "<<j<<std::endl;
198  // std::cout<<"step= "<<step<<std::endl;
199  // std::cout<<"fs1, fs2: "<<fs1<<" "<<fs2<<std::endl;
200  // std::cout<<"fs1-fs2: "<<fs1-fs2<<std::endl;
201  break;
202  }
203  }
204 
205  // vgrd(i) = grd;
206  // vgrd2(i) = g2;
207  // vgstp(i) = gstep;
208  }
209 
210 #ifndef _OPENMP
211  mpiproc.SyncVector(grd);
212  mpiproc.SyncVector(g2);
213  mpiproc.SyncVector(gstep);
214 #endif
215 
216  // print after parallel processing to avoid synchronization issues
217  print.Debug([&](std::ostream &os) {
218  const int pr = os.precision(13);
219  os << std::endl;
220  os << std::setw(14) << "Parameter" << std::setw(14) << "Gradient" << std::setw(14) << "g2 " << std::setw(14)
221  << "step" << std::endl;
222  for (int i = 0; i < int(n); i++) {
223  const int iext = Trafo().ExtOfInt(i);
224  os << std::setw(14) << Trafo().Name(iext) << " " << grd(i) << " " << g2(i) << " " << gstep(i) << std::endl;
225  }
226  os.precision(pr);
227  });
228 
229  return FunctionGradient(grd, g2, gstep);
230 }
231 
233 {
234  // return global precision (set in transformation)
235  return fTransformation.Precision();
236 }
237 
239 {
240  // return number of cycles for gradient calculation (set in strategy object)
241  return Strategy().GradientNCycles();
242 }
243 
245 {
246  // return gradient step tolerance (set in strategy object)
247  return Strategy().GradientStepTolerance();
248 }
249 
251 {
252  // return gradient tolerance (set in strategy object)
253  return Strategy().GradientTolerance();
254 }
255 
256 } // namespace Minuit2
257 
258 } // namespace ROOT
ROOT::Minuit2::MnMachinePrecision::Eps2
double Eps2() const
eps2 returns 2*sqrt(eps)
Definition: MnMachinePrecision.h:41
n
const Int_t n
Definition: legend1.C:16
ROOT::Minuit2::Numerical2PGradientCalculator::Strategy
const MnStrategy & Strategy() const
Definition: Numerical2PGradientCalculator.h:51
ROOT::Minuit2::Numerical2PGradientCalculator::Trafo
const MnUserTransformation & Trafo() const
Definition: Numerical2PGradientCalculator.h:49
ROOT::Minuit2::MPIProcess
Definition: MPIProcess.h:42
MnMachinePrecision.h
ROOT::Minuit2::LAVector
Definition: LAVector.h:32
ROOT::Minuit2::MnPrint::Debug
void Debug(const Ts &... args)
Definition: MnPrint.h:138
MnStrategy.h
ROOT::Minuit2::MnStrategy::GradientNCycles
unsigned int GradientNCycles() const
Definition: MnStrategy.h:40
ROOT::Minuit2::Numerical2PGradientCalculator::StepTolerance
double StepTolerance() const
Definition: Numerical2PGradientCalculator.cxx:244
MPIProcess.h
ROOT::Minuit2::MnUserTransformation::Precision
const MnMachinePrecision & Precision() const
forwarded interface
Definition: MnUserTransformation.h:117
ROOT::Minuit2::MPIProcess::StartElementIndex
unsigned int StartElementIndex() const
Definition: MPIProcess.h:55
MnUserTransformation.h
ROOT::Minuit2::Numerical2PGradientCalculator::Fcn
const MnFcn & Fcn() const
Definition: Numerical2PGradientCalculator.h:48
ROOT::Minuit2::MinimumParameters::Vec
const MnAlgebraicVector & Vec() const
Definition: MinimumParameters.h:38
ROOT::Minuit2::Numerical2PGradientCalculator::fFcn
const MnFcn & fFcn
Definition: Numerical2PGradientCalculator.h:58
ROOT::Minuit2::MnStrategy::GradientTolerance
double GradientTolerance() const
Definition: MnStrategy.h:42
x
Double_t x[n]
Definition: legend1.C:17
ROOT::Minuit2::InitialGradientCalculator
Class to calculate an initial estimate of the gradient.
Definition: InitialGradientCalculator.h:27
ROOT::Minuit2::FunctionGradient
Definition: FunctionGradient.h:21
ROOT::Minuit2::MnStrategy::GradientStepTolerance
double GradientStepTolerance() const
Definition: MnStrategy.h:41
MnFcn.h
ROOT::Minuit2::FunctionGradient::G2
const MnAlgebraicVector & G2() const
Definition: FunctionGradient.h:45
ROOT::Minuit2::MinimumParameters::Fval
double Fval() const
Definition: MinimumParameters.h:40
ROOT::Minuit2::MnUserTransformation::Name
const char * Name(unsigned int) const
Definition: MnUserTransformation.cxx:493
ROOT::Minuit2::MnMachinePrecision::Eps
double Eps() const
eps returns the smallest possible number so that 1.+eps > 1.
Definition: MnMachinePrecision.h:38
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:21
FunctionGradient.h
ROOT::Minuit2::MnFcn::Up
double Up() const
Definition: MnFcn.cxx:39
ROOT::Minuit2::MPIProcess::SyncVector
bool SyncVector(ROOT::Minuit2::MnAlgebraicVector &mnvector)
Definition: MPIProcess.cxx:153
size
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
ROOT::Minuit2::Numerical2PGradientCalculator::fStrategy
const MnStrategy & fStrategy
Definition: Numerical2PGradientCalculator.h:60
ROOT::Minuit2::Numerical2PGradientCalculator::operator()
virtual FunctionGradient operator()(const MinimumParameters &) const
Definition: Numerical2PGradientCalculator.cxx:34
Parameter
CPyCppyy::Parameter Parameter
Definition: clingwrapper.cxx:48
ROOT::Minuit2::MnMachinePrecision
Sets the relative floating point (double) arithmetic precision.
Definition: MnMachinePrecision.h:32
sqrt
double sqrt(double)
ROOT::Minuit2::FunctionGradient::Grad
const MnAlgebraicVector & Grad() const
Definition: FunctionGradient.h:40
ROOT::Minuit2::MinimumParameters::IsValid
bool IsValid() const
Definition: MinimumParameters.h:41
ROOT::Minuit2::Numerical2PGradientCalculator::Ncycle
unsigned int Ncycle() const
Definition: Numerical2PGradientCalculator.cxx:238
ROOT::Minuit2::FunctionGradient::Gstep
const MnAlgebraicVector & Gstep() const
Definition: FunctionGradient.h:46
ROOT::Minuit2::Numerical2PGradientCalculator::Precision
const MnMachinePrecision & Precision() const
Definition: Numerical2PGradientCalculator.cxx:232
ROOT::Minuit2::Numerical2PGradientCalculator::GradTolerance
double GradTolerance() const
Definition: Numerical2PGradientCalculator.cxx:250
InitialGradientCalculator.h
ROOT::Minuit2::MPIProcess::EndElementIndex
unsigned int EndElementIndex() const
Definition: MPIProcess.h:61
MinimumParameters.h
ROOT::Minuit2::Numerical2PGradientCalculator::fTransformation
const MnUserTransformation & fTransformation
Definition: Numerical2PGradientCalculator.h:59
MnPrint.h
ROOT::Minuit2::MnUserTransformation::ExtOfInt
unsigned int ExtOfInt(unsigned int internal) const
Definition: MnUserTransformation.h:102
ROOT
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
Definition: EExecutionPolicy.hxx:4
ROOT::Minuit2::MnPrint
Definition: MnPrint.h:73
int