Logo ROOT   6.19/01
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  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) {
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
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
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 
678  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
679  // do not use GSL integrator which is not thread safe
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) {
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  {
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 
1425  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
1426  // do not use GSL integrator which is not thread safe
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  {
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 
1609  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
1610  // do not use GSL integrator which is not thread safe
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  {
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
UInt_t GetImplicitMTPoolSize()
Returns the size of the pool used for implicit multi-threading.
Definition: TROOT.cxx:618
float xmin
Definition: THbookFile.cxx:93
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:528
VSD Structures.
Definition: StringConv.hxx:21
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:1128
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:1270
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:545
#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:1776
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:59
STL namespace.
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:560
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:1586
double RefVolume() const
retrieve the reference volume used to normalize the data when the option bin volume is set ...
Definition: BinData.h:535
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:1373
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:643
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:415
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:536
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
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:64
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:134
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:860
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:554
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:567
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
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.