Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
FitUtil.cxx
Go to the documentation of this file.
1 // @(#)root/mathcore:$Id$
2 // Author: L. Moneta Tue Nov 28 10:52:47 2006
3 
4 /**********************************************************************
5  * *
6  * Copyright (c) 2006 LCG ROOT Math Team, CERN/PH-SFT *
7  * *
8  * *
9  **********************************************************************/
10 
11 // Implementation file for class FitUtil
12 
13 #include "Fit/FitUtil.h"
14 
15 #include "Fit/BinData.h"
16 #include "Fit/UnBinData.h"
17 
18 #include "Math/IFunctionfwd.h"
19 #include "Math/IParamFunction.h"
20 #include "Math/Integrator.h"
21 #include "Math/IntegratorMultiDim.h"
22 #include "Math/WrappedFunction.h"
23 #include "Math/OneDimFunctionAdapter.h"
24 #include "Math/RichardsonDerivator.h"
25 
26 #include "Math/Error.h"
27 #include "Math/Util.h" // for safe log(x)
28 
29 #include <limits>
30 #include <cmath>
31 #include <cassert>
32 #include <algorithm>
33 #include <numeric>
34 //#include <memory>
35 
36 #include "TROOT.h"
37 
38 //#define DEBUG
39 #ifdef DEBUG
40 #define NSAMPLE 10
41 #include <iostream>
42 #endif
43 
44 // need to implement integral option
45 
46 namespace ROOT {
47 
48  namespace Fit {
49 
50  namespace FitUtil {
51 
52  // derivative with respect of the parameter to be integrated
53  template<class GradFunc = IGradModelFunction>
54  struct ParamDerivFunc {
55  ParamDerivFunc(const GradFunc & f) : fFunc(f), fIpar(0) {}
56  void SetDerivComponent(unsigned int ipar) { fIpar = ipar; }
57  double operator() (const double *x, const double *p) const {
58  return fFunc.ParameterDerivative( x, p, fIpar );
59  }
60  unsigned int NDim() const { return fFunc.NDim(); }
61  const GradFunc & fFunc;
62  unsigned int fIpar;
63  };
64 
65 // simple gradient calculator using the 2 points rule
66 
67  class SimpleGradientCalculator {
68 
69  public:
70  // construct from function and gradient dimension gdim
71  // gdim = npar for parameter gradient
72  // gdim = ndim for coordinate gradients
73  // construct (the param values will be passed later)
74  // one can choose between 2 points rule (1 extra evaluation) istrat=1
75  // or two point rule (2 extra evaluation)
76  // (found 2 points rule does not work correctly - minuit2FitBench fails)
77  SimpleGradientCalculator(int gdim, const IModelFunction & func,double eps = 2.E-8, int istrat = 1) :
78  fEps(eps),
79  fPrecision(1.E-8 ), // sqrt(epsilon)
80  fStrategy(istrat),
81  fN(gdim ),
82  fFunc(func),
83  fVec(std::vector<double>(gdim) ) // this can be probably optimized
84  {}
85 
86  // internal method to calculate single partial derivative
87  // assume cached vector fVec is already set
88  double DoParameterDerivative(const double *x, const double *p, double f0, int k) const {
89  double p0 = p[k];
90  double h = std::max( fEps* std::abs(p0), 8.0*fPrecision*(std::abs(p0) + fPrecision) );
91  fVec[k] += h;
92  double deriv = 0;
93  // t.b.d : treat case of infinities
94  //if (fval > - std::numeric_limits<double>::max() && fval < std::numeric_limits<double>::max() )
95  double f1 = fFunc(x, &fVec.front() );
96  if (fStrategy > 1) {
97  fVec[k] = p0 - h;
98  double f2 = fFunc(x, &fVec.front() );
99  deriv = 0.5 * ( f2 - f1 )/h;
100  }
101  else
102  deriv = ( f1 - f0 )/h;
103 
104  fVec[k] = p[k]; // restore original p value
105  return deriv;
106  }
107  // number of dimension in x (needed when calculating the integrals)
108  unsigned int NDim() const {
109  return fFunc.NDim();
110  }
111  // number of parameters (needed for grad ccalculation)
112  unsigned int NPar() const {
113  return fFunc.NPar();
114  }
115 
116  double ParameterDerivative(const double *x, const double *p, int ipar) const {
117  // fVec are the cached parameter values
118  std::copy(p, p+fN, fVec.begin());
119  double f0 = fFunc(x, p);
120  return DoParameterDerivative(x,p,f0,ipar);
121  }
122 
123  // calculate all gradient at point (x,p) knnowing already value f0 (we gain a function eval.)
124  void ParameterGradient(const double * x, const double * p, double f0, double * g) {
125  // fVec are the cached parameter values
126  std::copy(p, p+fN, fVec.begin());
127  for (unsigned int k = 0; k < fN; ++k) {
128  g[k] = DoParameterDerivative(x,p,f0,k);
129  }
130  }
131 
132  // calculate gradient w.r coordinate values
133  void Gradient(const double * x, const double * p, double f0, double * g) {
134  // fVec are the cached coordinate values
135  std::copy(x, x+fN, fVec.begin());
136  for (unsigned int k = 0; k < fN; ++k) {
137  double x0 = x[k];
138  double h = std::max( fEps* std::abs(x0), 8.0*fPrecision*(std::abs(x0) + fPrecision) );
139  fVec[k] += h;
140  // t.b.d : treat case of infinities
141  //if (fval > - std::numeric_limits<double>::max() && fval < std::numeric_limits<double>::max() )
142  double f1 = fFunc( &fVec.front(), p );
143  if (fStrategy > 1) {
144  fVec[k] = x0 - h;
145  double f2 = fFunc( &fVec.front(), p );
146  g[k] = 0.5 * ( f2 - f1 )/h;
147  }
148  else
149  g[k] = ( f1 - f0 )/h;
150 
151  fVec[k] = x[k]; // restore original x value
152  }
153  }
154 
155  private:
156 
157  double fEps;
158  double fPrecision;
159  int fStrategy; // strategy in calculation ( =1 use 2 point rule( 1 extra func) , = 2 use r point rule)
160  unsigned int fN; // gradient dimension
161  const IModelFunction & fFunc;
162  mutable std::vector<double> fVec; // cached coordinates (or parameter values in case of gradientpar)
163  };
164 
165 
166  // function to avoid infinities or nan
167  double CorrectValue(double rval) {
168  // avoid infinities or nan in rval
169  if (rval > - std::numeric_limits<double>::max() && rval < std::numeric_limits<double>::max() )
170  return rval;
171  else if (rval < 0)
172  // case -inf
173  return -std::numeric_limits<double>::max();
174  else
175  // case + inf or nan
176  return + std::numeric_limits<double>::max();
177  }
178 
179  // Check if the value is a finite number. The argument rval is updated if it is infinite or NaN,
180  // setting it to the maximum finite value (preserving the sign).
181  bool CheckInfNaNValue(double &rval)
182  {
183  if (rval > - std::numeric_limits<double>::max() && rval < std::numeric_limits<double>::max() )
184  return true;
185  else if (rval < 0) {
186  // case -inf
187  rval = -std::numeric_limits<double>::max();
188  return false;
189  }
190  else {
191  // case + inf or nan
192  rval = + std::numeric_limits<double>::max();
193  return false;
194  }
195  }
196 
197 
198  // calculation of the integral of the gradient functions
199  // for a function providing derivative w.r.t parameters
200  // x1 and x2 defines the integration interval , p the parameters
201  template <class GFunc>
202  void CalculateGradientIntegral(const GFunc & gfunc,
203  const double *x1, const double * x2, const double * p, double *g) {
204 
205  // needs to calculate the integral for each partial derivative
206  ParamDerivFunc<GFunc> pfunc( gfunc);
207  IntegralEvaluator<ParamDerivFunc<GFunc> > igDerEval( pfunc, p, true);
208  // loop on the parameters
209  unsigned int npar = gfunc.NPar();
210  for (unsigned int k = 0; k < npar; ++k ) {
211  pfunc.SetDerivComponent(k);
212  g[k] = igDerEval( x1, x2 );
213  }
214  }
215 
216 
217 
218  } // end namespace FitUtil
219 
220 
221 
222 //___________________________________________________________________________________________________________________________
223 // for chi2 functions
224 //___________________________________________________________________________________________________________________________
225 
226  double FitUtil::EvaluateChi2(const IModelFunction &func, const BinData &data, const double *p, unsigned int &,
227  ROOT::Fit::ExecutionPolicy executionPolicy, unsigned nChunks)
228  {
229  // evaluate the chi2 given a function reference , the data and returns the value and also in nPoints
230  // the actual number of used points
231  // normal chi2 using only error on values (from fitting histogram)
232  // optionally the integral of function in the bin is used
233 
234  unsigned int n = data.Size();
235 
236  // set parameters of the function to cache integral value
237 #ifdef USE_PARAMCACHE
238  (const_cast<IModelFunction &>(func)).SetParameters(p);
239 #endif
240  // do not cache parameter values (it is not thread safe)
241  //func.SetParameters(p);
242 
243 
244  // get fit option and check case if using integral of bins
245  const DataOptions & fitOpt = data.Opt();
246  bool useBinIntegral = fitOpt.fIntegral && data.HasBinEdges();
247  bool useBinVolume = (fitOpt.fBinVolume && data.HasBinEdges());
248  bool useExpErrors = (fitOpt.fExpErrors);
249  bool isWeighted = data.IsWeighted();
250 
251 #ifdef DEBUG
252  std::cout << "\n\nFit data size = " << n << std::endl;
253  std::cout << "evaluate chi2 using function " << &func << " " << p << std::endl;
254  std::cout << "use empty bins " << fitOpt.fUseEmpty << std::endl;
255  std::cout << "use integral " << fitOpt.fIntegral << std::endl;
256  std::cout << "use all error=1 " << fitOpt.fErrors1 << std::endl;
257  if (isWeighted) std::cout << "Weighted data set - sumw = " << data.SumOfContent() << " sumw2 = " << data.SumOfError2() << std::endl;
258 #endif
259 
260  ROOT::Math::IntegrationOneDim::Type igType = ROOT::Math::IntegrationOneDim::kDEFAULT;
261  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
262  // do not use GSL integrator which is not thread safe
263  igType = ROOT::Math::IntegrationOneDim::kGAUSS;
264  }
265 #ifdef USE_PARAMCACHE
266  IntegralEvaluator<> igEval( func, 0, useBinIntegral, igType);
267 #else
268  IntegralEvaluator<> igEval( func, p, useBinIntegral, igType);
269 #endif
270  double maxResValue = std::numeric_limits<double>::max() /n;
271  double wrefVolume = 1.0;
272  if (useBinVolume) {
273  if (fitOpt.fNormBinVolume) wrefVolume /= data.RefVolume();
274  }
275 
276  (const_cast<IModelFunction &>(func)).SetParameters(p);
277 
278  auto mapFunction = [&](const unsigned i){
279 
280  double chi2{};
281  double fval{};
282 
283  const auto x1 = data.GetCoordComponent(i, 0);
284  const auto y = data.Value(i);
285  auto invError = data.InvError(i);
286 
287  //invError = (invError!= 0.0) ? 1.0/invError :1;
288 
289  const double * x = nullptr;
290  std::vector<double> xc;
291  double binVolume = 1.0;
292  if (useBinVolume) {
293  unsigned int ndim = data.NDim();
294  const double * x2 = data.BinUpEdge(i);
295  xc.resize(data.NDim());
296  for (unsigned int j = 0; j < ndim; ++j) {
297  auto xx = *data.GetCoordComponent(i, j);
298  binVolume *= std::abs(x2[j]- xx);
299  xc[j] = 0.5*(x2[j]+ xx);
300  }
301  x = xc.data();
302  // normalize the bin volume using a reference value
303  binVolume *= wrefVolume;
304  } else if(data.NDim() > 1) {
305  xc.resize(data.NDim());
306  xc[0] = *x1;
307  for (unsigned int j = 1; j < data.NDim(); ++j)
308  xc[j] = *data.GetCoordComponent(i, j);
309  x = xc.data();
310  } else {
311  x = x1;
312  }
313 
314 
315  if (!useBinIntegral) {
316 #ifdef USE_PARAMCACHE
317  fval = func ( x );
318 #else
319  fval = func ( x, p );
320 #endif
321  }
322  else {
323  // calculate integral normalized by bin volume
324  // need to set function and parameters here in case loop is parallelized
325  fval = igEval( x, data.BinUpEdge(i)) ;
326  }
327  // normalize result if requested according to bin volume
328  if (useBinVolume) fval *= binVolume;
329 
330  // expected errors
331  if (useExpErrors) {
332  double invWeight = 1.0;
333  if (isWeighted) {
334  // we need first to check if a weight factor needs to be applied
335  // weight = sumw2/sumw = error**2/content
336  //invWeight = y * invError * invError;
337  // we use always the global weight and not the observed one in the bin
338  // for empty bins use global weight (if it is weighted data.SumError2() is not zero)
339  invWeight = data.SumOfContent()/ data.SumOfError2();
340  //if (invError > 0) invWeight = y * invError * invError;
341  }
342 
343  // if (invError == 0) invWeight = (data.SumOfError2() > 0) ? data.SumOfContent()/ data.SumOfError2() : 1.0;
344  // compute expected error as f(x) / weight
345  double invError2 = (fval > 0) ? invWeight / fval : 0.0;
346  invError = std::sqrt(invError2);
347  //std::cout << "using Pearson chi2 " << x[0] << " " << 1./invError2 << " " << fval << std::endl;
348  }
349 
350 //#define DEBUG
351 #ifdef DEBUG
352  std::cout << x[0] << " " << y << " " << 1./invError << " params : ";
353  for (unsigned int ipar = 0; ipar < func.NPar(); ++ipar)
354  std::cout << p[ipar] << "\t";
355  std::cout << "\tfval = " << fval << " bin volume " << binVolume << " ref " << wrefVolume << std::endl;
356 #endif
357 //#undef DEBUG
358 
359  if (invError > 0) {
360 
361  double tmp = ( y -fval )* invError;
362  double resval = tmp * tmp;
363 
364 
365  // avoid inifinity or nan in chi2 values due to wrong function values
366  if ( resval < maxResValue )
367  chi2 += resval;
368  else {
369  //nRejected++;
370  chi2 += maxResValue;
371  }
372  }
373  return chi2;
374  };
375 
376 #ifdef R__USE_IMT
377  auto redFunction = [](const std::vector<double> & objs){
378  return std::accumulate(objs.begin(), objs.end(), double{});
379  };
380 #else
381  (void)nChunks;
382 
383  // If IMT is disabled, force the execution policy to the serial case
384  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
385  Warning("FitUtil::EvaluateChi2", "Multithread execution policy requires IMT, which is disabled. Changing "
386  "to ROOT::Fit::ExecutionPolicy::kSerial.");
387  executionPolicy = ROOT::Fit::ExecutionPolicy::kSerial;
388  }
389 #endif
390 
391  double res{};
392  if(executionPolicy == ROOT::Fit::ExecutionPolicy::kSerial){
393  for (unsigned int i=0; i<n; ++i) {
394  res += mapFunction(i);
395  }
396 #ifdef R__USE_IMT
397  } else if(executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
398  ROOT::TThreadExecutor pool;
399  auto chunks = nChunks !=0? nChunks: setAutomaticChunking(data.Size());
400  res = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, n), redFunction, chunks);
401 #endif
402 // } else if(executionPolicy == ROOT::Fit::kMultitProcess){
403  // ROOT::TProcessExecutor pool;
404  // res = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, n), redFunction);
405  } else{
406  Error("FitUtil::EvaluateChi2","Execution policy unknown. Avalaible choices:\n ROOT::Fit::ExecutionPolicy::kSerial (default)\n ROOT::Fit::ExecutionPolicy::kMultithread (requires IMT)\n");
407  }
408 
409  return res;
410 }
411 
412 
413 //___________________________________________________________________________________________________________________________
414 
415 double FitUtil::EvaluateChi2Effective(const IModelFunction & func, const BinData & data, const double * p, unsigned int & nPoints) {
416  // evaluate the chi2 given a function reference , the data and returns the value and also in nPoints
417  // the actual number of used points
418  // method using the error in the coordinates
419  // integral of bin does not make sense in this case
420 
421  unsigned int n = data.Size();
422 
423 #ifdef DEBUG
424  std::cout << "\n\nFit data size = " << n << std::endl;
425  std::cout << "evaluate effective chi2 using function " << &func << " " << p << std::endl;
426 #endif
427 
428  assert(data.HaveCoordErrors() || data.HaveAsymErrors());
429 
430  double chi2 = 0;
431  //int nRejected = 0;
432 
433 
434  //func.SetParameters(p);
435 
436  unsigned int ndim = func.NDim();
437 
438  // use Richardson derivator
439  ROOT::Math::RichardsonDerivator derivator;
440 
441  double maxResValue = std::numeric_limits<double>::max() /n;
442 
443 
444 
445  for (unsigned int i = 0; i < n; ++ i) {
446 
447 
448  double y = 0;
449  const double * x = data.GetPoint(i,y);
450 
451  double fval = func( x, p );
452 
453  double delta_y_func = y - fval;
454 
455 
456  double ey = 0;
457  const double * ex = 0;
458  if (!data.HaveAsymErrors() )
459  ex = data.GetPointError(i, ey);
460  else {
461  double eylow, eyhigh = 0;
462  ex = data.GetPointError(i, eylow, eyhigh);
463  if ( delta_y_func < 0)
464  ey = eyhigh; // function is higher than points
465  else
466  ey = eylow;
467  }
468  double e2 = ey * ey;
469  // before calculating the gradient check that all error in x are not zero
470  unsigned int j = 0;
471  while ( j < ndim && ex[j] == 0.) { j++; }
472  // if j is less ndim some elements are not zero
473  if (j < ndim) {
474  // need an adapter from a multi-dim function to a one-dimensional
475  ROOT::Math::OneDimMultiFunctionAdapter<const IModelFunction &> f1D(func,x,0,p);
476  // select optimal step size (use 10--2 by default as was done in TF1:
477  double kEps = 0.01;
478  double kPrecision = 1.E-8;
479  for (unsigned int icoord = 0; icoord < ndim; ++icoord) {
480  // calculate derivative for each coordinate
481  if (ex[icoord] > 0) {
482  //gradCalc.Gradient(x, p, fval, &grad[0]);
483  f1D.SetCoord(icoord);
484  // optimal spep size (take ex[] as scale for the points and 1% of it
485  double x0= x[icoord];
486  double h = std::max( kEps* std::abs(ex[icoord]), 8.0*kPrecision*(std::abs(x0) + kPrecision) );
487  double deriv = derivator.Derivative1(f1D, x[icoord], h);
488  double edx = ex[icoord] * deriv;
489  e2 += edx * edx;
490 #ifdef DEBUG
491  std::cout << "error for coord " << icoord << " = " << ex[icoord] << " deriv " << deriv << std::endl;
492 #endif
493  }
494  }
495  }
496  double w2 = (e2 > 0) ? 1.0/e2 : 0;
497  double resval = w2 * ( y - fval ) * ( y - fval);
498 
499 #ifdef DEBUG
500  std::cout << x[0] << " " << y << " ex " << ex[0] << " ey " << ey << " params : ";
501  for (unsigned int ipar = 0; ipar < func.NPar(); ++ipar)
502  std::cout << p[ipar] << "\t";
503  std::cout << "\tfval = " << fval << "\tresval = " << resval << std::endl;
504 #endif
505 
506  // avoid (infinity and nan ) in the chi2 sum
507  // eventually add possibility of excluding some points (like singularity)
508  if ( resval < maxResValue )
509  chi2 += resval;
510  else
511  chi2 += maxResValue;
512  //nRejected++;
513 
514  }
515 
516  // reset the number of fitting data points
517  nPoints = n; // no points are rejected
518  //if (nRejected != 0) nPoints = n - nRejected;
519 
520 #ifdef DEBUG
521  std::cout << "chi2 = " << chi2 << " n = " << nPoints << std::endl;
522 #endif
523 
524  return chi2;
525 
526 }
527 
528 
529 ////////////////////////////////////////////////////////////////////////////////
530 /// evaluate the chi2 contribution (residual term) only for data with no coord-errors
531 /// This function is used in the specialized least square algorithms like FUMILI or L.M.
532 /// if we have error on the coordinates the method is not yet implemented
533 /// integral option is also not yet implemented
534 /// one can use in that case normal chi2 method
535 
536 double FitUtil::EvaluateChi2Residual(const IModelFunction & func, const BinData & data, const double * p, unsigned int i, double * g) {
537  if (data.GetErrorType() == BinData::kCoordError && data.Opt().fCoordErrors ) {
538  MATH_ERROR_MSG("FitUtil::EvaluateChi2Residual","Error on the coordinates are not used in calculating Chi2 residual");
539  return 0; // it will assert otherwise later in GetPoint
540  }
541 
542 
543  //func.SetParameters(p);
544 
545  double y, invError = 0;
546 
547  const DataOptions & fitOpt = data.Opt();
548  bool useBinIntegral = fitOpt.fIntegral && data.HasBinEdges();
549  bool useBinVolume = (fitOpt.fBinVolume && data.HasBinEdges());
550  bool useExpErrors = (fitOpt.fExpErrors);
551 
552  const double * x1 = data.GetPoint(i,y, invError);
553 
554  IntegralEvaluator<> igEval( func, p, useBinIntegral);
555  double fval = 0;
556  unsigned int ndim = data.NDim();
557  double binVolume = 1.0;
558  const double * x2 = 0;
559  if (useBinVolume || useBinIntegral) x2 = data.BinUpEdge(i);
560 
561  double * xc = 0;
562 
563  if (useBinVolume) {
564  xc = new double[ndim];
565  for (unsigned int j = 0; j < ndim; ++j) {
566  binVolume *= std::abs( x2[j]-x1[j] );
567  xc[j] = 0.5*(x2[j]+ x1[j]);
568  }
569  // normalize the bin volume using a reference value
570  binVolume /= data.RefVolume();
571  }
572 
573  const double * x = (useBinVolume) ? xc : x1;
574 
575  if (!useBinIntegral) {
576  fval = func ( x, p );
577  }
578  else {
579  // calculate integral (normalized by bin volume)
580  // need to set function and parameters here in case loop is parallelized
581  fval = igEval( x1, x2) ;
582  }
583  // normalize result if requested according to bin volume
584  if (useBinVolume) fval *= binVolume;
585 
586  // expected errors
587  if (useExpErrors) {
588  // we need first to check if a weight factor needs to be applied
589  // weight = sumw2/sumw = error**2/content
590  //NOTE: assume histogram is not weighted
591  // don't know how to do with bins with weight = 0
592  //double invWeight = y * invError * invError;
593  // if (invError == 0) invWeight = (data.SumOfError2() > 0) ? data.SumOfContent()/ data.SumOfError2() : 1.0;
594  // compute expected error as f(x) / weight
595  double invError2 = (fval > 0) ? 1.0 / fval : 0.0;
596  invError = std::sqrt(invError2);
597  }
598 
599 
600  double resval = ( y -fval )* invError;
601 
602  // avoid infinities or nan in resval
603  resval = CorrectValue(resval);
604 
605  // estimate gradient
606  if (g != 0) {
607 
608  unsigned int npar = func.NPar();
609  const IGradModelFunction * gfunc = dynamic_cast<const IGradModelFunction *>( &func);
610 
611  if (gfunc != 0) {
612  //case function provides gradient
613  if (!useBinIntegral ) {
614  gfunc->ParameterGradient( x , p, g );
615  }
616  else {
617  // needs to calculate the integral for each partial derivative
618  CalculateGradientIntegral( *gfunc, x1, x2, p, g);
619  }
620  }
621  else {
622  SimpleGradientCalculator gc( npar, func);
623  if (!useBinIntegral )
624  gc.ParameterGradient(x, p, fval, g);
625  else {
626  // needs to calculate the integral for each partial derivative
627  CalculateGradientIntegral( gc, x1, x2, p, g);
628  }
629  }
630  // mutiply by - 1 * weight
631  for (unsigned int k = 0; k < npar; ++k) {
632  g[k] *= - invError;
633  if (useBinVolume) g[k] *= binVolume;
634  }
635  }
636 
637  if (useBinVolume) delete [] xc;
638 
639  return resval;
640 
641 }
642 
643 void FitUtil::EvaluateChi2Gradient(const IModelFunction &f, const BinData &data, const double *p, double *grad,
644  unsigned int &nPoints, ROOT::Fit::ExecutionPolicy executionPolicy, unsigned nChunks)
645 {
646  // evaluate the gradient of the chi2 function
647  // this function is used when the model function knows how to calculate the derivative and we can
648  // avoid that the minimizer re-computes them
649  //
650  // case of chi2 effective (errors on coordinate) is not supported
651 
652  if (data.HaveCoordErrors()) {
653  MATH_ERROR_MSG("FitUtil::EvaluateChi2Gradient",
654  "Error on the coordinates are not used in calculating Chi2 gradient");
655  return; // it will assert otherwise later in GetPoint
656  }
657 
658  const IGradModelFunction *fg = dynamic_cast<const IGradModelFunction *>(&f);
659  assert(fg != nullptr); // must be called by a gradient function
660 
661  const IGradModelFunction &func = *fg;
662 
663 #ifdef DEBUG
664  std::cout << "\n\nFit data size = " << n << std::endl;
665  std::cout << "evaluate chi2 using function gradient " << &func << " " << p << std::endl;
666 #endif
667 
668  const DataOptions &fitOpt = data.Opt();
669  bool useBinIntegral = fitOpt.fIntegral && data.HasBinEdges();
670  bool useBinVolume = (fitOpt.fBinVolume && data.HasBinEdges());
671 
672  double wrefVolume = 1.0;
673  if (useBinVolume) {
674  if (fitOpt.fNormBinVolume) wrefVolume /= data.RefVolume();
675  }
676 
677  ROOT::Math::IntegrationOneDim::Type igType = ROOT::Math::IntegrationOneDim::kDEFAULT;
678  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
679  // do not use GSL integrator which is not thread safe
680  igType = ROOT::Math::IntegrationOneDim::kGAUSS;
681  }
682  IntegralEvaluator<> igEval(func, p, useBinIntegral,igType);
683 
684  unsigned int npar = func.NPar();
685  unsigned initialNPoints = data.Size();
686 
687  std::vector<bool> isPointRejected(initialNPoints);
688 
689  auto mapFunction = [&](const unsigned int i) {
690  // set all vector values to zero
691  std::vector<double> gradFunc(npar);
692  std::vector<double> pointContribution(npar);
693 
694  const auto x1 = data.GetCoordComponent(i, 0);
695  const auto y = data.Value(i);
696  auto invError = data.Error(i);
697 
698  invError = (invError != 0.0) ? 1.0 / invError : 1;
699 
700  double fval = 0;
701 
702  const double *x = nullptr;
703  std::vector<double> xc;
704 
705  unsigned int ndim = data.NDim();
706  double binVolume = 1;
707  if (useBinVolume) {
708  const double *x2 = data.BinUpEdge(i);
709 
710  xc.resize(ndim);
711  for (unsigned int j = 0; j < ndim; ++j) {
712  auto x1_j = *data.GetCoordComponent(i, j);
713  binVolume *= std::abs(x2[j] - x1_j);
714  xc[j] = 0.5 * (x2[j] + x1_j);
715  }
716 
717  x = xc.data();
718 
719  // normalize the bin volume using a reference value
720  binVolume *= wrefVolume;
721  } else if (ndim > 1) {
722  xc.resize(ndim);
723  xc[0] = *x1;
724  for (unsigned int j = 1; j < ndim; ++j)
725  xc[j] = *data.GetCoordComponent(i, j);
726  x = xc.data();
727  } else {
728  x = x1;
729  }
730 
731  if (!useBinIntegral) {
732  fval = func(x, p);
733  func.ParameterGradient(x, p, &gradFunc[0]);
734  } else {
735  auto x2 = data.BinUpEdge(i);
736  // calculate normalized integral and gradient (divided by bin volume)
737  // need to set function and parameters here in case loop is parallelized
738  fval = igEval(x, x2);
739  CalculateGradientIntegral(func, x, x2, p, &gradFunc[0]);
740  }
741  if (useBinVolume)
742  fval *= binVolume;
743 
744 #ifdef DEBUG
745  std::cout << x[0] << " " << y << " " << 1. / invError << " params : ";
746  for (unsigned int ipar = 0; ipar < npar; ++ipar)
747  std::cout << p[ipar] << "\t";
748  std::cout << "\tfval = " << fval << std::endl;
749 #endif
750  if (!CheckInfNaNValue(fval)) {
751  isPointRejected[i] = true;
752  // Return a zero contribution to all partial derivatives on behalf of the current point
753  return pointContribution;
754  }
755 
756  // loop on the parameters
757  unsigned int ipar = 0;
758  for (; ipar < npar; ++ipar) {
759 
760  // correct gradient for bin volumes
761  if (useBinVolume)
762  gradFunc[ipar] *= binVolume;
763 
764  // avoid singularity in the function (infinity and nan ) in the chi2 sum
765  // eventually add possibility of excluding some points (like singularity)
766  double dfval = gradFunc[ipar];
767  if (!CheckInfNaNValue(dfval)) {
768  break; // exit loop on parameters
769  }
770 
771  // calculate derivative point contribution
772  pointContribution[ipar] = -2.0 * (y - fval) * invError * invError * gradFunc[ipar];
773  }
774 
775  if (ipar < npar) {
776  // case loop was broken for an overflow in the gradient calculation
777  isPointRejected[i] = true;
778  }
779 
780  return pointContribution;
781  };
782 
783  // Vertically reduce the set of vectors by summing its equally-indexed components
784  auto redFunction = [&](const std::vector<std::vector<double>> &pointContributions) {
785  std::vector<double> result(npar);
786 
787  for (auto const &pointContribution : pointContributions) {
788  for (unsigned int parameterIndex = 0; parameterIndex < npar; parameterIndex++)
789  result[parameterIndex] += pointContribution[parameterIndex];
790  }
791 
792  return result;
793  };
794 
795  std::vector<double> g(npar);
796 
797 #ifndef R__USE_IMT
798  // If IMT is disabled, force the execution policy to the serial case
799  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
800  Warning("FitUtil::EvaluateChi2Gradient", "Multithread execution policy requires IMT, which is disabled. Changing "
801  "to ROOT::Fit::ExecutionPolicy::kSerial.");
802  executionPolicy = ROOT::Fit::ExecutionPolicy::kSerial;
803  }
804 #endif
805 
806  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kSerial) {
807  std::vector<std::vector<double>> allGradients(initialNPoints);
808  for (unsigned int i = 0; i < initialNPoints; ++i) {
809  allGradients[i] = mapFunction(i);
810  }
811  g = redFunction(allGradients);
812  }
813 #ifdef R__USE_IMT
814  else if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
815  ROOT::TThreadExecutor pool;
816  auto chunks = nChunks != 0 ? nChunks : setAutomaticChunking(initialNPoints);
817  g = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, initialNPoints), redFunction, chunks);
818  }
819 #endif
820  // else if(executionPolicy == ROOT::Fit::kMultiprocess){
821  // ROOT::TProcessExecutor pool;
822  // g = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, n), redFunction);
823  // }
824  else {
825  Error("FitUtil::EvaluateChi2Gradient",
826  "Execution policy unknown. Avalaible choices:\n 0: Serial (default)\n 1: MultiThread (requires IMT)\n");
827  }
828 
829 #ifndef R__USE_IMT
830  //to fix compiler warning
831  (void)nChunks;
832 #endif
833 
834  // correct the number of points
835  nPoints = initialNPoints;
836 
837  if (std::any_of(isPointRejected.begin(), isPointRejected.end(), [](bool point) { return point; })) {
838  unsigned nRejected = std::accumulate(isPointRejected.begin(), isPointRejected.end(), 0);
839  assert(nRejected <= initialNPoints);
840  nPoints = initialNPoints - nRejected;
841 
842  if (nPoints < npar)
843  MATH_ERROR_MSG("FitUtil::EvaluateChi2Gradient",
844  "Error - too many points rejected for overflow in gradient calculation");
845  }
846 
847  // copy result
848  std::copy(g.begin(), g.end(), grad);
849 }
850 
851 //______________________________________________________________________________________________________
852 //
853 // Log Likelihood functions
854 //_______________________________________________________________________________________________________
855 
856 // utility function used by the likelihoods
857 
858 // for LogLikelihood functions
859 
860 double FitUtil::EvaluatePdf(const IModelFunction & func, const UnBinData & data, const double * p, unsigned int i, double * g) {
861  // evaluate the pdf contribution to the generic logl function in case of bin data
862  // return actually the log of the pdf and its derivatives
863 
864 
865  //func.SetParameters(p);
866 
867 
868  const double * x = data.Coords(i);
869  double fval = func ( x, p );
870  double logPdf = ROOT::Math::Util::EvalLog(fval);
871  //return
872  if (g == 0) return logPdf;
873 
874  const IGradModelFunction * gfunc = dynamic_cast<const IGradModelFunction *>( &func);
875 
876  // gradient calculation
877  if (gfunc != 0) {
878  //case function provides gradient
879  gfunc->ParameterGradient( x , p, g );
880  }
881  else {
882  // estimate gradieant numerically with simple 2 point rule
883  // should probably calculate gradient of log(pdf) is more stable numerically
884  SimpleGradientCalculator gc(func.NPar(), func);
885  gc.ParameterGradient(x, p, fval, g );
886  }
887  // divide gradient by function value since returning the logs
888  for (unsigned int ipar = 0; ipar < func.NPar(); ++ipar) {
889  g[ipar] /= fval; // this should be checked against infinities
890  }
891 
892 #ifdef DEBUG
893  std::cout << x[i] << "\t";
894  std::cout << "\tpar = [ " << func.NPar() << " ] = ";
895  for (unsigned int ipar = 0; ipar < func.NPar(); ++ipar)
896  std::cout << p[ipar] << "\t";
897  std::cout << "\tfval = " << fval;
898  std::cout << "\tgrad = [ ";
899  for (unsigned int ipar = 0; ipar < func.NPar(); ++ipar)
900  std::cout << g[ipar] << "\t";
901  std::cout << " ] " << std::endl;
902 #endif
903 
904 
905  return logPdf;
906 }
907 
908 double FitUtil::EvaluateLogL(const IModelFunctionTempl<double> &func, const UnBinData &data, const double *p,
909  int iWeight, bool extended, unsigned int &nPoints,
910  ROOT::Fit::ExecutionPolicy executionPolicy, unsigned nChunks)
911 {
912  // evaluate the LogLikelihood
913 
914  unsigned int n = data.Size();
915 
916  //unsigned int nRejected = 0;
917 
918  bool normalizeFunc = false;
919 
920  // set parameters of the function to cache integral value
921 #ifdef USE_PARAMCACHE
922  (const_cast<IModelFunctionTempl<double> &>(func)).SetParameters(p);
923 #endif
924 
925  nPoints = data.Size(); // npoints
926 
927 #ifdef R__USE_IMT
928  // in case parameter needs to be propagated to user function use trick to set parameters by calling one time the function
929  // this will be done in sequential mode and parameters can be set in a thread safe manner
930  if (!normalizeFunc) {
931  if (data.NDim() == 1) {
932  const double * x = data.GetCoordComponent(0,0);
933  func( x, p);
934  }
935  else {
936  std::vector<double> x(data.NDim());
937  for (unsigned int j = 0; j < data.NDim(); ++j)
938  x[j] = *data.GetCoordComponent(0, j);
939  func( x.data(), p);
940  }
941  }
942 #endif
943 
944  double norm = 1.0;
945  if (normalizeFunc) {
946  // compute integral of the function
947  std::vector<double> xmin(data.NDim());
948  std::vector<double> xmax(data.NDim());
949  IntegralEvaluator<> igEval(func, p, true);
950  // compute integral in the ranges where is defined
951  if (data.Range().Size() > 0) {
952  norm = 0;
953  for (unsigned int ir = 0; ir < data.Range().Size(); ++ir) {
954  data.Range().GetRange(&xmin[0], &xmax[0], ir);
955  norm += igEval.Integral(xmin.data(), xmax.data());
956  }
957  } else {
958  // use (-inf +inf)
959  data.Range().GetRange(&xmin[0], &xmax[0]);
960  // check if funcition is zero at +- inf
961  if (func(xmin.data(), p) != 0 || func(xmax.data(), p) != 0) {
962  MATH_ERROR_MSG("FitUtil::EvaluateLogLikelihood",
963  "A range has not been set and the function is not zero at +/- inf");
964  return 0;
965  }
966  norm = igEval.Integral(&xmin[0], &xmax[0]);
967  }
968  }
969 
970  // needed to compue effective global weight in case of extended likelihood
971 
972  auto mapFunction = [&](const unsigned i) {
973  double W = 0;
974  double W2 = 0;
975  double fval = 0;
976 
977  if (data.NDim() > 1) {
978  std::vector<double> x(data.NDim());
979  for (unsigned int j = 0; j < data.NDim(); ++j)
980  x[j] = *data.GetCoordComponent(i, j);
981 #ifdef USE_PARAMCACHE
982  fval = func(x.data());
983 #else
984  fval = func(x.data(), p);
985 #endif
986 
987  // one -dim case
988  } else {
989  const auto x = data.GetCoordComponent(i, 0);
990 #ifdef USE_PARAMCACHE
991  fval = func(x);
992 #else
993  fval = func(x, p);
994 #endif
995  }
996 
997  if (normalizeFunc)
998  fval = fval * (1 / norm);
999 
1000  // function EvalLog protects against negative or too small values of fval
1001  double logval = ROOT::Math::Util::EvalLog(fval);
1002  if (iWeight > 0) {
1003  double weight = data.Weight(i);
1004  logval *= weight;
1005  if (iWeight == 2) {
1006  logval *= weight; // use square of weights in likelihood
1007  if (!extended) {
1008  // needed sum of weights and sum of weight square if likelkihood is extended
1009  W = weight;
1010  W2 = weight * weight;
1011  }
1012  }
1013  }
1014  return LikelihoodAux<double>(logval, W, W2);
1015  };
1016 
1017 #ifdef R__USE_IMT
1018  // auto redFunction = [](const std::vector<LikelihoodAux<double>> & objs){
1019  // return std::accumulate(objs.begin(), objs.end(), LikelihoodAux<double>(0.0,0.0,0.0),
1020  // [](const LikelihoodAux<double> &l1, const LikelihoodAux<double> &l2){
1021  // return l1+l2;
1022  // });
1023  // };
1024  // do not use std::accumulate to be sure to maintain always the same order
1025  auto redFunction = [](const std::vector<LikelihoodAux<double>> & objs){
1026  auto l0 = LikelihoodAux<double>(0.0,0.0,0.0);
1027  for ( auto & l : objs ) {
1028  l0 = l0 + l;
1029  }
1030  return l0;
1031  };
1032 #else
1033  (void)nChunks;
1034 
1035  // If IMT is disabled, force the execution policy to the serial case
1036  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
1037  Warning("FitUtil::EvaluateLogL", "Multithread execution policy requires IMT, which is disabled. Changing "
1038  "to ROOT::Fit::ExecutionPolicy::kSerial.");
1039  executionPolicy = ROOT::Fit::ExecutionPolicy::kSerial;
1040  }
1041 #endif
1042 
1043  double logl{};
1044  double sumW{};
1045  double sumW2{};
1046  if(executionPolicy == ROOT::Fit::ExecutionPolicy::kSerial){
1047  for (unsigned int i=0; i<n; ++i) {
1048  auto resArray = mapFunction(i);
1049  logl+=resArray.logvalue;
1050  sumW+=resArray.weight;
1051  sumW2+=resArray.weight2;
1052  }
1053 #ifdef R__USE_IMT
1054  } else if(executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
1055  ROOT::TThreadExecutor pool;
1056  auto chunks = nChunks !=0? nChunks: setAutomaticChunking(data.Size());
1057  auto resArray = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, n), redFunction, chunks);
1058  logl=resArray.logvalue;
1059  sumW=resArray.weight;
1060  sumW2=resArray.weight2;
1061 #endif
1062 // } else if(executionPolicy == ROOT::Fit::kMultitProcess){
1063  // ROOT::TProcessExecutor pool;
1064  // res = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, n), redFunction);
1065  } else{
1066  Error("FitUtil::EvaluateLogL","Execution policy unknown. Avalaible choices:\n ROOT::Fit::ExecutionPolicy::kSerial (default)\n ROOT::Fit::ExecutionPolicy::kMultithread (requires IMT)\n");
1067  }
1068 
1069  if (extended) {
1070  // add Poisson extended term
1071  double extendedTerm = 0; // extended term in likelihood
1072  double nuTot = 0;
1073  // nuTot is integral of function in the range
1074  // if function has been normalized integral has been already computed
1075  if (!normalizeFunc) {
1076  IntegralEvaluator<> igEval( func, p, true);
1077  std::vector<double> xmin(data.NDim());
1078  std::vector<double> xmax(data.NDim());
1079 
1080  // compute integral in the ranges where is defined
1081  if (data.Range().Size() > 0 ) {
1082  nuTot = 0;
1083  for (unsigned int ir = 0; ir < data.Range().Size(); ++ir) {
1084  data.Range().GetRange(&xmin[0],&xmax[0],ir);
1085  nuTot += igEval.Integral(xmin.data(),xmax.data());
1086  }
1087  } else {
1088  // use (-inf +inf)
1089  data.Range().GetRange(&xmin[0],&xmax[0]);
1090  // check if funcition is zero at +- inf
1091  if (func(xmin.data(), p) != 0 || func(xmax.data(), p) != 0) {
1092  MATH_ERROR_MSG("FitUtil::EvaluateLogLikelihood","A range has not been set and the function is not zero at +/- inf");
1093  return 0;
1094  }
1095  nuTot = igEval.Integral(&xmin[0],&xmax[0]);
1096  }
1097 
1098  // force to be last parameter value
1099  //nutot = p[func.NDim()-1];
1100  if (iWeight != 2)
1101  extendedTerm = - nuTot; // no need to add in this case n log(nu) since is already computed before
1102  else {
1103  // case use weight square in likelihood : compute total effective weight = sw2/sw
1104  // ignore for the moment case when sumW is zero
1105  extendedTerm = - (sumW2 / sumW) * nuTot;
1106  }
1107 
1108  }
1109  else {
1110  nuTot = norm;
1111  extendedTerm = - nuTot + double(n) * ROOT::Math::Util::EvalLog( nuTot);
1112  // in case of weights need to use here sum of weights (to be done)
1113  }
1114  logl += extendedTerm;
1115 
1116  }
1117 
1118 #ifdef DEBUG
1119  std::cout << "Evaluated log L for parameters (";
1120  for (unsigned int ip = 0; ip < func.NPar(); ++ip)
1121  std::cout << " " << p[ip];
1122  std::cout << ") fval = " << -logl << std::endl;
1123 #endif
1124 
1125  return -logl;
1126 }
1127 
1128 void FitUtil::EvaluateLogLGradient(const IModelFunction &f, const UnBinData &data, const double *p, double *grad,
1129  unsigned int &, ROOT::Fit::ExecutionPolicy executionPolicy, unsigned nChunks)
1130 {
1131  // evaluate the gradient of the log likelihood function
1132 
1133  const IGradModelFunction *fg = dynamic_cast<const IGradModelFunction *>(&f);
1134  assert(fg != nullptr); // must be called by a grad function
1135 
1136  const IGradModelFunction &func = *fg;
1137 
1138  unsigned int npar = func.NPar();
1139  unsigned initialNPoints = data.Size();
1140 
1141  (const_cast<IGradModelFunction &>(func)).SetParameters(p);
1142 
1143 #ifdef DEBUG
1144  std::cout << "\n===> Evaluate Gradient for parameters ";
1145  for (unsigned int ip = 0; ip < npar; ++ip)
1146  std::cout << " " << p[ip];
1147  std::cout << "\n";
1148 #endif
1149 
1150  const double kdmax1 = std::sqrt(std::numeric_limits<double>::max());
1151  const double kdmax2 = std::numeric_limits<double>::max() / (4 * initialNPoints);
1152 
1153  auto mapFunction = [&](const unsigned int i) {
1154  std::vector<double> gradFunc(npar);
1155  std::vector<double> pointContribution(npar);
1156 
1157 
1158  const double * x = nullptr;
1159  std::vector<double> xc;
1160  if (data.NDim() > 1) {
1161  xc.resize(data.NDim() );
1162  for (unsigned int j = 0; j < data.NDim(); ++j)
1163  xc[j] = *data.GetCoordComponent(i, j);
1164  x = xc.data();
1165  } else {
1166  x = data.GetCoordComponent(i, 0);
1167  }
1168 
1169  double fval = func(x, p);
1170  func.ParameterGradient(x, p, &gradFunc[0]);
1171 
1172 #ifdef DEBUG
1173  {
1174  R__LOCKGUARD(gROOTMutex);
1175  if (i < 5 || (i > data.Size()-5) ) {
1176  if (data.NDim() > 1) std::cout << i << " x " << x[0] << " y " << x[1] << " func " << fval
1177  << " gradient " << gradFunc[0] << " " << gradFunc[1] << " " << gradFunc[3] << std::endl;
1178  else std::cout << i << " x " << x[0] << " gradient " << gradFunc[0] << " " << gradFunc[1] << " " << gradFunc[3] << std::endl;
1179  }
1180  }
1181 #endif
1182 
1183  for (unsigned int kpar = 0; kpar < npar; ++kpar) {
1184  if (fval > 0)
1185  pointContribution[kpar] = -1. / fval * gradFunc[kpar];
1186  else if (gradFunc[kpar] != 0) {
1187  double gg = kdmax1 * gradFunc[kpar];
1188  if (gg > 0)
1189  gg = std::min(gg, kdmax2);
1190  else
1191  gg = std::max(gg, -kdmax2);
1192  pointContribution[kpar] = -gg;
1193  }
1194  // if func derivative is zero term is also zero so do not add in g[kpar]
1195  }
1196 
1197  return pointContribution;
1198  };
1199 
1200  // Vertically reduce the set of vectors by summing its equally-indexed components
1201  auto redFunction = [&](const std::vector<std::vector<double>> &pointContributions) {
1202  std::vector<double> result(npar);
1203 
1204  for (auto const &pointContribution : pointContributions) {
1205  for (unsigned int parameterIndex = 0; parameterIndex < npar; parameterIndex++)
1206  result[parameterIndex] += pointContribution[parameterIndex];
1207  }
1208 
1209  return result;
1210  };
1211 
1212  std::vector<double> g(npar);
1213 
1214 #ifndef R__USE_IMT
1215  // If IMT is disabled, force the execution policy to the serial case
1216  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
1217  Warning("FitUtil::EvaluateLogLGradient", "Multithread execution policy requires IMT, which is disabled. Changing "
1218  "to ROOT::Fit::ExecutionPolicy::kSerial.");
1219  executionPolicy = ROOT::Fit::ExecutionPolicy::kSerial;
1220  }
1221 #endif
1222 
1223  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kSerial) {
1224  std::vector<std::vector<double>> allGradients(initialNPoints);
1225  for (unsigned int i = 0; i < initialNPoints; ++i) {
1226  allGradients[i] = mapFunction(i);
1227  }
1228  g = redFunction(allGradients);
1229  }
1230 #ifdef R__USE_IMT
1231  else if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
1232  ROOT::TThreadExecutor pool;
1233  auto chunks = nChunks != 0 ? nChunks : setAutomaticChunking(initialNPoints);
1234  g = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, initialNPoints), redFunction, chunks);
1235  }
1236 #endif
1237 
1238  // else if(executionPolicy == ROOT::Fit::ExecutionPolicy::kMultiprocess){
1239  // ROOT::TProcessExecutor pool;
1240  // g = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, n), redFunction);
1241  // }
1242  else {
1243  Error("FitUtil::EvaluateLogLGradient", "Execution policy unknown. Avalaible choices:\n "
1244  "ROOT::Fit::ExecutionPolicy::kSerial (default)\n "
1245  "ROOT::Fit::ExecutionPolicy::kMultithread (requires IMT)\n");
1246  }
1247 
1248 #ifndef R__USE_IMT
1249  // to fix compiler warning
1250  (void)nChunks;
1251 #endif
1252 
1253  // copy result
1254  std::copy(g.begin(), g.end(), grad);
1255 
1256 #ifdef DEBUG
1257  std::cout << "FitUtil.cxx : Final gradient ";
1258  for (unsigned int param = 0; param < npar; param++) {
1259  std::cout << " " << grad[param];
1260  }
1261  std::cout << "\n";
1262 #endif
1263 }
1264 //_________________________________________________________________________________________________
1265 // for binned log likelihood functions
1266 ////////////////////////////////////////////////////////////////////////////////
1267 /// evaluate the pdf (Poisson) contribution to the logl (return actually log of pdf)
1268 /// and its gradient
1269 
1270 double FitUtil::EvaluatePoissonBinPdf(const IModelFunction & func, const BinData & data, const double * p, unsigned int i, double * g ) {
1271  double y = 0;
1272  const double * x1 = data.GetPoint(i,y);
1273 
1274  const DataOptions & fitOpt = data.Opt();
1275  bool useBinIntegral = fitOpt.fIntegral && data.HasBinEdges();
1276  bool useBinVolume = (fitOpt.fBinVolume && data.HasBinEdges());
1277 
1278  IntegralEvaluator<> igEval( func, p, useBinIntegral);
1279  double fval = 0;
1280  const double * x2 = 0;
1281  // calculate the bin volume
1282  double binVolume = 1;
1283  std::vector<double> xc;
1284  if (useBinVolume) {
1285  unsigned int ndim = data.NDim();
1286  xc.resize(ndim);
1287  x2 = data.BinUpEdge(i);
1288  for (unsigned int j = 0; j < ndim; ++j) {
1289  binVolume *= std::abs( x2[j]-x1[j] );
1290  xc[j] = 0.5*(x2[j]+ x1[j]);
1291  }
1292  // normalize the bin volume using a reference value
1293  binVolume /= data.RefVolume();
1294  }
1295 
1296  const double * x = (useBinVolume) ? &xc.front() : x1;
1297 
1298  if (!useBinIntegral ) {
1299  fval = func ( x, p );
1300  }
1301  else {
1302  // calculate integral normalized (divided by bin volume)
1303  x2 = data.BinUpEdge(i);
1304  fval = igEval( x1, x2 ) ;
1305  }
1306  if (useBinVolume) fval *= binVolume;
1307 
1308  // logPdf for Poisson: ignore constant term depending on N
1309  fval = std::max(fval, 0.0); // avoid negative or too small values
1310  double logPdf = - fval;
1311  if (y > 0.0) {
1312  // include also constants due to saturate model (see Baker-Cousins paper)
1313  logPdf += y * ROOT::Math::Util::EvalLog( fval / y) + y;
1314  }
1315  // need to return the pdf contribution (not the log)
1316 
1317  //double pdfval = std::exp(logPdf);
1318 
1319  //if (g == 0) return pdfval;
1320  if (g == 0) return logPdf;
1321 
1322  unsigned int npar = func.NPar();
1323  const IGradModelFunction * gfunc = dynamic_cast<const IGradModelFunction *>( &func);
1324 
1325  // gradient calculation
1326  if (gfunc != 0) {
1327  //case function provides gradient
1328  if (!useBinIntegral )
1329  gfunc->ParameterGradient( x , p, g );
1330  else {
1331  // needs to calculate the integral for each partial derivative
1332  CalculateGradientIntegral( *gfunc, x1, x2, p, g);
1333  }
1334 
1335  }
1336  else {
1337  SimpleGradientCalculator gc(func.NPar(), func);
1338  if (!useBinIntegral )
1339  gc.ParameterGradient(x, p, fval, g);
1340  else {
1341  // needs to calculate the integral for each partial derivative
1342  CalculateGradientIntegral( gc, x1, x2, p, g);
1343  }
1344  }
1345  // correct g[] do be derivative of poisson term
1346  for (unsigned int k = 0; k < npar; ++k) {
1347  // apply bin volume correction
1348  if (useBinVolume) g[k] *= binVolume;
1349 
1350  // correct for Poisson term
1351  if ( fval > 0)
1352  g[k] *= ( y/fval - 1.) ;//* pdfval;
1353  else if (y > 0) {
1354  const double kdmax1 = std::sqrt( std::numeric_limits<double>::max() );
1355  g[k] *= kdmax1;
1356  }
1357  else // y == 0 cannot have negative y
1358  g[k] *= -1;
1359  }
1360 
1361 
1362 #ifdef DEBUG
1363  std::cout << "x = " << x[0] << " logPdf = " << logPdf << " grad";
1364  for (unsigned int ipar = 0; ipar < npar; ++ipar)
1365  std::cout << g[ipar] << "\t";
1366  std::cout << std::endl;
1367 #endif
1368 
1369 // return pdfval;
1370  return logPdf;
1371 }
1372 
1373 double FitUtil::EvaluatePoissonLogL(const IModelFunction &func, const BinData &data, const double *p, int iWeight,
1374  bool extended, unsigned int &nPoints, ROOT::Fit::ExecutionPolicy executionPolicy,
1375  unsigned nChunks)
1376 {
1377  // evaluate the Poisson Log Likelihood
1378  // for binned likelihood fits
1379  // this is Sum ( f(x_i) - y_i * log( f (x_i) ) )
1380  // add as well constant term for saturated model to make it like a Chi2/2
1381  // by default is etended. If extended is false the fit is not extended and
1382  // the global poisson term is removed (i.e is a binomial fit)
1383  // (remember that in this case one needs to have a function with a fixed normalization
1384  // like in a non extended unbinned fit)
1385  //
1386  // if use Weight use a weighted dataset
1387  // iWeight = 1 ==> logL = Sum( w f(x_i) )
1388  // case of iWeight==1 is actually identical to weight==0
1389  // iWeight = 2 ==> logL = Sum( w*w * f(x_i) )
1390  //
1391  // nPoints returns the points where bin content is not zero
1392 
1393 
1394  unsigned int n = data.Size();
1395 
1396 #ifdef USE_PARAMCACHE
1397  (const_cast<IModelFunction &>(func)).SetParameters(p);
1398 #endif
1399 
1400  nPoints = data.Size(); // npoints
1401 
1402 
1403  // get fit option and check case of using integral of bins
1404  const DataOptions &fitOpt = data.Opt();
1405  bool useBinIntegral = fitOpt.fIntegral && data.HasBinEdges();
1406  bool useBinVolume = (fitOpt.fBinVolume && data.HasBinEdges());
1407  bool useW2 = (iWeight == 2);
1408 
1409  // normalize if needed by a reference volume value
1410  double wrefVolume = 1.0;
1411  if (useBinVolume) {
1412  if (fitOpt.fNormBinVolume) wrefVolume /= data.RefVolume();
1413  }
1414 
1415 //#define DEBUG
1416 #ifdef DEBUG
1417  std::cout << "Evaluate PoissonLogL for params = [ ";
1418  for (unsigned int j = 0; j < func.NPar(); ++j) std::cout << p[j] << " , ";
1419  std::cout << "] - data size = " << n << " useBinIntegral " << useBinIntegral << " useBinVolume "
1420  << useBinVolume << " useW2 " << useW2 << " wrefVolume = " << wrefVolume << std::endl;
1421 #endif
1422 
1423 
1424  ROOT::Math::IntegrationOneDim::Type igType = ROOT::Math::IntegrationOneDim::kDEFAULT;
1425  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
1426  // do not use GSL integrator which is not thread safe
1427  igType = ROOT::Math::IntegrationOneDim::kGAUSS;
1428  }
1429 #ifdef USE_PARAMCACHE
1430  IntegralEvaluator<> igEval(func, 0, useBinIntegral, igType);
1431 #else
1432  IntegralEvaluator<> igEval(func, p, useBinIntegral, igType);
1433 #endif
1434 
1435  auto mapFunction = [&](const unsigned i) {
1436  auto x1 = data.GetCoordComponent(i, 0);
1437  auto y = *data.ValuePtr(i);
1438 
1439  const double *x = nullptr;
1440  std::vector<double> xc;
1441  double fval = 0;
1442  double binVolume = 1.0;
1443 
1444  if (useBinVolume) {
1445  unsigned int ndim = data.NDim();
1446  const double *x2 = data.BinUpEdge(i);
1447  xc.resize(data.NDim());
1448  for (unsigned int j = 0; j < ndim; ++j) {
1449  auto xx = *data.GetCoordComponent(i, j);
1450  binVolume *= std::abs(x2[j] - xx);
1451  xc[j] = 0.5 * (x2[j] + xx);
1452  }
1453  x = xc.data();
1454  // normalize the bin volume using a reference value
1455  binVolume *= wrefVolume;
1456  } else if (data.NDim() > 1) {
1457  xc.resize(data.NDim());
1458  xc[0] = *x1;
1459  for (unsigned int j = 1; j < data.NDim(); ++j) {
1460  xc[j] = *data.GetCoordComponent(i, j);
1461  }
1462  x = xc.data();
1463  } else {
1464  x = x1;
1465  }
1466 
1467  if (!useBinIntegral) {
1468 #ifdef USE_PARAMCACHE
1469  fval = func(x);
1470 #else
1471  fval = func(x, p);
1472 #endif
1473  } else {
1474  // calculate integral (normalized by bin volume)
1475  // need to set function and parameters here in case loop is parallelized
1476  fval = igEval(x, data.BinUpEdge(i));
1477  }
1478  if (useBinVolume) fval *= binVolume;
1479 
1480 
1481 
1482 #ifdef DEBUG
1483  int NSAMPLE = 100;
1484  if (i % NSAMPLE == 0) {
1485  std::cout << "evt " << i << " x = [ ";
1486  for (unsigned int j = 0; j < func.NDim(); ++j) std::cout << x[j] << " , ";
1487  std::cout << "] ";
1488  if (fitOpt.fIntegral) {
1489  std::cout << "x2 = [ ";
1490  for (unsigned int j = 0; j < func.NDim(); ++j) std::cout << data.BinUpEdge(i)[j] << " , ";
1491  std::cout << "] ";
1492  }
1493  std::cout << " y = " << y << " fval = " << fval << std::endl;
1494  }
1495 #endif
1496 
1497 
1498  // EvalLog protects against 0 values of fval but don't want to add in the -log sum
1499  // negative values of fval
1500  fval = std::max(fval, 0.0);
1501 
1502  double nloglike = 0; // negative loglikelihood
1503  if (useW2) {
1504  // apply weight correction . Effective weight is error^2/ y
1505  // and expected events in bins is fval/weight
1506  // can apply correction only when y is not zero otherwise weight is undefined
1507  // (in case of weighted likelihood I don't care about the constant term due to
1508  // the saturated model)
1509 
1510  // use for the empty bins the global weight
1511  double weight = 1.0;
1512  if (y != 0) {
1513  double error = data.Error(i);
1514  weight = (error * error) / y; // this is the bin effective weight
1515  nloglike += weight * y * ( ROOT::Math::Util::EvalLog(y) - ROOT::Math::Util::EvalLog(fval) );
1516  }
1517  else {
1518  // for empty bin use the average weight computed from the total data weight
1519  weight = data.SumOfError2()/ data.SumOfContent();
1520  }
1521  if (extended) {
1522  nloglike += weight * ( fval - y);
1523  }
1524 
1525  } else {
1526  // standard case no weights or iWeight=1
1527  // this is needed for Poisson likelihood (which are extened and not for multinomial)
1528  // the formula below include constant term due to likelihood of saturated model (f(x) = y)
1529  // (same formula as in Baker-Cousins paper, page 439 except a factor of 2
1530  if (extended) nloglike = fval - y;
1531 
1532  if (y > 0) {
1533  nloglike += y * (ROOT::Math::Util::EvalLog(y) - ROOT::Math::Util::EvalLog(fval));
1534  }
1535  }
1536 #ifdef DEBUG
1537  {
1538  R__LOCKGUARD(gROOTMutex);
1539  std::cout << " nll = " << nloglike << std::endl;
1540  }
1541 #endif
1542  return nloglike;
1543  };
1544 
1545 #ifdef R__USE_IMT
1546  auto redFunction = [](const std::vector<double> &objs) {
1547  return std::accumulate(objs.begin(), objs.end(), double{});
1548  };
1549 #else
1550  (void)nChunks;
1551 
1552  // If IMT is disabled, force the execution policy to the serial case
1553  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
1554  Warning("FitUtil::EvaluatePoissonLogL", "Multithread execution policy requires IMT, which is disabled. Changing "
1555  "to ROOT::Fit::ExecutionPolicy::kSerial.");
1556  executionPolicy = ROOT::Fit::ExecutionPolicy::kSerial;
1557  }
1558 #endif
1559 
1560  double res{};
1561  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kSerial) {
1562  for (unsigned int i = 0; i < n; ++i) {
1563  res += mapFunction(i);
1564  }
1565 #ifdef R__USE_IMT
1566  } else if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
1567  ROOT::TThreadExecutor pool;
1568  auto chunks = nChunks != 0 ? nChunks : setAutomaticChunking(data.Size());
1569  res = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, n), redFunction, chunks);
1570 #endif
1571  // } else if(executionPolicy == ROOT::Fit::kMultitProcess){
1572  // ROOT::TProcessExecutor pool;
1573  // res = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, n), redFunction);
1574  } else {
1575  Error("FitUtil::EvaluatePoissonLogL",
1576  "Execution policy unknown. Avalaible choices:\n ROOT::Fit::ExecutionPolicy::kSerial (default)\n ROOT::Fit::ExecutionPolicy::kMultithread (requires IMT)\n");
1577  }
1578 
1579 #ifdef DEBUG
1580  std::cout << "Loglikelihood = " << res << std::endl;
1581 #endif
1582 
1583  return res;
1584 }
1585 
1586 void FitUtil::EvaluatePoissonLogLGradient(const IModelFunction &f, const BinData &data, const double *p, double *grad,
1587  unsigned int &, ROOT::Fit::ExecutionPolicy executionPolicy, unsigned nChunks)
1588 {
1589  // evaluate the gradient of the Poisson log likelihood function
1590 
1591  const IGradModelFunction *fg = dynamic_cast<const IGradModelFunction *>(&f);
1592  assert(fg != nullptr); // must be called by a grad function
1593 
1594  const IGradModelFunction &func = *fg;
1595 
1596 #ifdef USE_PARAMCACHE
1597  (const_cast<IGradModelFunction &>(func)).SetParameters(p);
1598 #endif
1599 
1600  const DataOptions &fitOpt = data.Opt();
1601  bool useBinIntegral = fitOpt.fIntegral && data.HasBinEdges();
1602  bool useBinVolume = (fitOpt.fBinVolume && data.HasBinEdges());
1603 
1604  double wrefVolume = 1.0;
1605  if (useBinVolume && fitOpt.fNormBinVolume)
1606  wrefVolume /= data.RefVolume();
1607 
1608  ROOT::Math::IntegrationOneDim::Type igType = ROOT::Math::IntegrationOneDim::kDEFAULT;
1609  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
1610  // do not use GSL integrator which is not thread safe
1611  igType = ROOT::Math::IntegrationOneDim::kGAUSS;
1612  }
1613 
1614  IntegralEvaluator<> igEval(func, p, useBinIntegral, igType);
1615 
1616  unsigned int npar = func.NPar();
1617  unsigned initialNPoints = data.Size();
1618 
1619  auto mapFunction = [&](const unsigned int i) {
1620  // set all vector values to zero
1621  std::vector<double> gradFunc(npar);
1622  std::vector<double> pointContribution(npar);
1623 
1624  const auto x1 = data.GetCoordComponent(i, 0);
1625  const auto y = data.Value(i);
1626  auto invError = data.Error(i);
1627 
1628  invError = (invError != 0.0) ? 1.0 / invError : 1;
1629 
1630  double fval = 0;
1631 
1632  const double *x = nullptr;
1633  std::vector<double> xc;
1634 
1635  unsigned ndim = data.NDim();
1636  double binVolume = 1.0;
1637  if (useBinVolume) {
1638  const double *x2 = data.BinUpEdge(i);
1639 
1640  xc.resize(ndim);
1641  for (unsigned int j = 0; j < ndim; ++j) {
1642  auto x1_j = *data.GetCoordComponent(i, j);
1643  binVolume *= std::abs(x2[j] - x1_j);
1644  xc[j] = 0.5 * (x2[j] + x1_j);
1645  }
1646 
1647  x = xc.data();
1648 
1649  // normalize the bin volume using a reference value
1650  binVolume *= wrefVolume;
1651  } else if (ndim > 1) {
1652  xc.resize(ndim);
1653  xc[0] = *x1;
1654  for (unsigned int j = 1; j < ndim; ++j)
1655  xc[j] = *data.GetCoordComponent(i, j);
1656  x = xc.data();
1657  } else {
1658  x = x1;
1659  }
1660 
1661  if (!useBinIntegral) {
1662  fval = func(x, p);
1663  func.ParameterGradient(x, p, &gradFunc[0]);
1664  } else {
1665  // calculate integral (normalized by bin volume)
1666  // need to set function and parameters here in case loop is parallelized
1667  auto x2 = data.BinUpEdge(i);
1668  fval = igEval(x, x2);
1669  CalculateGradientIntegral(func, x, x2, p, &gradFunc[0]);
1670  }
1671  if (useBinVolume)
1672  fval *= binVolume;
1673 
1674 #ifdef DEBUG
1675  {
1676  R__LOCKGUARD(gROOTMutex);
1677  if (i < 5 || (i > data.Size()-5) ) {
1678  if (data.NDim() > 1) std::cout << i << " x " << x[0] << " y " << x[1] << " func " << fval
1679  << " gradient " << gradFunc[0] << " " << gradFunc[1] << " " << gradFunc[3] << std::endl;
1680  else std::cout << i << " x " << x[0] << " gradient " << gradFunc[0] << " " << gradFunc[1] << " " << gradFunc[3] << std::endl;
1681  }
1682  }
1683 #endif
1684 
1685  // correct the gradient
1686  for (unsigned int ipar = 0; ipar < npar; ++ipar) {
1687 
1688  // correct gradient for bin volumes
1689  if (useBinVolume)
1690  gradFunc[ipar] *= binVolume;
1691 
1692  // df/dp * (1. - y/f )
1693  if (fval > 0)
1694  pointContribution[ipar] = gradFunc[ipar] * (1. - y / fval);
1695  else if (gradFunc[ipar] != 0) {
1696  const double kdmax1 = std::sqrt(std::numeric_limits<double>::max());
1697  const double kdmax2 = std::numeric_limits<double>::max() / (4 * initialNPoints);
1698  double gg = kdmax1 * gradFunc[ipar];
1699  if (gg > 0)
1700  gg = std::min(gg, kdmax2);
1701  else
1702  gg = std::max(gg, -kdmax2);
1703  pointContribution[ipar] = -gg;
1704  }
1705  }
1706 
1707 
1708  return pointContribution;
1709  };
1710 
1711  // Vertically reduce the set of vectors by summing its equally-indexed components
1712  auto redFunction = [&](const std::vector<std::vector<double>> &pointContributions) {
1713  std::vector<double> result(npar);
1714 
1715  for (auto const &pointContribution : pointContributions) {
1716  for (unsigned int parameterIndex = 0; parameterIndex < npar; parameterIndex++)
1717  result[parameterIndex] += pointContribution[parameterIndex];
1718  }
1719 
1720  return result;
1721  };
1722 
1723  std::vector<double> g(npar);
1724 
1725 #ifndef R__USE_IMT
1726  // If IMT is disabled, force the execution policy to the serial case
1727  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
1728  Warning("FitUtil::EvaluatePoissonLogLGradient",
1729  "Multithread execution policy requires IMT, which is disabled. Changing "
1730  "to ROOT::Fit::ExecutionPolicy::kSerial.");
1731  executionPolicy = ROOT::Fit::ExecutionPolicy::kSerial;
1732  }
1733 #endif
1734 
1735  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kSerial) {
1736  std::vector<std::vector<double>> allGradients(initialNPoints);
1737  for (unsigned int i = 0; i < initialNPoints; ++i) {
1738  allGradients[i] = mapFunction(i);
1739  }
1740  g = redFunction(allGradients);
1741  }
1742 #ifdef R__USE_IMT
1743  else if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
1744  ROOT::TThreadExecutor pool;
1745  auto chunks = nChunks != 0 ? nChunks : setAutomaticChunking(initialNPoints);
1746  g = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, initialNPoints), redFunction, chunks);
1747  }
1748 #endif
1749 
1750  // else if(executionPolicy == ROOT::Fit::kMultiprocess){
1751  // ROOT::TProcessExecutor pool;
1752  // g = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, n), redFunction);
1753  // }
1754  else {
1755  Error("FitUtil::EvaluatePoissonLogLGradient",
1756  "Execution policy unknown. Avalaible choices:\n 0: Serial (default)\n 1: MultiThread (requires IMT)\n");
1757  }
1758 
1759 #ifndef R__USE_IMT
1760  //to fix compiler warning
1761  (void)nChunks;
1762 #endif
1763 
1764  // copy result
1765  std::copy(g.begin(), g.end(), grad);
1766 
1767 #ifdef DEBUG
1768  std::cout << "***** Final gradient : ";
1769  for (unsigned int ii = 0; ii< npar; ++ii) std::cout << grad[ii] << " ";
1770  std::cout << "\n";
1771 #endif
1772 
1773 }
1774 
1775 
1776 unsigned FitUtil::setAutomaticChunking(unsigned nEvents){
1777  auto ncpu = ROOT::GetImplicitMTPoolSize();
1778  if (nEvents/ncpu < 1000) return ncpu;
1779  return nEvents/1000;
1780  //return ((nEvents/ncpu + 1) % 1000) *40 ; //arbitrary formula
1781 }
1782 
1783 }
1784 
1785 } // end namespace ROOT