Logo ROOT  

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