Logo ROOT   master
Reference Guide
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"
22 #include "Math/WrappedFunction.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 
261  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
262  // do not use GSL integrator which is not thread safe
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  xc.resize(data.NDim());
295  for (unsigned int j = 0; j < ndim; ++j) {
296  double xx = *data.GetCoordComponent(i, j);
297  double x2 = data.GetBinUpEdgeComponent(i, j);
298  binVolume *= std::abs(x2 - xx);
299  xc[j] = 0.5*(x2 + 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  // multi-dim case (no bin volume)
306  xc.resize(data.NDim());
307  xc[0] = *x1;
308  for (unsigned int j = 1; j < data.NDim(); ++j)
309  xc[j] = *data.GetCoordComponent(i, j);
310  x = xc.data();
311  } else {
312  x = x1;
313  }
314 
315 
316  if (!useBinIntegral) {
317 #ifdef USE_PARAMCACHE
318  fval = func ( x );
319 #else
320  fval = func ( x, p );
321 #endif
322  }
323  else {
324  // calculate integral normalized by bin volume
325  // need to set function and parameters here in case loop is parallelized
326  std::vector<double> x2(data.NDim());
327  data.GetBinUpEdgeCoordinates(i, x2.data());
328  fval = igEval(x, x2.data());
329  }
330  // normalize result if requested according to bin volume
331  if (useBinVolume) fval *= binVolume;
332 
333  // expected errors
334  if (useExpErrors) {
335  double invWeight = 1.0;
336  if (isWeighted) {
337  // we need first to check if a weight factor needs to be applied
338  // weight = sumw2/sumw = error**2/content
339  //invWeight = y * invError * invError;
340  // we use always the global weight and not the observed one in the bin
341  // for empty bins use global weight (if it is weighted data.SumError2() is not zero)
342  invWeight = data.SumOfContent()/ data.SumOfError2();
343  //if (invError > 0) invWeight = y * invError * invError;
344  }
345 
346  // if (invError == 0) invWeight = (data.SumOfError2() > 0) ? data.SumOfContent()/ data.SumOfError2() : 1.0;
347  // compute expected error as f(x) / weight
348  double invError2 = (fval > 0) ? invWeight / fval : 0.0;
349  invError = std::sqrt(invError2);
350  //std::cout << "using Pearson chi2 " << x[0] << " " << 1./invError2 << " " << fval << std::endl;
351  }
352 
353 //#define DEBUG
354 #ifdef DEBUG
355  std::cout << x[0] << " " << y << " " << 1./invError << " params : ";
356  for (unsigned int ipar = 0; ipar < func.NPar(); ++ipar)
357  std::cout << p[ipar] << "\t";
358  std::cout << "\tfval = " << fval << " bin volume " << binVolume << " ref " << wrefVolume << std::endl;
359 #endif
360 //#undef DEBUG
361 
362  if (invError > 0) {
363 
364  double tmp = ( y -fval )* invError;
365  double resval = tmp * tmp;
366 
367 
368  // avoid inifinity or nan in chi2 values due to wrong function values
369  if ( resval < maxResValue )
370  chi2 += resval;
371  else {
372  //nRejected++;
373  chi2 += maxResValue;
374  }
375  }
376  return chi2;
377  };
378 
379 #ifdef R__USE_IMT
380  auto redFunction = [](const std::vector<double> & objs){
381  return std::accumulate(objs.begin(), objs.end(), double{});
382  };
383 #else
384  (void)nChunks;
385 
386  // If IMT is disabled, force the execution policy to the serial case
387  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
388  Warning("FitUtil::EvaluateChi2", "Multithread execution policy requires IMT, which is disabled. Changing "
389  "to ROOT::Fit::ExecutionPolicy::kSerial.");
390  executionPolicy = ROOT::Fit::ExecutionPolicy::kSerial;
391  }
392 #endif
393 
394  double res{};
395  if(executionPolicy == ROOT::Fit::ExecutionPolicy::kSerial){
396  for (unsigned int i=0; i<n; ++i) {
397  res += mapFunction(i);
398  }
399 #ifdef R__USE_IMT
400  } else if(executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
402  auto chunks = nChunks !=0? nChunks: setAutomaticChunking(data.Size());
403  res = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, n), redFunction, chunks);
404 #endif
405 // } else if(executionPolicy == ROOT::Fit::kMultitProcess){
406  // ROOT::TProcessExecutor pool;
407  // res = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, n), redFunction);
408  } else{
409  Error("FitUtil::EvaluateChi2","Execution policy unknown. Avalaible choices:\n ROOT::Fit::ExecutionPolicy::kSerial (default)\n ROOT::Fit::ExecutionPolicy::kMultithread (requires IMT)\n");
410  }
411 
412  return res;
413 }
414 
415 
416 //___________________________________________________________________________________________________________________________
417 
418 double FitUtil::EvaluateChi2Effective(const IModelFunction & func, const BinData & data, const double * p, unsigned int & nPoints) {
419  // evaluate the chi2 given a function reference , the data and returns the value and also in nPoints
420  // the actual number of used points
421  // method using the error in the coordinates
422  // integral of bin does not make sense in this case
423 
424  unsigned int n = data.Size();
425 
426 #ifdef DEBUG
427  std::cout << "\n\nFit data size = " << n << std::endl;
428  std::cout << "evaluate effective chi2 using function " << &func << " " << p << std::endl;
429 #endif
430 
431  assert(data.HaveCoordErrors() || data.HaveAsymErrors());
432 
433  double chi2 = 0;
434  //int nRejected = 0;
435 
436 
437  //func.SetParameters(p);
438 
439  unsigned int ndim = func.NDim();
440 
441  // use Richardson derivator
443 
444  double maxResValue = std::numeric_limits<double>::max() /n;
445 
446 
447 
448  for (unsigned int i = 0; i < n; ++ i) {
449 
450 
451  double y = 0;
452  const double * x = data.GetPoint(i,y);
453 
454  double fval = func( x, p );
455 
456  double delta_y_func = y - fval;
457 
458 
459  double ey = 0;
460  const double * ex = 0;
461  if (!data.HaveAsymErrors() )
462  ex = data.GetPointError(i, ey);
463  else {
464  double eylow, eyhigh = 0;
465  ex = data.GetPointError(i, eylow, eyhigh);
466  if ( delta_y_func < 0)
467  ey = eyhigh; // function is higher than points
468  else
469  ey = eylow;
470  }
471  double e2 = ey * ey;
472  // before calculating the gradient check that all error in x are not zero
473  unsigned int j = 0;
474  while ( j < ndim && ex[j] == 0.) { j++; }
475  // if j is less ndim some elements are not zero
476  if (j < ndim) {
477  // need an adapter from a multi-dim function to a one-dimensional
479  // select optimal step size (use 10--2 by default as was done in TF1:
480  double kEps = 0.01;
481  double kPrecision = 1.E-8;
482  for (unsigned int icoord = 0; icoord < ndim; ++icoord) {
483  // calculate derivative for each coordinate
484  if (ex[icoord] > 0) {
485  //gradCalc.Gradient(x, p, fval, &grad[0]);
486  f1D.SetCoord(icoord);
487  // optimal spep size (take ex[] as scale for the points and 1% of it
488  double x0= x[icoord];
489  double h = std::max( kEps* std::abs(ex[icoord]), 8.0*kPrecision*(std::abs(x0) + kPrecision) );
490  double deriv = derivator.Derivative1(f1D, x[icoord], h);
491  double edx = ex[icoord] * deriv;
492  e2 += edx * edx;
493 #ifdef DEBUG
494  std::cout << "error for coord " << icoord << " = " << ex[icoord] << " deriv " << deriv << std::endl;
495 #endif
496  }
497  }
498  }
499  double w2 = (e2 > 0) ? 1.0/e2 : 0;
500  double resval = w2 * ( y - fval ) * ( y - fval);
501 
502 #ifdef DEBUG
503  std::cout << x[0] << " " << y << " ex " << ex[0] << " ey " << ey << " params : ";
504  for (unsigned int ipar = 0; ipar < func.NPar(); ++ipar)
505  std::cout << p[ipar] << "\t";
506  std::cout << "\tfval = " << fval << "\tresval = " << resval << std::endl;
507 #endif
508 
509  // avoid (infinity and nan ) in the chi2 sum
510  // eventually add possibility of excluding some points (like singularity)
511  if ( resval < maxResValue )
512  chi2 += resval;
513  else
514  chi2 += maxResValue;
515  //nRejected++;
516 
517  }
518 
519  // reset the number of fitting data points
520  nPoints = n; // no points are rejected
521  //if (nRejected != 0) nPoints = n - nRejected;
522 
523 #ifdef DEBUG
524  std::cout << "chi2 = " << chi2 << " n = " << nPoints << std::endl;
525 #endif
526 
527  return chi2;
528 
529 }
530 
531 
532 ////////////////////////////////////////////////////////////////////////////////
533 /// evaluate the chi2 contribution (residual term) only for data with no coord-errors
534 /// This function is used in the specialized least square algorithms like FUMILI or L.M.
535 /// if we have error on the coordinates the method is not yet implemented
536 /// integral option is also not yet implemented
537 /// one can use in that case normal chi2 method
538 
539 double FitUtil::EvaluateChi2Residual(const IModelFunction & func, const BinData & data, const double * p, unsigned int i, double * g) {
540  if (data.GetErrorType() == BinData::kCoordError && data.Opt().fCoordErrors ) {
541  MATH_ERROR_MSG("FitUtil::EvaluateChi2Residual","Error on the coordinates are not used in calculating Chi2 residual");
542  return 0; // it will assert otherwise later in GetPoint
543  }
544 
545 
546  //func.SetParameters(p);
547 
548  double y, invError = 0;
549 
550  const DataOptions & fitOpt = data.Opt();
551  bool useBinIntegral = fitOpt.fIntegral && data.HasBinEdges();
552  bool useBinVolume = (fitOpt.fBinVolume && data.HasBinEdges());
553  bool useExpErrors = (fitOpt.fExpErrors);
554 
555  const double * x1 = data.GetPoint(i,y, invError);
556 
557  IntegralEvaluator<> igEval( func, p, useBinIntegral);
558  double fval = 0;
559  unsigned int ndim = data.NDim();
560  double binVolume = 1.0;
561  const double * x2 = 0;
562  if (useBinVolume || useBinIntegral) x2 = data.BinUpEdge(i);
563 
564  double * xc = 0;
565 
566  if (useBinVolume) {
567  xc = new double[ndim];
568  for (unsigned int j = 0; j < ndim; ++j) {
569  binVolume *= std::abs( x2[j]-x1[j] );
570  xc[j] = 0.5*(x2[j]+ x1[j]);
571  }
572  // normalize the bin volume using a reference value
573  binVolume /= data.RefVolume();
574  }
575 
576  const double * x = (useBinVolume) ? xc : x1;
577 
578  if (!useBinIntegral) {
579  fval = func ( x, p );
580  }
581  else {
582  // calculate integral (normalized by bin volume)
583  // need to set function and parameters here in case loop is parallelized
584  fval = igEval( x1, x2) ;
585  }
586  // normalize result if requested according to bin volume
587  if (useBinVolume) fval *= binVolume;
588 
589  // expected errors
590  if (useExpErrors) {
591  // we need first to check if a weight factor needs to be applied
592  // weight = sumw2/sumw = error**2/content
593  //NOTE: assume histogram is not weighted
594  // don't know how to do with bins with weight = 0
595  //double invWeight = y * invError * invError;
596  // if (invError == 0) invWeight = (data.SumOfError2() > 0) ? data.SumOfContent()/ data.SumOfError2() : 1.0;
597  // compute expected error as f(x) / weight
598  double invError2 = (fval > 0) ? 1.0 / fval : 0.0;
599  invError = std::sqrt(invError2);
600  }
601 
602 
603  double resval = ( y -fval )* invError;
604 
605  // avoid infinities or nan in resval
606  resval = CorrectValue(resval);
607 
608  // estimate gradient
609  if (g != 0) {
610 
611  unsigned int npar = func.NPar();
612  const IGradModelFunction * gfunc = dynamic_cast<const IGradModelFunction *>( &func);
613 
614  if (gfunc != 0) {
615  //case function provides gradient
616  if (!useBinIntegral ) {
617  gfunc->ParameterGradient( x , p, g );
618  }
619  else {
620  // needs to calculate the integral for each partial derivative
621  CalculateGradientIntegral( *gfunc, x1, x2, p, g);
622  }
623  }
624  else {
625  SimpleGradientCalculator gc( npar, func);
626  if (!useBinIntegral )
627  gc.ParameterGradient(x, p, fval, g);
628  else {
629  // needs to calculate the integral for each partial derivative
630  CalculateGradientIntegral( gc, x1, x2, p, g);
631  }
632  }
633  // mutiply by - 1 * weight
634  for (unsigned int k = 0; k < npar; ++k) {
635  g[k] *= - invError;
636  if (useBinVolume) g[k] *= binVolume;
637  }
638  }
639 
640  if (useBinVolume) delete [] xc;
641 
642  return resval;
643 
644 }
645 
646 void FitUtil::EvaluateChi2Gradient(const IModelFunction &f, const BinData &data, const double *p, double *grad,
647  unsigned int &nPoints, ROOT::Fit::ExecutionPolicy executionPolicy, unsigned nChunks)
648 {
649  // evaluate the gradient of the chi2 function
650  // this function is used when the model function knows how to calculate the derivative and we can
651  // avoid that the minimizer re-computes them
652  //
653  // case of chi2 effective (errors on coordinate) is not supported
654 
655  if (data.HaveCoordErrors()) {
656  MATH_ERROR_MSG("FitUtil::EvaluateChi2Gradient",
657  "Error on the coordinates are not used in calculating Chi2 gradient");
658  return; // it will assert otherwise later in GetPoint
659  }
660 
661  const IGradModelFunction *fg = dynamic_cast<const IGradModelFunction *>(&f);
662  assert(fg != nullptr); // must be called by a gradient function
663 
664  const IGradModelFunction &func = *fg;
665 
666 #ifdef DEBUG
667  std::cout << "\n\nFit data size = " << n << std::endl;
668  std::cout << "evaluate chi2 using function gradient " << &func << " " << p << std::endl;
669 #endif
670 
671  const DataOptions &fitOpt = data.Opt();
672  bool useBinIntegral = fitOpt.fIntegral && data.HasBinEdges();
673  bool useBinVolume = (fitOpt.fBinVolume && data.HasBinEdges());
674 
675  double wrefVolume = 1.0;
676  if (useBinVolume) {
677  if (fitOpt.fNormBinVolume) wrefVolume /= data.RefVolume();
678  }
679 
681  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
682  // do not use GSL integrator which is not thread safe
684  }
685  IntegralEvaluator<> igEval(func, p, useBinIntegral,igType);
686 
687  unsigned int npar = func.NPar();
688  unsigned initialNPoints = data.Size();
689 
690  std::vector<bool> isPointRejected(initialNPoints);
691 
692  auto mapFunction = [&](const unsigned int i) {
693  // set all vector values to zero
694  std::vector<double> gradFunc(npar);
695  std::vector<double> pointContribution(npar);
696 
697  const auto x1 = data.GetCoordComponent(i, 0);
698  const auto y = data.Value(i);
699  auto invError = data.Error(i);
700 
701  invError = (invError != 0.0) ? 1.0 / invError : 1;
702 
703  double fval = 0;
704 
705  const double *x = nullptr;
706  std::vector<double> xc;
707 
708  unsigned int ndim = data.NDim();
709  double binVolume = 1;
710  if (useBinVolume) {
711  xc.resize(ndim);
712  for (unsigned int j = 0; j < ndim; ++j) {
713  double x1_j = *data.GetCoordComponent(i, j);
714  double x2_j = data.GetBinUpEdgeComponent(i, j);
715  binVolume *= std::abs(x2_j - x1_j);
716  xc[j] = 0.5 * (x2_j + x1_j);
717  }
718 
719  x = xc.data();
720 
721  // normalize the bin volume using a reference value
722  binVolume *= wrefVolume;
723  } else if (ndim > 1) {
724  xc.resize(ndim);
725  xc[0] = *x1;
726  for (unsigned int j = 1; j < ndim; ++j)
727  xc[j] = *data.GetCoordComponent(i, j);
728  x = xc.data();
729  } else {
730  x = x1;
731  }
732 
733  if (!useBinIntegral) {
734  fval = func(x, p);
735  func.ParameterGradient(x, p, &gradFunc[0]);
736  } else {
737  std::vector<double> x2(data.NDim());
738  data.GetBinUpEdgeCoordinates(i, x2.data());
739  // calculate normalized integral and gradient (divided by bin volume)
740  // need to set function and parameters here in case loop is parallelized
741  fval = igEval(x, x2.data());
742  CalculateGradientIntegral(func, x, x2.data(), p, &gradFunc[0]);
743  }
744  if (useBinVolume)
745  fval *= binVolume;
746 
747 #ifdef DEBUG
748  std::cout << x[0] << " " << y << " " << 1. / invError << " params : ";
749  for (unsigned int ipar = 0; ipar < npar; ++ipar)
750  std::cout << p[ipar] << "\t";
751  std::cout << "\tfval = " << fval << std::endl;
752 #endif
753  if (!CheckInfNaNValue(fval)) {
754  isPointRejected[i] = true;
755  // Return a zero contribution to all partial derivatives on behalf of the current point
756  return pointContribution;
757  }
758 
759  // loop on the parameters
760  unsigned int ipar = 0;
761  for (; ipar < npar; ++ipar) {
762 
763  // correct gradient for bin volumes
764  if (useBinVolume)
765  gradFunc[ipar] *= binVolume;
766 
767  // avoid singularity in the function (infinity and nan ) in the chi2 sum
768  // eventually add possibility of excluding some points (like singularity)
769  double dfval = gradFunc[ipar];
770  if (!CheckInfNaNValue(dfval)) {
771  break; // exit loop on parameters
772  }
773 
774  // calculate derivative point contribution
775  pointContribution[ipar] = -2.0 * (y - fval) * invError * invError * gradFunc[ipar];
776  }
777 
778  if (ipar < npar) {
779  // case loop was broken for an overflow in the gradient calculation
780  isPointRejected[i] = true;
781  }
782 
783  return pointContribution;
784  };
785 
786  // Vertically reduce the set of vectors by summing its equally-indexed components
787  auto redFunction = [&](const std::vector<std::vector<double>> &pointContributions) {
788  std::vector<double> result(npar);
789 
790  for (auto const &pointContribution : pointContributions) {
791  for (unsigned int parameterIndex = 0; parameterIndex < npar; parameterIndex++)
792  result[parameterIndex] += pointContribution[parameterIndex];
793  }
794 
795  return result;
796  };
797 
798  std::vector<double> g(npar);
799 
800 #ifndef R__USE_IMT
801  // If IMT is disabled, force the execution policy to the serial case
802  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
803  Warning("FitUtil::EvaluateChi2Gradient", "Multithread execution policy requires IMT, which is disabled. Changing "
804  "to ROOT::Fit::ExecutionPolicy::kSerial.");
805  executionPolicy = ROOT::Fit::ExecutionPolicy::kSerial;
806  }
807 #endif
808 
809  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kSerial) {
810  std::vector<std::vector<double>> allGradients(initialNPoints);
811  for (unsigned int i = 0; i < initialNPoints; ++i) {
812  allGradients[i] = mapFunction(i);
813  }
814  g = redFunction(allGradients);
815  }
816 #ifdef R__USE_IMT
817  else if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
819  auto chunks = nChunks != 0 ? nChunks : setAutomaticChunking(initialNPoints);
820  g = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, initialNPoints), redFunction, chunks);
821  }
822 #endif
823  // else if(executionPolicy == ROOT::Fit::kMultiprocess){
824  // ROOT::TProcessExecutor pool;
825  // g = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, n), redFunction);
826  // }
827  else {
828  Error("FitUtil::EvaluateChi2Gradient",
829  "Execution policy unknown. Avalaible choices:\n 0: Serial (default)\n 1: MultiThread (requires IMT)\n");
830  }
831 
832 #ifndef R__USE_IMT
833  //to fix compiler warning
834  (void)nChunks;
835 #endif
836 
837  // correct the number of points
838  nPoints = initialNPoints;
839 
840  if (std::any_of(isPointRejected.begin(), isPointRejected.end(), [](bool point) { return point; })) {
841  unsigned nRejected = std::accumulate(isPointRejected.begin(), isPointRejected.end(), 0);
842  assert(nRejected <= initialNPoints);
843  nPoints = initialNPoints - nRejected;
844 
845  if (nPoints < npar)
846  MATH_ERROR_MSG("FitUtil::EvaluateChi2Gradient",
847  "Error - too many points rejected for overflow in gradient calculation");
848  }
849 
850  // copy result
851  std::copy(g.begin(), g.end(), grad);
852 }
853 
854 //______________________________________________________________________________________________________
855 //
856 // Log Likelihood functions
857 //_______________________________________________________________________________________________________
858 
859 // utility function used by the likelihoods
860 
861 // for LogLikelihood functions
862 
863 double FitUtil::EvaluatePdf(const IModelFunction & func, const UnBinData & data, const double * p, unsigned int i, double * g) {
864  // evaluate the pdf contribution to the generic logl function in case of bin data
865  // return actually the log of the pdf and its derivatives
866 
867 
868  //func.SetParameters(p);
869 
870 
871  const double * x = data.Coords(i);
872  double fval = func ( x, p );
873  double logPdf = ROOT::Math::Util::EvalLog(fval);
874  //return
875  if (g == 0) return logPdf;
876 
877  const IGradModelFunction * gfunc = dynamic_cast<const IGradModelFunction *>( &func);
878 
879  // gradient calculation
880  if (gfunc != 0) {
881  //case function provides gradient
882  gfunc->ParameterGradient( x , p, g );
883  }
884  else {
885  // estimate gradieant numerically with simple 2 point rule
886  // should probably calculate gradient of log(pdf) is more stable numerically
887  SimpleGradientCalculator gc(func.NPar(), func);
888  gc.ParameterGradient(x, p, fval, g );
889  }
890  // divide gradient by function value since returning the logs
891  for (unsigned int ipar = 0; ipar < func.NPar(); ++ipar) {
892  g[ipar] /= fval; // this should be checked against infinities
893  }
894 
895 #ifdef DEBUG
896  std::cout << x[i] << "\t";
897  std::cout << "\tpar = [ " << func.NPar() << " ] = ";
898  for (unsigned int ipar = 0; ipar < func.NPar(); ++ipar)
899  std::cout << p[ipar] << "\t";
900  std::cout << "\tfval = " << fval;
901  std::cout << "\tgrad = [ ";
902  for (unsigned int ipar = 0; ipar < func.NPar(); ++ipar)
903  std::cout << g[ipar] << "\t";
904  std::cout << " ] " << std::endl;
905 #endif
906 
907 
908  return logPdf;
909 }
910 
911 double FitUtil::EvaluateLogL(const IModelFunctionTempl<double> &func, const UnBinData &data, const double *p,
912  int iWeight, bool extended, unsigned int &nPoints,
913  ROOT::Fit::ExecutionPolicy executionPolicy, unsigned nChunks)
914 {
915  // evaluate the LogLikelihood
916 
917  unsigned int n = data.Size();
918 
919  //unsigned int nRejected = 0;
920 
921  bool normalizeFunc = false;
922 
923  // set parameters of the function to cache integral value
924 #ifdef USE_PARAMCACHE
925  (const_cast<IModelFunctionTempl<double> &>(func)).SetParameters(p);
926 #endif
927 
928  nPoints = data.Size(); // npoints
929 
930 #ifdef R__USE_IMT
931  // in case parameter needs to be propagated to user function use trick to set parameters by calling one time the function
932  // this will be done in sequential mode and parameters can be set in a thread safe manner
933  if (!normalizeFunc) {
934  if (data.NDim() == 1) {
935  const double * x = data.GetCoordComponent(0,0);
936  func( x, p);
937  }
938  else {
939  std::vector<double> x(data.NDim());
940  for (unsigned int j = 0; j < data.NDim(); ++j)
941  x[j] = *data.GetCoordComponent(0, j);
942  func( x.data(), p);
943  }
944  }
945 #endif
946 
947  double norm = 1.0;
948  if (normalizeFunc) {
949  // compute integral of the function
950  std::vector<double> xmin(data.NDim());
951  std::vector<double> xmax(data.NDim());
952  IntegralEvaluator<> igEval(func, p, true);
953  // compute integral in the ranges where is defined
954  if (data.Range().Size() > 0) {
955  norm = 0;
956  for (unsigned int ir = 0; ir < data.Range().Size(); ++ir) {
957  data.Range().GetRange(&xmin[0], &xmax[0], ir);
958  norm += igEval.Integral(xmin.data(), xmax.data());
959  }
960  } else {
961  // use (-inf +inf)
962  data.Range().GetRange(&xmin[0], &xmax[0]);
963  // check if funcition is zero at +- inf
964  if (func(xmin.data(), p) != 0 || func(xmax.data(), p) != 0) {
965  MATH_ERROR_MSG("FitUtil::EvaluateLogLikelihood",
966  "A range has not been set and the function is not zero at +/- inf");
967  return 0;
968  }
969  norm = igEval.Integral(&xmin[0], &xmax[0]);
970  }
971  }
972 
973  // needed to compue effective global weight in case of extended likelihood
974 
975  auto mapFunction = [&](const unsigned i) {
976  double W = 0;
977  double W2 = 0;
978  double fval = 0;
979 
980  if (data.NDim() > 1) {
981  std::vector<double> x(data.NDim());
982  for (unsigned int j = 0; j < data.NDim(); ++j)
983  x[j] = *data.GetCoordComponent(i, j);
984 #ifdef USE_PARAMCACHE
985  fval = func(x.data());
986 #else
987  fval = func(x.data(), p);
988 #endif
989 
990  // one -dim case
991  } else {
992  const auto x = data.GetCoordComponent(i, 0);
993 #ifdef USE_PARAMCACHE
994  fval = func(x);
995 #else
996  fval = func(x, p);
997 #endif
998  }
999 
1000  if (normalizeFunc)
1001  fval = fval * (1 / norm);
1002 
1003  // function EvalLog protects against negative or too small values of fval
1004  double logval = ROOT::Math::Util::EvalLog(fval);
1005  if (iWeight > 0) {
1006  double weight = data.Weight(i);
1007  logval *= weight;
1008  if (iWeight == 2) {
1009  logval *= weight; // use square of weights in likelihood
1010  if (!extended) {
1011  // needed sum of weights and sum of weight square if likelkihood is extended
1012  W = weight;
1013  W2 = weight * weight;
1014  }
1015  }
1016  }
1017  return LikelihoodAux<double>(logval, W, W2);
1018  };
1019 
1020 #ifdef R__USE_IMT
1021  // auto redFunction = [](const std::vector<LikelihoodAux<double>> & objs){
1022  // return std::accumulate(objs.begin(), objs.end(), LikelihoodAux<double>(0.0,0.0,0.0),
1023  // [](const LikelihoodAux<double> &l1, const LikelihoodAux<double> &l2){
1024  // return l1+l2;
1025  // });
1026  // };
1027  // do not use std::accumulate to be sure to maintain always the same order
1028  auto redFunction = [](const std::vector<LikelihoodAux<double>> & objs){
1029  auto l0 = LikelihoodAux<double>(0.0,0.0,0.0);
1030  for ( auto & l : objs ) {
1031  l0 = l0 + l;
1032  }
1033  return l0;
1034  };
1035 #else
1036  (void)nChunks;
1037 
1038  // If IMT is disabled, force the execution policy to the serial case
1039  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
1040  Warning("FitUtil::EvaluateLogL", "Multithread execution policy requires IMT, which is disabled. Changing "
1041  "to ROOT::Fit::ExecutionPolicy::kSerial.");
1042  executionPolicy = ROOT::Fit::ExecutionPolicy::kSerial;
1043  }
1044 #endif
1045 
1046  double logl{};
1047  double sumW{};
1048  double sumW2{};
1049  if(executionPolicy == ROOT::Fit::ExecutionPolicy::kSerial){
1050  for (unsigned int i=0; i<n; ++i) {
1051  auto resArray = mapFunction(i);
1052  logl+=resArray.logvalue;
1053  sumW+=resArray.weight;
1054  sumW2+=resArray.weight2;
1055  }
1056 #ifdef R__USE_IMT
1057  } else if(executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
1058  ROOT::TThreadExecutor pool;
1059  auto chunks = nChunks !=0? nChunks: setAutomaticChunking(data.Size());
1060  auto resArray = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, n), redFunction, chunks);
1061  logl=resArray.logvalue;
1062  sumW=resArray.weight;
1063  sumW2=resArray.weight2;
1064 #endif
1065 // } else if(executionPolicy == ROOT::Fit::kMultitProcess){
1066  // ROOT::TProcessExecutor pool;
1067  // res = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, n), redFunction);
1068  } else{
1069  Error("FitUtil::EvaluateLogL","Execution policy unknown. Avalaible choices:\n ROOT::Fit::ExecutionPolicy::kSerial (default)\n ROOT::Fit::ExecutionPolicy::kMultithread (requires IMT)\n");
1070  }
1071 
1072  if (extended) {
1073  // add Poisson extended term
1074  double extendedTerm = 0; // extended term in likelihood
1075  double nuTot = 0;
1076  // nuTot is integral of function in the range
1077  // if function has been normalized integral has been already computed
1078  if (!normalizeFunc) {
1079  IntegralEvaluator<> igEval( func, p, true);
1080  std::vector<double> xmin(data.NDim());
1081  std::vector<double> xmax(data.NDim());
1082 
1083  // compute integral in the ranges where is defined
1084  if (data.Range().Size() > 0 ) {
1085  nuTot = 0;
1086  for (unsigned int ir = 0; ir < data.Range().Size(); ++ir) {
1087  data.Range().GetRange(&xmin[0],&xmax[0],ir);
1088  nuTot += igEval.Integral(xmin.data(),xmax.data());
1089  }
1090  } else {
1091  // use (-inf +inf)
1092  data.Range().GetRange(&xmin[0],&xmax[0]);
1093  // check if funcition is zero at +- inf
1094  if (func(xmin.data(), p) != 0 || func(xmax.data(), p) != 0) {
1095  MATH_ERROR_MSG("FitUtil::EvaluateLogLikelihood","A range has not been set and the function is not zero at +/- inf");
1096  return 0;
1097  }
1098  nuTot = igEval.Integral(&xmin[0],&xmax[0]);
1099  }
1100 
1101  // force to be last parameter value
1102  //nutot = p[func.NDim()-1];
1103  if (iWeight != 2)
1104  extendedTerm = - nuTot; // no need to add in this case n log(nu) since is already computed before
1105  else {
1106  // case use weight square in likelihood : compute total effective weight = sw2/sw
1107  // ignore for the moment case when sumW is zero
1108  extendedTerm = - (sumW2 / sumW) * nuTot;
1109  }
1110 
1111  }
1112  else {
1113  nuTot = norm;
1114  extendedTerm = - nuTot + double(n) * ROOT::Math::Util::EvalLog( nuTot);
1115  // in case of weights need to use here sum of weights (to be done)
1116  }
1117  logl += extendedTerm;
1118 
1119  }
1120 
1121 #ifdef DEBUG
1122  std::cout << "Evaluated log L for parameters (";
1123  for (unsigned int ip = 0; ip < func.NPar(); ++ip)
1124  std::cout << " " << p[ip];
1125  std::cout << ") fval = " << -logl << std::endl;
1126 #endif
1127 
1128  return -logl;
1129 }
1130 
1131 void FitUtil::EvaluateLogLGradient(const IModelFunction &f, const UnBinData &data, const double *p, double *grad,
1132  unsigned int &, ROOT::Fit::ExecutionPolicy executionPolicy, unsigned nChunks)
1133 {
1134  // evaluate the gradient of the log likelihood function
1135 
1136  const IGradModelFunction *fg = dynamic_cast<const IGradModelFunction *>(&f);
1137  assert(fg != nullptr); // must be called by a grad function
1138 
1139  const IGradModelFunction &func = *fg;
1140 
1141  unsigned int npar = func.NPar();
1142  unsigned initialNPoints = data.Size();
1143 
1144  (const_cast<IGradModelFunction &>(func)).SetParameters(p);
1145 
1146 #ifdef DEBUG
1147  std::cout << "\n===> Evaluate Gradient for parameters ";
1148  for (unsigned int ip = 0; ip < npar; ++ip)
1149  std::cout << " " << p[ip];
1150  std::cout << "\n";
1151 #endif
1152 
1153  const double kdmax1 = std::sqrt(std::numeric_limits<double>::max());
1154  const double kdmax2 = std::numeric_limits<double>::max() / (4 * initialNPoints);
1155 
1156  auto mapFunction = [&](const unsigned int i) {
1157  std::vector<double> gradFunc(npar);
1158  std::vector<double> pointContribution(npar);
1159 
1160 
1161  const double * x = nullptr;
1162  std::vector<double> xc;
1163  if (data.NDim() > 1) {
1164  xc.resize(data.NDim() );
1165  for (unsigned int j = 0; j < data.NDim(); ++j)
1166  xc[j] = *data.GetCoordComponent(i, j);
1167  x = xc.data();
1168  } else {
1169  x = data.GetCoordComponent(i, 0);
1170  }
1171 
1172  double fval = func(x, p);
1173  func.ParameterGradient(x, p, &gradFunc[0]);
1174 
1175 #ifdef DEBUG
1176  {
1178  if (i < 5 || (i > data.Size()-5) ) {
1179  if (data.NDim() > 1) std::cout << i << " x " << x[0] << " y " << x[1] << " func " << fval
1180  << " gradient " << gradFunc[0] << " " << gradFunc[1] << " " << gradFunc[3] << std::endl;
1181  else std::cout << i << " x " << x[0] << " gradient " << gradFunc[0] << " " << gradFunc[1] << " " << gradFunc[3] << std::endl;
1182  }
1183  }
1184 #endif
1185 
1186  for (unsigned int kpar = 0; kpar < npar; ++kpar) {
1187  if (fval > 0)
1188  pointContribution[kpar] = -1. / fval * gradFunc[kpar];
1189  else if (gradFunc[kpar] != 0) {
1190  double gg = kdmax1 * gradFunc[kpar];
1191  if (gg > 0)
1192  gg = std::min(gg, kdmax2);
1193  else
1194  gg = std::max(gg, -kdmax2);
1195  pointContribution[kpar] = -gg;
1196  }
1197  // if func derivative is zero term is also zero so do not add in g[kpar]
1198  }
1199 
1200  return pointContribution;
1201  };
1202 
1203  // Vertically reduce the set of vectors by summing its equally-indexed components
1204  auto redFunction = [&](const std::vector<std::vector<double>> &pointContributions) {
1205  std::vector<double> result(npar);
1206 
1207  for (auto const &pointContribution : pointContributions) {
1208  for (unsigned int parameterIndex = 0; parameterIndex < npar; parameterIndex++)
1209  result[parameterIndex] += pointContribution[parameterIndex];
1210  }
1211 
1212  return result;
1213  };
1214 
1215  std::vector<double> g(npar);
1216 
1217 #ifndef R__USE_IMT
1218  // If IMT is disabled, force the execution policy to the serial case
1219  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
1220  Warning("FitUtil::EvaluateLogLGradient", "Multithread execution policy requires IMT, which is disabled. Changing "
1221  "to ROOT::Fit::ExecutionPolicy::kSerial.");
1222  executionPolicy = ROOT::Fit::ExecutionPolicy::kSerial;
1223  }
1224 #endif
1225 
1226  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kSerial) {
1227  std::vector<std::vector<double>> allGradients(initialNPoints);
1228  for (unsigned int i = 0; i < initialNPoints; ++i) {
1229  allGradients[i] = mapFunction(i);
1230  }
1231  g = redFunction(allGradients);
1232  }
1233 #ifdef R__USE_IMT
1234  else if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
1235  ROOT::TThreadExecutor pool;
1236  auto chunks = nChunks != 0 ? nChunks : setAutomaticChunking(initialNPoints);
1237  g = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, initialNPoints), redFunction, chunks);
1238  }
1239 #endif
1240 
1241  // else if(executionPolicy == ROOT::Fit::ExecutionPolicy::kMultiprocess){
1242  // ROOT::TProcessExecutor pool;
1243  // g = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, n), redFunction);
1244  // }
1245  else {
1246  Error("FitUtil::EvaluateLogLGradient", "Execution policy unknown. Avalaible choices:\n "
1247  "ROOT::Fit::ExecutionPolicy::kSerial (default)\n "
1248  "ROOT::Fit::ExecutionPolicy::kMultithread (requires IMT)\n");
1249  }
1250 
1251 #ifndef R__USE_IMT
1252  // to fix compiler warning
1253  (void)nChunks;
1254 #endif
1255 
1256  // copy result
1257  std::copy(g.begin(), g.end(), grad);
1258 
1259 #ifdef DEBUG
1260  std::cout << "FitUtil.cxx : Final gradient ";
1261  for (unsigned int param = 0; param < npar; param++) {
1262  std::cout << " " << grad[param];
1263  }
1264  std::cout << "\n";
1265 #endif
1266 }
1267 //_________________________________________________________________________________________________
1268 // for binned log likelihood functions
1269 ////////////////////////////////////////////////////////////////////////////////
1270 /// evaluate the pdf (Poisson) contribution to the logl (return actually log of pdf)
1271 /// and its gradient
1272 
1273 double FitUtil::EvaluatePoissonBinPdf(const IModelFunction & func, const BinData & data, const double * p, unsigned int i, double * g ) {
1274  double y = 0;
1275  const double * x1 = data.GetPoint(i,y);
1276 
1277  const DataOptions & fitOpt = data.Opt();
1278  bool useBinIntegral = fitOpt.fIntegral && data.HasBinEdges();
1279  bool useBinVolume = (fitOpt.fBinVolume && data.HasBinEdges());
1280 
1281  IntegralEvaluator<> igEval( func, p, useBinIntegral);
1282  double fval = 0;
1283  const double * x2 = 0;
1284  // calculate the bin volume
1285  double binVolume = 1;
1286  std::vector<double> xc;
1287  if (useBinVolume) {
1288  unsigned int ndim = data.NDim();
1289  xc.resize(ndim);
1290  for (unsigned int j = 0; j < ndim; ++j) {
1291  double x2j = data.GetBinUpEdgeComponent(i, j);
1292  binVolume *= std::abs( x2j-x1[j] );
1293  xc[j] = 0.5*(x2j+ x1[j]);
1294  }
1295  // normalize the bin volume using a reference value
1296  binVolume /= data.RefVolume();
1297  }
1298 
1299  const double * x = (useBinVolume) ? &xc.front() : x1;
1300 
1301  if (!useBinIntegral ) {
1302  fval = func ( x, p );
1303  }
1304  else {
1305  // calculate integral normalized (divided by bin volume)
1306  std::vector<double> vx2(data.NDim());
1307  data.GetBinUpEdgeCoordinates(i, vx2.data());
1308  fval = igEval( x1, vx2.data() ) ;
1309  }
1310  if (useBinVolume) fval *= binVolume;
1311 
1312  // logPdf for Poisson: ignore constant term depending on N
1313  fval = std::max(fval, 0.0); // avoid negative or too small values
1314  double logPdf = - fval;
1315  if (y > 0.0) {
1316  // include also constants due to saturate model (see Baker-Cousins paper)
1317  logPdf += y * ROOT::Math::Util::EvalLog( fval / y) + y;
1318  }
1319  // need to return the pdf contribution (not the log)
1320 
1321  //double pdfval = std::exp(logPdf);
1322 
1323  //if (g == 0) return pdfval;
1324  if (g == 0) return logPdf;
1325 
1326  unsigned int npar = func.NPar();
1327  const IGradModelFunction * gfunc = dynamic_cast<const IGradModelFunction *>( &func);
1328 
1329  // gradient calculation
1330  if (gfunc != 0) {
1331  //case function provides gradient
1332  if (!useBinIntegral )
1333  gfunc->ParameterGradient( x , p, g );
1334  else {
1335  // needs to calculate the integral for each partial derivative
1336  CalculateGradientIntegral( *gfunc, x1, x2, p, g);
1337  }
1338 
1339  }
1340  else {
1341  SimpleGradientCalculator gc(func.NPar(), func);
1342  if (!useBinIntegral )
1343  gc.ParameterGradient(x, p, fval, g);
1344  else {
1345  // needs to calculate the integral for each partial derivative
1346  CalculateGradientIntegral( gc, x1, x2, p, g);
1347  }
1348  }
1349  // correct g[] do be derivative of poisson term
1350  for (unsigned int k = 0; k < npar; ++k) {
1351  // apply bin volume correction
1352  if (useBinVolume) g[k] *= binVolume;
1353 
1354  // correct for Poisson term
1355  if ( fval > 0)
1356  g[k] *= ( y/fval - 1.) ;//* pdfval;
1357  else if (y > 0) {
1358  const double kdmax1 = std::sqrt( std::numeric_limits<double>::max() );
1359  g[k] *= kdmax1;
1360  }
1361  else // y == 0 cannot have negative y
1362  g[k] *= -1;
1363  }
1364 
1365 
1366 #ifdef DEBUG
1367  std::cout << "x = " << x[0] << " logPdf = " << logPdf << " grad";
1368  for (unsigned int ipar = 0; ipar < npar; ++ipar)
1369  std::cout << g[ipar] << "\t";
1370  std::cout << std::endl;
1371 #endif
1372 
1373 // return pdfval;
1374  return logPdf;
1375 }
1376 
1377 double FitUtil::EvaluatePoissonLogL(const IModelFunction &func, const BinData &data, const double *p, int iWeight,
1378  bool extended, unsigned int &nPoints, ROOT::Fit::ExecutionPolicy executionPolicy,
1379  unsigned nChunks)
1380 {
1381  // evaluate the Poisson Log Likelihood
1382  // for binned likelihood fits
1383  // this is Sum ( f(x_i) - y_i * log( f (x_i) ) )
1384  // add as well constant term for saturated model to make it like a Chi2/2
1385  // by default is etended. If extended is false the fit is not extended and
1386  // the global poisson term is removed (i.e is a binomial fit)
1387  // (remember that in this case one needs to have a function with a fixed normalization
1388  // like in a non extended unbinned fit)
1389  //
1390  // if use Weight use a weighted dataset
1391  // iWeight = 1 ==> logL = Sum( w f(x_i) )
1392  // case of iWeight==1 is actually identical to weight==0
1393  // iWeight = 2 ==> logL = Sum( w*w * f(x_i) )
1394  //
1395  // nPoints returns the points where bin content is not zero
1396 
1397 
1398  unsigned int n = data.Size();
1399 
1400 #ifdef USE_PARAMCACHE
1401  (const_cast<IModelFunction &>(func)).SetParameters(p);
1402 #endif
1403 
1404  nPoints = data.Size(); // npoints
1405 
1406 
1407  // get fit option and check case of using integral of bins
1408  const DataOptions &fitOpt = data.Opt();
1409  bool useBinIntegral = fitOpt.fIntegral && data.HasBinEdges();
1410  bool useBinVolume = (fitOpt.fBinVolume && data.HasBinEdges());
1411  bool useW2 = (iWeight == 2);
1412 
1413  // normalize if needed by a reference volume value
1414  double wrefVolume = 1.0;
1415  if (useBinVolume) {
1416  if (fitOpt.fNormBinVolume) wrefVolume /= data.RefVolume();
1417  }
1418 
1419 //#define DEBUG
1420 #ifdef DEBUG
1421  std::cout << "Evaluate PoissonLogL for params = [ ";
1422  for (unsigned int j = 0; j < func.NPar(); ++j) std::cout << p[j] << " , ";
1423  std::cout << "] - data size = " << n << " useBinIntegral " << useBinIntegral << " useBinVolume "
1424  << useBinVolume << " useW2 " << useW2 << " wrefVolume = " << wrefVolume << std::endl;
1425 #endif
1426 
1427 
1429  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
1430  // do not use GSL integrator which is not thread safe
1432  }
1433 #ifdef USE_PARAMCACHE
1434  IntegralEvaluator<> igEval(func, 0, useBinIntegral, igType);
1435 #else
1436  IntegralEvaluator<> igEval(func, p, useBinIntegral, igType);
1437 #endif
1438 
1439  auto mapFunction = [&](const unsigned i) {
1440  auto x1 = data.GetCoordComponent(i, 0);
1441  auto y = *data.ValuePtr(i);
1442 
1443  const double *x = nullptr;
1444  std::vector<double> xc;
1445  double fval = 0;
1446  double binVolume = 1.0;
1447 
1448  if (useBinVolume) {
1449  unsigned int ndim = data.NDim();
1450  xc.resize(data.NDim());
1451  for (unsigned int j = 0; j < ndim; ++j) {
1452  double xx = *data.GetCoordComponent(i, j);
1453  double x2 = data.GetBinUpEdgeComponent(i, j);
1454  binVolume *= std::abs(x2 - xx);
1455  xc[j] = 0.5 * (x2 + xx);
1456  }
1457  x = xc.data();
1458  // normalize the bin volume using a reference value
1459  binVolume *= wrefVolume;
1460  } else if (data.NDim() > 1) {
1461  xc.resize(data.NDim());
1462  xc[0] = *x1;
1463  for (unsigned int j = 1; j < data.NDim(); ++j) {
1464  xc[j] = *data.GetCoordComponent(i, j);
1465  }
1466  x = xc.data();
1467  } else {
1468  x = x1;
1469  }
1470 
1471  if (!useBinIntegral) {
1472 #ifdef USE_PARAMCACHE
1473  fval = func(x);
1474 #else
1475  fval = func(x, p);
1476 #endif
1477  } else {
1478  // calculate integral (normalized by bin volume)
1479  // need to set function and parameters here in case loop is parallelized
1480  std::vector<double> x2(data.NDim());
1481  data.GetBinUpEdgeCoordinates(i, x2.data());
1482  fval = igEval(x, x2.data());
1483  }
1484  if (useBinVolume) fval *= binVolume;
1485 
1486 
1487 
1488 #ifdef DEBUG
1489  int NSAMPLE = 100;
1490  if (i % NSAMPLE == 0) {
1491  std::cout << "evt " << i << " x = [ ";
1492  for (unsigned int j = 0; j < func.NDim(); ++j) std::cout << x[j] << " , ";
1493  std::cout << "] ";
1494  if (fitOpt.fIntegral) {
1495  std::cout << "x2 = [ ";
1496  for (unsigned int j = 0; j < func.NDim(); ++j) std::cout << data.GetBinUpEdgeComponent(i, j) << " , ";
1497  std::cout << "] ";
1498  }
1499  std::cout << " y = " << y << " fval = " << fval << std::endl;
1500  }
1501 #endif
1502 
1503 
1504  // EvalLog protects against 0 values of fval but don't want to add in the -log sum
1505  // negative values of fval
1506  fval = std::max(fval, 0.0);
1507 
1508  double nloglike = 0; // negative loglikelihood
1509  if (useW2) {
1510  // apply weight correction . Effective weight is error^2/ y
1511  // and expected events in bins is fval/weight
1512  // can apply correction only when y is not zero otherwise weight is undefined
1513  // (in case of weighted likelihood I don't care about the constant term due to
1514  // the saturated model)
1515 
1516  // use for the empty bins the global weight
1517  double weight = 1.0;
1518  if (y != 0) {
1519  double error = data.Error(i);
1520  weight = (error * error) / y; // this is the bin effective weight
1521  nloglike += weight * y * ( ROOT::Math::Util::EvalLog(y) - ROOT::Math::Util::EvalLog(fval) );
1522  }
1523  else {
1524  // for empty bin use the average weight computed from the total data weight
1525  weight = data.SumOfError2()/ data.SumOfContent();
1526  }
1527  if (extended) {
1528  nloglike += weight * ( fval - y);
1529  }
1530 
1531  } else {
1532  // standard case no weights or iWeight=1
1533  // this is needed for Poisson likelihood (which are extened and not for multinomial)
1534  // the formula below include constant term due to likelihood of saturated model (f(x) = y)
1535  // (same formula as in Baker-Cousins paper, page 439 except a factor of 2
1536  if (extended) nloglike = fval - y;
1537 
1538  if (y > 0) {
1539  nloglike += y * (ROOT::Math::Util::EvalLog(y) - ROOT::Math::Util::EvalLog(fval));
1540  }
1541  }
1542 #ifdef DEBUG
1543  {
1545  std::cout << " nll = " << nloglike << std::endl;
1546  }
1547 #endif
1548  return nloglike;
1549  };
1550 
1551 #ifdef R__USE_IMT
1552  auto redFunction = [](const std::vector<double> &objs) {
1553  return std::accumulate(objs.begin(), objs.end(), double{});
1554  };
1555 #else
1556  (void)nChunks;
1557 
1558  // If IMT is disabled, force the execution policy to the serial case
1559  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
1560  Warning("FitUtil::EvaluatePoissonLogL", "Multithread execution policy requires IMT, which is disabled. Changing "
1561  "to ROOT::Fit::ExecutionPolicy::kSerial.");
1562  executionPolicy = ROOT::Fit::ExecutionPolicy::kSerial;
1563  }
1564 #endif
1565 
1566  double res{};
1567  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kSerial) {
1568  for (unsigned int i = 0; i < n; ++i) {
1569  res += mapFunction(i);
1570  }
1571 #ifdef R__USE_IMT
1572  } else if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
1573  ROOT::TThreadExecutor pool;
1574  auto chunks = nChunks != 0 ? nChunks : setAutomaticChunking(data.Size());
1575  res = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, n), redFunction, chunks);
1576 #endif
1577  // } else if(executionPolicy == ROOT::Fit::kMultitProcess){
1578  // ROOT::TProcessExecutor pool;
1579  // res = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, n), redFunction);
1580  } else {
1581  Error("FitUtil::EvaluatePoissonLogL",
1582  "Execution policy unknown. Avalaible choices:\n ROOT::Fit::ExecutionPolicy::kSerial (default)\n ROOT::Fit::ExecutionPolicy::kMultithread (requires IMT)\n");
1583  }
1584 
1585 #ifdef DEBUG
1586  std::cout << "Loglikelihood = " << res << std::endl;
1587 #endif
1588 
1589  return res;
1590 }
1591 
1592 void FitUtil::EvaluatePoissonLogLGradient(const IModelFunction &f, const BinData &data, const double *p, double *grad,
1593  unsigned int &, ROOT::Fit::ExecutionPolicy executionPolicy, unsigned nChunks)
1594 {
1595  // evaluate the gradient of the Poisson log likelihood function
1596 
1597  const IGradModelFunction *fg = dynamic_cast<const IGradModelFunction *>(&f);
1598  assert(fg != nullptr); // must be called by a grad function
1599 
1600  const IGradModelFunction &func = *fg;
1601 
1602 #ifdef USE_PARAMCACHE
1603  (const_cast<IGradModelFunction &>(func)).SetParameters(p);
1604 #endif
1605 
1606  const DataOptions &fitOpt = data.Opt();
1607  bool useBinIntegral = fitOpt.fIntegral && data.HasBinEdges();
1608  bool useBinVolume = (fitOpt.fBinVolume && data.HasBinEdges());
1609 
1610  double wrefVolume = 1.0;
1611  if (useBinVolume && fitOpt.fNormBinVolume)
1612  wrefVolume /= data.RefVolume();
1613 
1615  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
1616  // do not use GSL integrator which is not thread safe
1618  }
1619 
1620  IntegralEvaluator<> igEval(func, p, useBinIntegral, igType);
1621 
1622  unsigned int npar = func.NPar();
1623  unsigned initialNPoints = data.Size();
1624 
1625  auto mapFunction = [&](const unsigned int i) {
1626  // set all vector values to zero
1627  std::vector<double> gradFunc(npar);
1628  std::vector<double> pointContribution(npar);
1629 
1630  const auto x1 = data.GetCoordComponent(i, 0);
1631  const auto y = data.Value(i);
1632  auto invError = data.Error(i);
1633 
1634  invError = (invError != 0.0) ? 1.0 / invError : 1;
1635 
1636  double fval = 0;
1637 
1638  const double *x = nullptr;
1639  std::vector<double> xc;
1640 
1641  unsigned ndim = data.NDim();
1642  double binVolume = 1.0;
1643  if (useBinVolume) {
1644 
1645  xc.resize(ndim);
1646 
1647  for (unsigned int j = 0; j < ndim; ++j) {
1648  double x1_j = *data.GetCoordComponent(i, j);
1649  double x2_j = data.GetBinUpEdgeComponent(i, j);
1650  binVolume *= std::abs(x2_j - x1_j);
1651  xc[j] = 0.5 * (x2_j + x1_j);
1652  }
1653 
1654  x = xc.data();
1655 
1656  // normalize the bin volume using a reference value
1657  binVolume *= wrefVolume;
1658  } else if (ndim > 1) {
1659  xc.resize(ndim);
1660  xc[0] = *x1;
1661  for (unsigned int j = 1; j < ndim; ++j)
1662  xc[j] = *data.GetCoordComponent(i, j);
1663  x = xc.data();
1664  } else {
1665  x = x1;
1666  }
1667 
1668  if (!useBinIntegral) {
1669  fval = func(x, p);
1670  func.ParameterGradient(x, p, &gradFunc[0]);
1671  } else {
1672  // calculate integral (normalized by bin volume)
1673  // need to set function and parameters here in case loop is parallelized
1674  std::vector<double> x2(data.NDim());
1675  data.GetBinUpEdgeCoordinates(i, x2.data());
1676  fval = igEval(x, x2.data());
1677  CalculateGradientIntegral(func, x, x2.data(), p, &gradFunc[0]);
1678  }
1679  if (useBinVolume)
1680  fval *= binVolume;
1681 
1682 #ifdef DEBUG
1683  {
1685  if (i < 5 || (i > data.Size()-5) ) {
1686  if (data.NDim() > 1) std::cout << i << " x " << x[0] << " y " << x[1] << " func " << fval
1687  << " gradient " << gradFunc[0] << " " << gradFunc[1] << " " << gradFunc[3] << std::endl;
1688  else std::cout << i << " x " << x[0] << " gradient " << gradFunc[0] << " " << gradFunc[1] << " " << gradFunc[3] << std::endl;
1689  }
1690  }
1691 #endif
1692 
1693  // correct the gradient
1694  for (unsigned int ipar = 0; ipar < npar; ++ipar) {
1695 
1696  // correct gradient for bin volumes
1697  if (useBinVolume)
1698  gradFunc[ipar] *= binVolume;
1699 
1700  // df/dp * (1. - y/f )
1701  if (fval > 0)
1702  pointContribution[ipar] = gradFunc[ipar] * (1. - y / fval);
1703  else if (gradFunc[ipar] != 0) {
1704  const double kdmax1 = std::sqrt(std::numeric_limits<double>::max());
1705  const double kdmax2 = std::numeric_limits<double>::max() / (4 * initialNPoints);
1706  double gg = kdmax1 * gradFunc[ipar];
1707  if (gg > 0)
1708  gg = std::min(gg, kdmax2);
1709  else
1710  gg = std::max(gg, -kdmax2);
1711  pointContribution[ipar] = -gg;
1712  }
1713  }
1714 
1715 
1716  return pointContribution;
1717  };
1718 
1719  // Vertically reduce the set of vectors by summing its equally-indexed components
1720  auto redFunction = [&](const std::vector<std::vector<double>> &pointContributions) {
1721  std::vector<double> result(npar);
1722 
1723  for (auto const &pointContribution : pointContributions) {
1724  for (unsigned int parameterIndex = 0; parameterIndex < npar; parameterIndex++)
1725  result[parameterIndex] += pointContribution[parameterIndex];
1726  }
1727 
1728  return result;
1729  };
1730 
1731  std::vector<double> g(npar);
1732 
1733 #ifndef R__USE_IMT
1734  // If IMT is disabled, force the execution policy to the serial case
1735  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
1736  Warning("FitUtil::EvaluatePoissonLogLGradient",
1737  "Multithread execution policy requires IMT, which is disabled. Changing "
1738  "to ROOT::Fit::ExecutionPolicy::kSerial.");
1739  executionPolicy = ROOT::Fit::ExecutionPolicy::kSerial;
1740  }
1741 #endif
1742 
1743  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kSerial) {
1744  std::vector<std::vector<double>> allGradients(initialNPoints);
1745  for (unsigned int i = 0; i < initialNPoints; ++i) {
1746  allGradients[i] = mapFunction(i);
1747  }
1748  g = redFunction(allGradients);
1749  }
1750 #ifdef R__USE_IMT
1751  else if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
1752  ROOT::TThreadExecutor pool;
1753  auto chunks = nChunks != 0 ? nChunks : setAutomaticChunking(initialNPoints);
1754  g = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, initialNPoints), redFunction, chunks);
1755  }
1756 #endif
1757 
1758  // else if(executionPolicy == ROOT::Fit::kMultiprocess){
1759  // ROOT::TProcessExecutor pool;
1760  // g = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, n), redFunction);
1761  // }
1762  else {
1763  Error("FitUtil::EvaluatePoissonLogLGradient",
1764  "Execution policy unknown. Avalaible choices:\n 0: Serial (default)\n 1: MultiThread (requires IMT)\n");
1765  }
1766 
1767 #ifndef R__USE_IMT
1768  //to fix compiler warning
1769  (void)nChunks;
1770 #endif
1771 
1772  // copy result
1773  std::copy(g.begin(), g.end(), grad);
1774 
1775 #ifdef DEBUG
1776  std::cout << "***** Final gradient : ";
1777  for (unsigned int ii = 0; ii< npar; ++ii) std::cout << grad[ii] << " ";
1778  std::cout << "\n";
1779 #endif
1780 
1781 }
1782 
1783 
1784 unsigned FitUtil::setAutomaticChunking(unsigned nEvents){
1785  auto ncpu = ROOT::GetThreadPoolSize();
1786  if (nEvents/ncpu < 1000) return ncpu;
1787  return nEvents/1000;
1788  //return ((nEvents/ncpu + 1) % 1000) *40 ; //arbitrary formula
1789 }
1790 
1791 }
1792 
1793 } // end namespace ROOT
float xmin
Definition: THbookFile.cxx:93
double
Definition: Converters.cxx:921
bool HaveAsymErrors() const
flag to control if data provides asymmetric errors on the value
Definition: BinData.h:143
unsigned int NDim() const
return coordinate data dimension
Definition: UnBinData.h:280
bool HasBinEdges() const
query if the data store the bin edges instead of the center
Definition: BinData.h:534
Returns the available number of logical cores.
Definition: RNumpyDS.hxx:30
double EvaluateChi2(const IModelFunction &func, const BinData &data, const double *x, unsigned int &nPoints, ROOT::Fit::ExecutionPolicy executionPolicy, unsigned nChunks=0)
Chi2 Functions.
Definition: FitUtil.cxx:226
void EvaluateLogLGradient(const IModelFunction &func, const UnBinData &data, const double *x, double *grad, unsigned int &nPoints, ROOT::Fit::ExecutionPolicy executionPolicy=ROOT::Fit::ExecutionPolicy::kSerial, unsigned nChunks=0)
evaluate the LogL gradient given a model function and the data at the point x.
Definition: FitUtil.cxx:1131
double EvaluatePoissonBinPdf(const IModelFunction &func, const BinData &data, const double *x, unsigned int ipoint, double *g=0)
evaluate the pdf contribution to the Poisson LogL given a model function and the BinPoint data...
Definition: FitUtil.cxx:1273
void CalculateGradientIntegral(const GFunc &gfunc, const double *x1, const double *x2, const double *p, double *g)
Definition: FitUtil.cxx:202
ErrorType GetErrorType() const
retrieve the errortype
Definition: BinData.h:551
#define g(i)
Definition: RSha256.hxx:105
const double * ValuePtr(unsigned int ipoint) const
return a pointer to the value for the given fit point
Definition: BinData.h:229
unsigned setAutomaticChunking(unsigned nEvents)
Definition: FitUtil.cxx:1784
Class describing the unbinned data sets (just x coordinates values) of any dimensions.
Definition: UnBinData.h:42
#define f(i)
Definition: RSha256.hxx:104
Interface (abstract class) for parametric gradient multi-dimensional functions providing in addition ...
R__EXTERN TVirtualMutex * gROOTMutex
Definition: TROOT.h:60
STL namespace.
void GetBinUpEdgeCoordinates(unsigned int ipoint, double *x) const
Thread save version of function retrieving the bin up-edge in case of multidimensions.
Definition: BinData.h:521
bool CheckInfNaNValue(double &rval)
Definition: FitUtil.cxx:181
TRObject operator()(const T1 &t1) const
double SumOfError2() const
compute the total sum of the error square (sum of weight square in case of a weighted data set) ...
Definition: BinData.h:566
const double * Coords(unsigned int ipoint) const
return a pointer to the coordinates data for the given fit point
Definition: FitData.h:246
const double * GetPoint(unsigned int ipoint, double &value) const
retrieve at the same time a pointer to the coordinate data and the fit value More efficient than call...
Definition: BinData.h:370
double sqrt(double)
static const double x2[5]
Double_t x[n]
Definition: legend1.C:17
void GetRange(unsigned int irange, unsigned int icoord, double &xmin, double &xmax) const
get the i-th range for given coordinate.
Definition: DataRange.h:103
OneDimMultiFunctionAdapter class to wrap a multidimensional function in one dimensional one...
const double * GetCoordComponent(unsigned int ipoint, unsigned int icoord) const
returns a single coordinate component of a point.
Definition: FitData.h:229
void EvaluatePoissonLogLGradient(const IModelFunction &func, const BinData &data, const double *x, double *grad, unsigned int &nPoints, ROOT::Fit::ExecutionPolicy executionPolicy=ROOT::Fit::ExecutionPolicy::kSerial, unsigned nChunks=0)
evaluate the Poisson LogL given a model function and the data at the point x.
Definition: FitUtil.cxx:1592
double RefVolume() const
retrieve the reference volume used to normalize the data when the option bin volume is set ...
Definition: BinData.h:541
double Derivative1(double x)
Returns the first derivative of the function at point x, computed by Richardson&#39;s extrapolation metho...
double EvaluatePoissonLogL(const IModelFunction &func, const BinData &data, const double *x, int iWeight, bool extended, unsigned int &nPoints, ROOT::Fit::ExecutionPolicy executionPolicy, unsigned nChunks=0)
evaluate the Poisson LogL given a model function and the data at the point x.
Definition: FitUtil.cxx:1377
unsigned int Size() const
return number of fit points
Definition: FitData.h:303
double Error(unsigned int ipoint) const
Definition: BinData.h:251
#define MATH_ERROR_MSG(loc, str)
Definition: Error.h:82
void Error(const char *location, const char *msgfmt,...)
T EvalLog(T x)
safe evaluation of log(x) with a protections against negative or zero argument to the log smooth line...
Definition: Util.h:64
void EvaluateChi2Gradient(const IModelFunction &func, const BinData &data, const double *x, double *grad, unsigned int &nPoints, ROOT::Fit::ExecutionPolicy executionPolicy=ROOT::Fit::ExecutionPolicy::kSerial, unsigned nChunks=0)
evaluate the Chi2 gradient given a model function and the data at the point x.
Definition: FitUtil.cxx:646
This class provides a simple interface to execute the same task multiple times in parallel...
virtual unsigned int NPar() const =0
Return the number of Parameters.
double Weight(unsigned int ipoint) const
return weight
Definition: UnBinData.h:257
DataOptions : simple structure holding the options on how the data are filled.
Definition: DataOptions.h:28
double EvaluateChi2Effective(const IModelFunction &func, const BinData &data, const double *x, unsigned int &nPoints)
evaluate the effective Chi2 given a model function and the data at the point x.
Definition: FitUtil.cxx:418
const double * BinUpEdge(unsigned int ipoint) const
return an array containing the upper edge of the bin for coordinate i In case of empty bin they could...
Definition: BinData.h:508
unsigned int Size(unsigned int icoord=0) const
return range size for coordinate icoord (starts from zero) Size == 0 indicates no range is present [-...
Definition: DataRange.h:70
double EvaluateChi2Residual(const IModelFunction &func, const BinData &data, const double *x, unsigned int ipoint, double *g=0)
evaluate the residual contribution to the Chi2 given a model function and the BinPoint data and if th...
Definition: FitUtil.cxx:539
const DataOptions & Opt() const
access to options
Definition: FitData.h:319
virtual void ParameterGradient(const T *x, const double *p, T *grad) const
Evaluate the all the derivatives (gradient vector) of the function with respect to the parameters at ...
Class describing the binned data sets : vectors of x coordinates, y values and optionally error on y ...
Definition: BinData.h:53
double GetBinUpEdgeComponent(unsigned int ipoint, unsigned int icoord) const
returns a single coordinate error component of a point.
Definition: BinData.h:491
float xmax
Definition: THbookFile.cxx:93
bool HaveCoordErrors() const
flag to control if data provides error on the coordinates
Definition: BinData.h:131
void Warning(const char *location, const char *msgfmt,...)
constexpr Double_t E()
Base of natural log: .
Definition: TMath.h:97
Type
enumeration specifying the integration types.
#define h(i)
Definition: RSha256.hxx:106
ROOT::Math::IParamMultiFunction IModelFunction
Definition: FitUtil.h:63
double CorrectValue(double rval)
Definition: FitUtil.cxx:167
virtual void SetParameters(const double *p)=0
Set the parameter values.
static const double x1[5]
A pseudo container class which is a generator of indices.
Definition: TSeq.hxx:66
TFitResultPtr Fit(FitObject *h1, TF1 *f1, Foption_t &option, const ROOT::Math::MinimizerOptions &moption, const char *goption, ROOT::Fit::DataRange &range)
Definition: HFitImpl.cxx:133
double EvaluatePdf(const IModelFunction &func, const UnBinData &data, const double *x, unsigned int ipoint, double *g=0)
evaluate the pdf contribution to the LogL given a model function and the BinPoint data...
Definition: FitUtil.cxx:863
Double_t y[n]
Definition: legend1.C:17
Double_t ey[n]
Definition: legend1.C:17
#define R__LOCKGUARD(mutex)
double EvaluateLogL(const IModelFunction &func, const UnBinData &data, const double *p, int iWeight, bool extended, unsigned int &nPoints, ROOT::Fit::ExecutionPolicy executionPolicy, unsigned nChunks=0)
evaluate the LogL given a model function and the data at the point x.
double Value(unsigned int ipoint) const
return the value for the given fit point
Definition: BinData.h:217
double SumOfContent() const
compute the total sum of the data content (sum of weights in case of weighted data set) ...
Definition: BinData.h:560
typedef void((*Func_t)())
auto MapReduce(F func, unsigned nTimes, R redfunc) -> typename std::result_of< F()>::type
This method behaves just like Map, but an additional redfunc function must be provided.
bool IsWeighted() const
return true if the data set is weighted We cannot compute ourselfs because sometimes errors are fille...
Definition: BinData.h:573
auto * l
Definition: textangle.C:4
const double * GetPointError(unsigned int ipoint, double &errvalue) const
Retrieve the errors on the point (coordinate and value) for the given fit point It must be called onl...
Definition: BinData.h:451
TF1 * f1
Definition: legend1.C:11
double InvError(unsigned int ipoint) const
Return the inverse of error on the value for the given fit point useful when error in the coordinates...
Definition: BinData.h:314
unsigned int NDim() const
return coordinate data dimension
Definition: FitData.h:311
const DataRange & Range() const
access to range
Definition: FitData.h:331
UInt_t GetThreadPoolSize()
Returns the size of ROOT&#39;s thread pool.
Definition: TROOT.cxx:560
Double_t ex[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
User class for calculating the derivatives of a function.
virtual unsigned int NDim() const =0
Retrieve the dimension of the function.