Logo ROOT  
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"
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
46namespace 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>
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
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;
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
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>
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 &nPoints,
227 ::ROOT::EExecutionPolicy 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::EExecutionPolicy::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::EExecutionPolicy::kMultiThread) {
388 Warning("FitUtil::EvaluateChi2", "Multithread execution policy requires IMT, which is disabled. Changing "
389 "to ROOT::EExecutionPolicy::kSequential.");
390 executionPolicy = ROOT::EExecutionPolicy::kSequential;
391 }
392#endif
393
394 double res{};
395 if(executionPolicy == ROOT::EExecutionPolicy::kSequential){
396 for (unsigned int i=0; i<n; ++i) {
397 res += mapFunction(i);
398 }
399#ifdef R__USE_IMT
400 } else if(executionPolicy == ROOT::EExecutionPolicy::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::EExecutionPolicy::kSequential (default)\n ROOT::EExecutionPolicy::kMultiThread (requires IMT)\n");
410 }
411
412 // reset the number of fitting data points
413 nPoints = n; // no points are rejected
414 //if (nRejected != 0) nPoints = n - nRejected;
415
416 return res;
417}
418
419
420//___________________________________________________________________________________________________________________________
421
422double FitUtil::EvaluateChi2Effective(const IModelFunction & func, const BinData & data, const double * p, unsigned int & nPoints) {
423 // evaluate the chi2 given a function reference , the data and returns the value and also in nPoints
424 // the actual number of used points
425 // method using the error in the coordinates
426 // integral of bin does not make sense in this case
427
428 unsigned int n = data.Size();
429
430#ifdef DEBUG
431 std::cout << "\n\nFit data size = " << n << std::endl;
432 std::cout << "evaluate effective chi2 using function " << &func << " " << p << std::endl;
433#endif
434
435 assert(data.HaveCoordErrors() || data.HaveAsymErrors());
436
437 double chi2 = 0;
438 //int nRejected = 0;
439
440
441 //func.SetParameters(p);
442
443 unsigned int ndim = func.NDim();
444
445 // use Richardson derivator
447
448 double maxResValue = std::numeric_limits<double>::max() /n;
449
450
451
452 for (unsigned int i = 0; i < n; ++ i) {
453
454
455 double y = 0;
456 const double * x = data.GetPoint(i,y);
457
458 double fval = func( x, p );
459
460 double delta_y_func = y - fval;
461
462
463 double ey = 0;
464 const double * ex = 0;
465 if (!data.HaveAsymErrors() )
466 ex = data.GetPointError(i, ey);
467 else {
468 double eylow, eyhigh = 0;
469 ex = data.GetPointError(i, eylow, eyhigh);
470 if ( delta_y_func < 0)
471 ey = eyhigh; // function is higher than points
472 else
473 ey = eylow;
474 }
475 double e2 = ey * ey;
476 // before calculating the gradient check that all error in x are not zero
477 unsigned int j = 0;
478 while ( j < ndim && ex[j] == 0.) { j++; }
479 // if j is less ndim some elements are not zero
480 if (j < ndim) {
481 // need an adapter from a multi-dim function to a one-dimensional
483 // select optimal step size (use 10--2 by default as was done in TF1:
484 double kEps = 0.01;
485 double kPrecision = 1.E-8;
486 for (unsigned int icoord = 0; icoord < ndim; ++icoord) {
487 // calculate derivative for each coordinate
488 if (ex[icoord] > 0) {
489 //gradCalc.Gradient(x, p, fval, &grad[0]);
490 f1D.SetCoord(icoord);
491 // optimal spep size (take ex[] as scale for the points and 1% of it
492 double x0= x[icoord];
493 double h = std::max( kEps* std::abs(ex[icoord]), 8.0*kPrecision*(std::abs(x0) + kPrecision) );
494 double deriv = derivator.Derivative1(f1D, x[icoord], h);
495 double edx = ex[icoord] * deriv;
496 e2 += edx * edx;
497#ifdef DEBUG
498 std::cout << "error for coord " << icoord << " = " << ex[icoord] << " deriv " << deriv << std::endl;
499#endif
500 }
501 }
502 }
503 double w2 = (e2 > 0) ? 1.0/e2 : 0;
504 double resval = w2 * ( y - fval ) * ( y - fval);
505
506#ifdef DEBUG
507 std::cout << x[0] << " " << y << " ex " << ex[0] << " ey " << ey << " params : ";
508 for (unsigned int ipar = 0; ipar < func.NPar(); ++ipar)
509 std::cout << p[ipar] << "\t";
510 std::cout << "\tfval = " << fval << "\tresval = " << resval << std::endl;
511#endif
512
513 // avoid (infinity and nan ) in the chi2 sum
514 // eventually add possibility of excluding some points (like singularity)
515 if ( resval < maxResValue )
516 chi2 += resval;
517 else
518 chi2 += maxResValue;
519 //nRejected++;
520
521 }
522
523 // reset the number of fitting data points
524 nPoints = n; // no points are rejected
525 //if (nRejected != 0) nPoints = n - nRejected;
526
527#ifdef DEBUG
528 std::cout << "chi2 = " << chi2 << " n = " << nPoints << std::endl;
529#endif
530
531 return chi2;
532
533}
534
535
536////////////////////////////////////////////////////////////////////////////////
537/// evaluate the chi2 contribution (residual term) only for data with no coord-errors
538/// This function is used in the specialized least square algorithms like FUMILI or L.M.
539/// if we have error on the coordinates the method is not yet implemented
540/// integral option is also not yet implemented
541/// one can use in that case normal chi2 method
542
543double FitUtil::EvaluateChi2Residual(const IModelFunction & func, const BinData & data, const double * p, unsigned int i, double * g) {
544 if (data.GetErrorType() == BinData::kCoordError && data.Opt().fCoordErrors ) {
545 MATH_ERROR_MSG("FitUtil::EvaluateChi2Residual","Error on the coordinates are not used in calculating Chi2 residual");
546 return 0; // it will assert otherwise later in GetPoint
547 }
548
549
550 //func.SetParameters(p);
551
552 double y, invError = 0;
553
554 const DataOptions & fitOpt = data.Opt();
555 bool useBinIntegral = fitOpt.fIntegral && data.HasBinEdges();
556 bool useBinVolume = (fitOpt.fBinVolume && data.HasBinEdges());
557 bool useExpErrors = (fitOpt.fExpErrors);
558
559 const double * x1 = data.GetPoint(i,y, invError);
560
561 IntegralEvaluator<> igEval( func, p, useBinIntegral);
562 double fval = 0;
563 unsigned int ndim = data.NDim();
564 double binVolume = 1.0;
565 const double * x2 = 0;
566 if (useBinVolume || useBinIntegral) x2 = data.BinUpEdge(i);
567
568 double * xc = 0;
569
570 if (useBinVolume) {
571 xc = new double[ndim];
572 for (unsigned int j = 0; j < ndim; ++j) {
573 binVolume *= std::abs( x2[j]-x1[j] );
574 xc[j] = 0.5*(x2[j]+ x1[j]);
575 }
576 // normalize the bin volume using a reference value
577 binVolume /= data.RefVolume();
578 }
579
580 const double * x = (useBinVolume) ? xc : x1;
581
582 if (!useBinIntegral) {
583 fval = func ( x, p );
584 }
585 else {
586 // calculate integral (normalized by bin volume)
587 // need to set function and parameters here in case loop is parallelized
588 fval = igEval( x1, x2) ;
589 }
590 // normalize result if requested according to bin volume
591 if (useBinVolume) fval *= binVolume;
592
593 // expected errors
594 if (useExpErrors) {
595 // we need first to check if a weight factor needs to be applied
596 // weight = sumw2/sumw = error**2/content
597 //NOTE: assume histogram is not weighted
598 // don't know how to do with bins with weight = 0
599 //double invWeight = y * invError * invError;
600 // if (invError == 0) invWeight = (data.SumOfError2() > 0) ? data.SumOfContent()/ data.SumOfError2() : 1.0;
601 // compute expected error as f(x) / weight
602 double invError2 = (fval > 0) ? 1.0 / fval : 0.0;
603 invError = std::sqrt(invError2);
604 }
605
606
607 double resval = ( y -fval )* invError;
608
609 // avoid infinities or nan in resval
610 resval = CorrectValue(resval);
611
612 // estimate gradient
613 if (g != 0) {
614
615 unsigned int npar = func.NPar();
616 const IGradModelFunction * gfunc = dynamic_cast<const IGradModelFunction *>( &func);
617
618 if (gfunc != 0) {
619 //case function provides gradient
620 if (!useBinIntegral ) {
621 gfunc->ParameterGradient( x , p, g );
622 }
623 else {
624 // needs to calculate the integral for each partial derivative
625 CalculateGradientIntegral( *gfunc, x1, x2, p, g);
626 }
627 }
628 else {
629 SimpleGradientCalculator gc( npar, func);
630 if (!useBinIntegral )
631 gc.ParameterGradient(x, p, fval, g);
632 else {
633 // needs to calculate the integral for each partial derivative
634 CalculateGradientIntegral( gc, x1, x2, p, g);
635 }
636 }
637 // mutiply by - 1 * weight
638 for (unsigned int k = 0; k < npar; ++k) {
639 g[k] *= - invError;
640 if (useBinVolume) g[k] *= binVolume;
641 }
642 }
643
644 if (useBinVolume) delete [] xc;
645
646 return resval;
647
648}
649
650void FitUtil::EvaluateChi2Gradient(const IModelFunction &f, const BinData &data, const double *p, double *grad,
651 unsigned int &nPoints, ROOT::EExecutionPolicy executionPolicy, unsigned nChunks)
652{
653 // evaluate the gradient of the chi2 function
654 // this function is used when the model function knows how to calculate the derivative and we can
655 // avoid that the minimizer re-computes them
656 //
657 // case of chi2 effective (errors on coordinate) is not supported
658
659 if (data.HaveCoordErrors()) {
660 MATH_ERROR_MSG("FitUtil::EvaluateChi2Gradient",
661 "Error on the coordinates are not used in calculating Chi2 gradient");
662 return; // it will assert otherwise later in GetPoint
663 }
664
665 const IGradModelFunction *fg = dynamic_cast<const IGradModelFunction *>(&f);
666 assert(fg != nullptr); // must be called by a gradient function
667
668 const IGradModelFunction &func = *fg;
669
670#ifdef DEBUG
671 std::cout << "\n\nFit data size = " << n << std::endl;
672 std::cout << "evaluate chi2 using function gradient " << &func << " " << p << std::endl;
673#endif
674
675 const DataOptions &fitOpt = data.Opt();
676 bool useBinIntegral = fitOpt.fIntegral && data.HasBinEdges();
677 bool useBinVolume = (fitOpt.fBinVolume && data.HasBinEdges());
678
679 double wrefVolume = 1.0;
680 if (useBinVolume) {
681 if (fitOpt.fNormBinVolume) wrefVolume /= data.RefVolume();
682 }
683
685 if (executionPolicy == ROOT::EExecutionPolicy::kMultiThread) {
686 // do not use GSL integrator which is not thread safe
688 }
689 IntegralEvaluator<> igEval(func, p, useBinIntegral,igType);
690
691 unsigned int npar = func.NPar();
692 unsigned initialNPoints = data.Size();
693
694 std::vector<bool> isPointRejected(initialNPoints);
695
696 auto mapFunction = [&](const unsigned int i) {
697 // set all vector values to zero
698 std::vector<double> gradFunc(npar);
699 std::vector<double> pointContribution(npar);
700
701 const auto x1 = data.GetCoordComponent(i, 0);
702 const auto y = data.Value(i);
703 auto invError = data.Error(i);
704
705 invError = (invError != 0.0) ? 1.0 / invError : 1;
706
707 double fval = 0;
708
709 const double *x = nullptr;
710 std::vector<double> xc;
711
712 unsigned int ndim = data.NDim();
713 double binVolume = 1;
714 if (useBinVolume) {
715 xc.resize(ndim);
716 for (unsigned int j = 0; j < ndim; ++j) {
717 double x1_j = *data.GetCoordComponent(i, j);
718 double x2_j = data.GetBinUpEdgeComponent(i, j);
719 binVolume *= std::abs(x2_j - x1_j);
720 xc[j] = 0.5 * (x2_j + x1_j);
721 }
722
723 x = xc.data();
724
725 // normalize the bin volume using a reference value
726 binVolume *= wrefVolume;
727 } else if (ndim > 1) {
728 xc.resize(ndim);
729 xc[0] = *x1;
730 for (unsigned int j = 1; j < ndim; ++j)
731 xc[j] = *data.GetCoordComponent(i, j);
732 x = xc.data();
733 } else {
734 x = x1;
735 }
736
737 if (!useBinIntegral) {
738 fval = func(x, p);
739 func.ParameterGradient(x, p, &gradFunc[0]);
740 } else {
741 std::vector<double> x2(data.NDim());
742 data.GetBinUpEdgeCoordinates(i, x2.data());
743 // calculate normalized integral and gradient (divided by bin volume)
744 // need to set function and parameters here in case loop is parallelized
745 fval = igEval(x, x2.data());
746 CalculateGradientIntegral(func, x, x2.data(), p, &gradFunc[0]);
747 }
748 if (useBinVolume)
749 fval *= binVolume;
750
751#ifdef DEBUG
752 std::cout << x[0] << " " << y << " " << 1. / invError << " params : ";
753 for (unsigned int ipar = 0; ipar < npar; ++ipar)
754 std::cout << p[ipar] << "\t";
755 std::cout << "\tfval = " << fval << std::endl;
756#endif
757 if (!CheckInfNaNValue(fval)) {
758 isPointRejected[i] = true;
759 // Return a zero contribution to all partial derivatives on behalf of the current point
760 return pointContribution;
761 }
762
763 // loop on the parameters
764 unsigned int ipar = 0;
765 for (; ipar < npar; ++ipar) {
766
767 // correct gradient for bin volumes
768 if (useBinVolume)
769 gradFunc[ipar] *= binVolume;
770
771 // avoid singularity in the function (infinity and nan ) in the chi2 sum
772 // eventually add possibility of excluding some points (like singularity)
773 double dfval = gradFunc[ipar];
774 if (!CheckInfNaNValue(dfval)) {
775 break; // exit loop on parameters
776 }
777
778 // calculate derivative point contribution
779 pointContribution[ipar] = -2.0 * (y - fval) * invError * invError * gradFunc[ipar];
780 }
781
782 if (ipar < npar) {
783 // case loop was broken for an overflow in the gradient calculation
784 isPointRejected[i] = true;
785 }
786
787 return pointContribution;
788 };
789
790 // Vertically reduce the set of vectors by summing its equally-indexed components
791 auto redFunction = [&](const std::vector<std::vector<double>> &pointContributions) {
792 std::vector<double> result(npar);
793
794 for (auto const &pointContribution : pointContributions) {
795 for (unsigned int parameterIndex = 0; parameterIndex < npar; parameterIndex++)
796 result[parameterIndex] += pointContribution[parameterIndex];
797 }
798
799 return result;
800 };
801
802 std::vector<double> g(npar);
803
804#ifndef R__USE_IMT
805 // If IMT is disabled, force the execution policy to the serial case
806 if (executionPolicy == ROOT::EExecutionPolicy::kMultiThread) {
807 Warning("FitUtil::EvaluateChi2Gradient", "Multithread execution policy requires IMT, which is disabled. Changing "
808 "to ROOT::EExecutionPolicy::kSequential.");
809 executionPolicy = ROOT::EExecutionPolicy::kSequential;
810 }
811#endif
812
813 if (executionPolicy == ROOT::EExecutionPolicy::kSequential) {
814 std::vector<std::vector<double>> allGradients(initialNPoints);
815 for (unsigned int i = 0; i < initialNPoints; ++i) {
816 allGradients[i] = mapFunction(i);
817 }
818 g = redFunction(allGradients);
819 }
820#ifdef R__USE_IMT
821 else if (executionPolicy == ROOT::EExecutionPolicy::kMultiThread) {
823 auto chunks = nChunks != 0 ? nChunks : setAutomaticChunking(initialNPoints);
824 g = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, initialNPoints), redFunction, chunks);
825 }
826#endif
827 // else if(executionPolicy == ROOT::Fit::kMultiprocess){
828 // ROOT::TProcessExecutor pool;
829 // g = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, n), redFunction);
830 // }
831 else {
832 Error("FitUtil::EvaluateChi2Gradient",
833 "Execution policy unknown. Avalaible choices:\n 0: Serial (default)\n 1: MultiThread (requires IMT)\n");
834 }
835
836#ifndef R__USE_IMT
837 //to fix compiler warning
838 (void)nChunks;
839#endif
840
841 // correct the number of points
842 nPoints = initialNPoints;
843
844 if (std::any_of(isPointRejected.begin(), isPointRejected.end(), [](bool point) { return point; })) {
845 unsigned nRejected = std::accumulate(isPointRejected.begin(), isPointRejected.end(), 0);
846 assert(nRejected <= initialNPoints);
847 nPoints = initialNPoints - nRejected;
848
849 if (nPoints < npar)
850 MATH_ERROR_MSG("FitUtil::EvaluateChi2Gradient",
851 "Error - too many points rejected for overflow in gradient calculation");
852 }
853
854 // copy result
855 std::copy(g.begin(), g.end(), grad);
856}
857
858//______________________________________________________________________________________________________
859//
860// Log Likelihood functions
861//_______________________________________________________________________________________________________
862
863// utility function used by the likelihoods
864
865// for LogLikelihood functions
866
867double FitUtil::EvaluatePdf(const IModelFunction & func, const UnBinData & data, const double * p, unsigned int i, double * g) {
868 // evaluate the pdf contribution to the generic logl function in case of bin data
869 // return actually the log of the pdf and its derivatives
870
871
872 //func.SetParameters(p);
873
874
875 const double * x = data.Coords(i);
876 double fval = func ( x, p );
877 double logPdf = ROOT::Math::Util::EvalLog(fval);
878 //return
879 if (g == 0) return logPdf;
880
881 const IGradModelFunction * gfunc = dynamic_cast<const IGradModelFunction *>( &func);
882
883 // gradient calculation
884 if (gfunc != 0) {
885 //case function provides gradient
886 gfunc->ParameterGradient( x , p, g );
887 }
888 else {
889 // estimate gradieant numerically with simple 2 point rule
890 // should probably calculate gradient of log(pdf) is more stable numerically
891 SimpleGradientCalculator gc(func.NPar(), func);
892 gc.ParameterGradient(x, p, fval, g );
893 }
894 // divide gradient by function value since returning the logs
895 for (unsigned int ipar = 0; ipar < func.NPar(); ++ipar) {
896 g[ipar] /= fval; // this should be checked against infinities
897 }
898
899#ifdef DEBUG
900 std::cout << x[i] << "\t";
901 std::cout << "\tpar = [ " << func.NPar() << " ] = ";
902 for (unsigned int ipar = 0; ipar < func.NPar(); ++ipar)
903 std::cout << p[ipar] << "\t";
904 std::cout << "\tfval = " << fval;
905 std::cout << "\tgrad = [ ";
906 for (unsigned int ipar = 0; ipar < func.NPar(); ++ipar)
907 std::cout << g[ipar] << "\t";
908 std::cout << " ] " << std::endl;
909#endif
910
911
912 return logPdf;
913}
914
915double FitUtil::EvaluateLogL(const IModelFunction &func, const UnBinData &data, const double *p,
916 int iWeight, bool extended, unsigned int &nPoints,
917 ROOT::EExecutionPolicy executionPolicy, unsigned nChunks)
918{
919 // evaluate the LogLikelihood
920
921 unsigned int n = data.Size();
922
923 //unsigned int nRejected = 0;
924
925 bool normalizeFunc = false;
926
927 // set parameters of the function to cache integral value
928#ifdef USE_PARAMCACHE
929 (const_cast<IModelFunctionTempl<double> &>(func)).SetParameters(p);
930#endif
931
932 nPoints = data.Size(); // npoints
933
934#ifdef R__USE_IMT
935 // in case parameter needs to be propagated to user function use trick to set parameters by calling one time the function
936 // this will be done in sequential mode and parameters can be set in a thread safe manner
937 if (!normalizeFunc) {
938 if (data.NDim() == 1) {
939 const double * x = data.GetCoordComponent(0,0);
940 func( x, p);
941 }
942 else {
943 std::vector<double> x(data.NDim());
944 for (unsigned int j = 0; j < data.NDim(); ++j)
945 x[j] = *data.GetCoordComponent(0, j);
946 func( x.data(), p);
947 }
948 }
949#endif
950
951 double norm = 1.0;
952 if (normalizeFunc) {
953 // compute integral of the function
954 std::vector<double> xmin(data.NDim());
955 std::vector<double> xmax(data.NDim());
956 IntegralEvaluator<> igEval(func, p, true);
957 // compute integral in the ranges where is defined
958 if (data.Range().Size() > 0) {
959 norm = 0;
960 for (unsigned int ir = 0; ir < data.Range().Size(); ++ir) {
961 data.Range().GetRange(&xmin[0], &xmax[0], ir);
962 norm += igEval.Integral(xmin.data(), xmax.data());
963 }
964 } else {
965 // use (-inf +inf)
966 data.Range().GetRange(&xmin[0], &xmax[0]);
967 // check if funcition is zero at +- inf
968 if (func(xmin.data(), p) != 0 || func(xmax.data(), p) != 0) {
969 MATH_ERROR_MSG("FitUtil::EvaluateLogLikelihood",
970 "A range has not been set and the function is not zero at +/- inf");
971 return 0;
972 }
973 norm = igEval.Integral(&xmin[0], &xmax[0]);
974 }
975 }
976
977 // needed to compue effective global weight in case of extended likelihood
978
979 auto mapFunction = [&](const unsigned i) {
980 double W = 0;
981 double W2 = 0;
982 double fval = 0;
983
984 if (data.NDim() > 1) {
985 std::vector<double> x(data.NDim());
986 for (unsigned int j = 0; j < data.NDim(); ++j)
987 x[j] = *data.GetCoordComponent(i, j);
988#ifdef USE_PARAMCACHE
989 fval = func(x.data());
990#else
991 fval = func(x.data(), p);
992#endif
993
994 // one -dim case
995 } else {
996 const auto x = data.GetCoordComponent(i, 0);
997#ifdef USE_PARAMCACHE
998 fval = func(x);
999#else
1000 fval = func(x, p);
1001#endif
1002 }
1003
1004 if (normalizeFunc)
1005 fval = fval * (1 / norm);
1006
1007 // function EvalLog protects against negative or too small values of fval
1008 double logval = ROOT::Math::Util::EvalLog(fval);
1009 if (iWeight > 0) {
1010 double weight = data.Weight(i);
1011 logval *= weight;
1012 if (iWeight == 2) {
1013 logval *= weight; // use square of weights in likelihood
1014 if (!extended) {
1015 // needed sum of weights and sum of weight square if likelkihood is extended
1016 W = weight;
1017 W2 = weight * weight;
1018 }
1019 }
1020 }
1021 return LikelihoodAux<double>(logval, W, W2);
1022 };
1023
1024#ifdef R__USE_IMT
1025 // auto redFunction = [](const std::vector<LikelihoodAux<double>> & objs){
1026 // return std::accumulate(objs.begin(), objs.end(), LikelihoodAux<double>(0.0,0.0,0.0),
1027 // [](const LikelihoodAux<double> &l1, const LikelihoodAux<double> &l2){
1028 // return l1+l2;
1029 // });
1030 // };
1031 // do not use std::accumulate to be sure to maintain always the same order
1032 auto redFunction = [](const std::vector<LikelihoodAux<double>> & objs){
1033 auto l0 = LikelihoodAux<double>(0.0,0.0,0.0);
1034 for ( auto & l : objs ) {
1035 l0 = l0 + l;
1036 }
1037 return l0;
1038 };
1039#else
1040 (void)nChunks;
1041
1042 // If IMT is disabled, force the execution policy to the serial case
1043 if (executionPolicy == ROOT::EExecutionPolicy::kMultiThread) {
1044 Warning("FitUtil::EvaluateLogL", "Multithread execution policy requires IMT, which is disabled. Changing "
1045 "to ROOT::EExecutionPolicy::kSequential.");
1046 executionPolicy = ROOT::EExecutionPolicy::kSequential;
1047 }
1048#endif
1049
1050 double logl{};
1051 double sumW{};
1052 double sumW2{};
1053 if(executionPolicy == ROOT::EExecutionPolicy::kSequential){
1054 for (unsigned int i=0; i<n; ++i) {
1055 auto resArray = mapFunction(i);
1056 logl+=resArray.logvalue;
1057 sumW+=resArray.weight;
1058 sumW2+=resArray.weight2;
1059 }
1060#ifdef R__USE_IMT
1061 } else if(executionPolicy == ROOT::EExecutionPolicy::kMultiThread) {
1063 auto chunks = nChunks !=0? nChunks: setAutomaticChunking(data.Size());
1064 auto resArray = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, n), redFunction, chunks);
1065 logl=resArray.logvalue;
1066 sumW=resArray.weight;
1067 sumW2=resArray.weight2;
1068#endif
1069// } else if(executionPolicy == ROOT::Fit::kMultitProcess){
1070 // ROOT::TProcessExecutor pool;
1071 // res = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, n), redFunction);
1072 } else{
1073 Error("FitUtil::EvaluateLogL","Execution policy unknown. Avalaible choices:\n ROOT::EExecutionPolicy::kSequential (default)\n ROOT::EExecutionPolicy::kMultiThread (requires IMT)\n");
1074 }
1075
1076 if (extended) {
1077 // add Poisson extended term
1078 double extendedTerm = 0; // extended term in likelihood
1079 double nuTot = 0;
1080 // nuTot is integral of function in the range
1081 // if function has been normalized integral has been already computed
1082 if (!normalizeFunc) {
1083 IntegralEvaluator<> igEval( func, p, true);
1084 std::vector<double> xmin(data.NDim());
1085 std::vector<double> xmax(data.NDim());
1086
1087 // compute integral in the ranges where is defined
1088 if (data.Range().Size() > 0 ) {
1089 nuTot = 0;
1090 for (unsigned int ir = 0; ir < data.Range().Size(); ++ir) {
1091 data.Range().GetRange(&xmin[0],&xmax[0],ir);
1092 nuTot += igEval.Integral(xmin.data(),xmax.data());
1093 }
1094 } else {
1095 // use (-inf +inf)
1096 data.Range().GetRange(&xmin[0],&xmax[0]);
1097 // check if funcition is zero at +- inf
1098 if (func(xmin.data(), p) != 0 || func(xmax.data(), p) != 0) {
1099 MATH_ERROR_MSG("FitUtil::EvaluateLogLikelihood","A range has not been set and the function is not zero at +/- inf");
1100 return 0;
1101 }
1102 nuTot = igEval.Integral(&xmin[0],&xmax[0]);
1103 }
1104
1105 // force to be last parameter value
1106 //nutot = p[func.NDim()-1];
1107 if (iWeight != 2)
1108 extendedTerm = - nuTot; // no need to add in this case n log(nu) since is already computed before
1109 else {
1110 // case use weight square in likelihood : compute total effective weight = sw2/sw
1111 // ignore for the moment case when sumW is zero
1112 extendedTerm = - (sumW2 / sumW) * nuTot;
1113 }
1114
1115 }
1116 else {
1117 nuTot = norm;
1118 extendedTerm = - nuTot + double(n) * ROOT::Math::Util::EvalLog( nuTot);
1119 // in case of weights need to use here sum of weights (to be done)
1120 }
1121 logl += extendedTerm;
1122
1123 }
1124
1125#ifdef DEBUG
1126 std::cout << "Evaluated log L for parameters (";
1127 for (unsigned int ip = 0; ip < func.NPar(); ++ip)
1128 std::cout << " " << p[ip];
1129 std::cout << ") fval = " << -logl << std::endl;
1130#endif
1131
1132 return -logl;
1133}
1134
1135void FitUtil::EvaluateLogLGradient(const IModelFunction &f, const UnBinData &data, const double *p, double *grad,
1136 unsigned int &nPoints, ROOT::EExecutionPolicy executionPolicy, unsigned nChunks)
1137{
1138 // evaluate the gradient of the log likelihood function
1139
1140 const IGradModelFunction *fg = dynamic_cast<const IGradModelFunction *>(&f);
1141 assert(fg != nullptr); // must be called by a grad function
1142
1143 const IGradModelFunction &func = *fg;
1144
1145 unsigned int npar = func.NPar();
1146 unsigned initialNPoints = data.Size();
1147
1148 (const_cast<IGradModelFunction &>(func)).SetParameters(p);
1149
1150#ifdef DEBUG
1151 std::cout << "\n===> Evaluate Gradient for parameters ";
1152 for (unsigned int ip = 0; ip < npar; ++ip)
1153 std::cout << " " << p[ip];
1154 std::cout << "\n";
1155#endif
1156
1157 const double kdmax1 = std::sqrt(std::numeric_limits<double>::max());
1158 const double kdmax2 = std::numeric_limits<double>::max() / (4 * initialNPoints);
1159
1160 auto mapFunction = [&](const unsigned int i) {
1161 std::vector<double> gradFunc(npar);
1162 std::vector<double> pointContribution(npar);
1163
1164
1165 const double * x = nullptr;
1166 std::vector<double> xc;
1167 if (data.NDim() > 1) {
1168 xc.resize(data.NDim() );
1169 for (unsigned int j = 0; j < data.NDim(); ++j)
1170 xc[j] = *data.GetCoordComponent(i, j);
1171 x = xc.data();
1172 } else {
1173 x = data.GetCoordComponent(i, 0);
1174 }
1175
1176 double fval = func(x, p);
1177 func.ParameterGradient(x, p, &gradFunc[0]);
1178
1179#ifdef DEBUG
1180 {
1182 if (i < 5 || (i > data.Size()-5) ) {
1183 if (data.NDim() > 1) std::cout << i << " x " << x[0] << " y " << x[1] << " func " << fval
1184 << " gradient " << gradFunc[0] << " " << gradFunc[1] << " " << gradFunc[3] << std::endl;
1185 else std::cout << i << " x " << x[0] << " gradient " << gradFunc[0] << " " << gradFunc[1] << " " << gradFunc[3] << std::endl;
1186 }
1187 }
1188#endif
1189
1190 for (unsigned int kpar = 0; kpar < npar; ++kpar) {
1191 if (fval > 0)
1192 pointContribution[kpar] = -1. / fval * gradFunc[kpar];
1193 else if (gradFunc[kpar] != 0) {
1194 double gg = kdmax1 * gradFunc[kpar];
1195 if (gg > 0)
1196 gg = std::min(gg, kdmax2);
1197 else
1198 gg = std::max(gg, -kdmax2);
1199 pointContribution[kpar] = -gg;
1200 }
1201 // if func derivative is zero term is also zero so do not add in g[kpar]
1202 }
1203
1204 return pointContribution;
1205 };
1206
1207 // Vertically reduce the set of vectors by summing its equally-indexed components
1208 auto redFunction = [&](const std::vector<std::vector<double>> &pointContributions) {
1209 std::vector<double> result(npar);
1210
1211 for (auto const &pointContribution : pointContributions) {
1212 for (unsigned int parameterIndex = 0; parameterIndex < npar; parameterIndex++)
1213 result[parameterIndex] += pointContribution[parameterIndex];
1214 }
1215
1216 return result;
1217 };
1218
1219 std::vector<double> g(npar);
1220
1221#ifndef R__USE_IMT
1222 // If IMT is disabled, force the execution policy to the serial case
1223 if (executionPolicy == ROOT::EExecutionPolicy::kMultiThread) {
1224 Warning("FitUtil::EvaluateLogLGradient", "Multithread execution policy requires IMT, which is disabled. Changing "
1225 "to ROOT::EExecutionPolicy::kSequential.");
1226 executionPolicy = ROOT::EExecutionPolicy::kSequential;
1227 }
1228#endif
1229
1230 if (executionPolicy == ROOT::EExecutionPolicy::kSequential) {
1231 std::vector<std::vector<double>> allGradients(initialNPoints);
1232 for (unsigned int i = 0; i < initialNPoints; ++i) {
1233 allGradients[i] = mapFunction(i);
1234 }
1235 g = redFunction(allGradients);
1236 }
1237#ifdef R__USE_IMT
1238 else if (executionPolicy == ROOT::EExecutionPolicy::kMultiThread) {
1240 auto chunks = nChunks != 0 ? nChunks : setAutomaticChunking(initialNPoints);
1241 g = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, initialNPoints), redFunction, chunks);
1242 }
1243#endif
1244 else {
1245 Error("FitUtil::EvaluateLogLGradient", "Execution policy unknown. Avalaible choices:\n "
1246 "ROOT::EExecutionPolicy::kSequential (default)\n "
1247 "ROOT::EExecutionPolicy::kMultiThread (requires IMT)\n");
1248 }
1249
1250#ifndef R__USE_IMT
1251 // to fix compiler warning
1252 (void)nChunks;
1253#endif
1254
1255 // copy result
1256 std::copy(g.begin(), g.end(), grad);
1257 nPoints = data.Size(); // npoints
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
1273double 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
1377double FitUtil::EvaluatePoissonLogL(const IModelFunction &func, const BinData &data, const double *p, int iWeight,
1378 bool extended, unsigned int &nPoints, ROOT::EExecutionPolicy 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::EExecutionPolicy::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::EExecutionPolicy::kMultiThread) {
1560 Warning("FitUtil::EvaluatePoissonLogL", "Multithread execution policy requires IMT, which is disabled. Changing "
1561 "to ROOT::EExecutionPolicy::kSequential.");
1562 executionPolicy = ROOT::EExecutionPolicy::kSequential;
1563 }
1564#endif
1565
1566 double res{};
1567 if (executionPolicy == ROOT::EExecutionPolicy::kSequential) {
1568 for (unsigned int i = 0; i < n; ++i) {
1569 res += mapFunction(i);
1570 }
1571#ifdef R__USE_IMT
1572 } else if (executionPolicy == ROOT::EExecutionPolicy::kMultiThread) {
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::EExecutionPolicy::kSequential (default)\n ROOT::EExecutionPolicy::kMultiThread (requires IMT)\n");
1583 }
1584
1585#ifdef DEBUG
1586 std::cout << "Loglikelihood = " << res << std::endl;
1587#endif
1588
1589 return res;
1590}
1591
1592void FitUtil::EvaluatePoissonLogLGradient(const IModelFunction &f, const BinData &data, const double *p, double *grad,
1593 unsigned int &, ROOT::EExecutionPolicy 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::EExecutionPolicy::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::EExecutionPolicy::kMultiThread) {
1736 Warning("FitUtil::EvaluatePoissonLogLGradient",
1737 "Multithread execution policy requires IMT, which is disabled. Changing "
1738 "to ROOT::EExecutionPolicy::kSequential.");
1739 executionPolicy = ROOT::EExecutionPolicy::kSequential;
1740 }
1741#endif
1742
1743 if (executionPolicy == ROOT::EExecutionPolicy::kSequential) {
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::EExecutionPolicy::kMultiThread) {
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
1784unsigned 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
double
Definition: Converters.cxx:939
#define MATH_ERROR_MSG(loc, str)
Definition: Error.h:83
typedef void(GLAPIENTRYP _GLUfuncptr)(void)
#define f(i)
Definition: RSha256.hxx:104
#define g(i)
Definition: RSha256.hxx:105
#define h(i)
Definition: RSha256.hxx:106
static const double x2[5]
static const double x1[5]
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition: TError.cxx:187
void Warning(const char *location, const char *msgfmt,...)
Use this function in warning situations.
Definition: TError.cxx:231
float xmin
Definition: THbookFile.cxx:95
float xmax
Definition: THbookFile.cxx:95
R__EXTERN TVirtualMutex * gROOTMutex
Definition: TROOT.h:61
#define R__LOCKGUARD(mutex)
Definition: TF1.cxx:153
Class describing the binned data sets : vectors of x coordinates, y values and optionally error on y ...
Definition: BinData.h:52
double SumOfContent() const
compute the total sum of the data content (sum of weights in case of weighted data set)
Definition: BinData.h:559
double RefVolume() const
retrieve the reference volume used to normalize the data when the option bin volume is set
Definition: BinData.h:540
double GetBinUpEdgeComponent(unsigned int ipoint, unsigned int icoord) const
returns a single coordinate error component of a point.
Definition: BinData.h:490
bool HasBinEdges() const
query if the data store the bin edges instead of the center
Definition: BinData.h:533
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:450
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:369
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:507
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:313
bool IsWeighted() const
return true if the data set is weighted We cannot compute ourselfs because sometimes errors are fille...
Definition: BinData.h:572
double Value(unsigned int ipoint) const
return the value for the given fit point
Definition: BinData.h:216
ErrorType GetErrorType() const
retrieve the errortype
Definition: BinData.h:550
bool HaveCoordErrors() const
flag to control if data provides error on the coordinates
Definition: BinData.h:130
const double * ValuePtr(unsigned int ipoint) const
return a pointer to the value for the given fit point
Definition: BinData.h:228
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:520
double Error(unsigned int ipoint) const
Definition: BinData.h:250
bool HaveAsymErrors() const
flag to control if data provides asymmetric errors on the value
Definition: BinData.h:142
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:565
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:71
void GetRange(unsigned int irange, unsigned int icoord, double &xmin, double &xmax) const
get the i-th range for given coordinate.
Definition: DataRange.h:104
unsigned int Size() const
return number of fit points
Definition: FitData.h:303
const double * GetCoordComponent(unsigned int ipoint, unsigned int icoord) const
returns a single coordinate component of a point.
Definition: FitData.h:229
unsigned int NDim() const
return coordinate data dimension
Definition: FitData.h:311
const DataOptions & Opt() const
access to options
Definition: FitData.h:319
const double * Coords(unsigned int ipoint) const
return a pointer to the coordinates data for the given fit point
Definition: FitData.h:246
const DataRange & Range() const
access to range
Definition: FitData.h:331
SimpleGradientCalculator(int gdim, const IModelFunction &func, double eps=2.E-8, int istrat=1)
Definition: FitUtil.cxx:77
double ParameterDerivative(const double *x, const double *p, int ipar) const
Definition: FitUtil.cxx:116
void Gradient(const double *x, const double *p, double f0, double *g)
Definition: FitUtil.cxx:133
void ParameterGradient(const double *x, const double *p, double f0, double *g)
Definition: FitUtil.cxx:124
double DoParameterDerivative(const double *x, const double *p, double f0, int k) const
Definition: FitUtil.cxx:88
Class describing the unbinned data sets (just x coordinates values) of any dimensions.
Definition: UnBinData.h:42
double Weight(unsigned int ipoint) const
return weight
Definition: UnBinData.h:257
unsigned int NDim() const
return coordinate data dimension
Definition: UnBinData.h:280
virtual unsigned int NDim() const =0
Retrieve the dimension of the function.
virtual void SetParameters(const double *p)=0
Set the parameter values.
virtual unsigned int NPar() const =0
Return the number of Parameters.
Interface (abstract class) for parametric gradient multi-dimensional functions providing in addition ...
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 ...
OneDimMultiFunctionAdapter class to wrap a multidimensional function in one dimensional one.
User class for calculating the derivatives of a function.
double Derivative1(double x)
Returns the first derivative of the function at point x, computed by Richardson's extrapolation metho...
A pseudo container class which is a generator of indices.
Definition: TSeq.hxx:66
This class provides a simple interface to execute the same task multiple times in parallel threads,...
auto MapReduce(F func, unsigned nTimes, R redfunc) -> typename std::result_of< F()>::type
Execute a function nTimes in parallel (Map) and accumulate the results into a single value (Reduce).
Type
enumeration specifying the integration types.
@ kGAUSS
simple Gauss integration method with fixed rule
@ kDEFAULT
default type specified in the static options
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
Double_t ey[n]
Definition: legend1.C:17
Double_t ex[n]
Definition: legend1.C:17
TF1 * f1
Definition: legend1.C:11
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
ROOT::Math::IParamMultiGradFunction IGradModelFunction
Definition: FitUtil.h:65
double EvaluatePdf(const IModelFunction &func, const UnBinData &data, const double *p, unsigned int ipoint, double *g=0)
evaluate the pdf contribution to the LogL given a model function and the BinPoint data.
Definition: FitUtil.cxx:867
ROOT::Math::IParamMultiFunction IModelFunction
Definition: FitUtil.h:64
void CalculateGradientIntegral(const GFunc &gfunc, const double *x1, const double *x2, const double *p, double *g)
Definition: FitUtil.cxx:202
double CorrectValue(double rval)
Definition: FitUtil.cxx:167
void EvaluatePoissonLogLGradient(const IModelFunction &func, const BinData &data, const double *p, double *grad, unsigned int &nPoints, ::ROOT::EExecutionPolicy executionPolicy=::ROOT::EExecutionPolicy::kSequential, unsigned nChunks=0)
evaluate the Poisson LogL given a model function and the data at the point p.
double EvaluatePoissonLogL(const IModelFunction &func, const BinData &data, const double *p, int iWeight, bool extended, unsigned int &nPoints, ::ROOT::EExecutionPolicy executionPolicy, unsigned nChunks=0)
evaluate the Poisson LogL given a model function and the data at the point p.
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:422
void EvaluateLogLGradient(const IModelFunction &func, const UnBinData &data, const double *p, double *grad, unsigned int &nPoints, ::ROOT::EExecutionPolicy executionPolicy=::ROOT::EExecutionPolicy::kSequential, unsigned nChunks=0)
evaluate the LogL gradient given a model function and the data at the point p.
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
bool CheckInfNaNValue(double &rval)
Definition: FitUtil.cxx:181
unsigned setAutomaticChunking(unsigned nEvents)
Definition: FitUtil.cxx:1784
double EvaluateLogL(const IModelFunction &func, const UnBinData &data, const double *p, int iWeight, bool extended, unsigned int &nPoints, ::ROOT::EExecutionPolicy executionPolicy, unsigned nChunks=0)
evaluate the LogL given a model function and the data at the point x.
double EvaluateChi2(const IModelFunction &func, const BinData &data, const double *p, unsigned int &nPoints, ::ROOT::EExecutionPolicy executionPolicy, unsigned nChunks=0)
Chi2 Functions.
Definition: FitUtil.cxx:226
double EvaluateChi2Residual(const IModelFunction &func, const BinData &data, const double *p, 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:543
void EvaluateChi2Gradient(const IModelFunction &func, const BinData &data, const double *p, double *grad, unsigned int &nPoints, ::ROOT::EExecutionPolicy executionPolicy=::ROOT::EExecutionPolicy::kSequential, unsigned nChunks=0)
evaluate the Chi2 gradient given a model function and the data at the point p.
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
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.
UInt_t GetThreadPoolSize()
Returns the size of ROOT's thread pool.
Definition: TROOT.cxx:565
const char * Size
Definition: TXMLSetup.cxx:56
DataOptions : simple structure holding the options on how the data are filled.
Definition: DataOptions.h:28
bool fErrors1
use the function range when creating the fit data (default is false)
Definition: DataOptions.h:52
bool fNormBinVolume
normalize data by the bin volume (it is used in the Poisson likelihood fits)
Definition: DataOptions.h:49
bool fUseEmpty
normalize data by a normalized the bin volume (bin volume divided by a reference value)
Definition: DataOptions.h:50
bool fExpErrors
use all errors equal to 1, i.e. fit without errors (default is false)
Definition: DataOptions.h:53
bool fBinVolume
use integral of bin content instead of bin center (default is false)
Definition: DataOptions.h:48
bool fCoordErrors
use expected errors from the function and not from the data
Definition: DataOptions.h:54
double operator()(const double *x, const double *p) const
Definition: FitUtil.cxx:57
void SetDerivComponent(unsigned int ipar)
Definition: FitUtil.cxx:56
unsigned int NDim() const
Definition: FitUtil.cxx:60
ParamDerivFunc(const GradFunc &f)
Definition: FitUtil.cxx:55
auto * l
Definition: textangle.C:4