Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
FitUtil.h
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// Header file for class FitUtil
12
13#ifndef ROOT_Fit_FitUtil
14#define ROOT_Fit_FitUtil
15
17#include "Math/IParamFunction.h"
18
19#ifdef R__USE_IMT
21#endif
23
24#include "Fit/BinData.h"
25#include "Fit/UnBinData.h"
27
28#include "Math/Integrator.h"
30
31#include "TError.h"
32#include <vector>
33
34// using parameter cache is not thread safe but needed for normalizing the functions
35#define USE_PARAMCACHE
36
37//#define DEBUG_FITUTIL
38
39#ifdef R__HAS_VECCORE
40namespace vecCore {
41template <class T>
42vecCore::Mask<T> Int2Mask(unsigned i)
43{
44 T x;
45 for (unsigned j = 0; j < vecCore::VectorSize<T>(); j++)
46 vecCore::Set<T>(x, j, j);
47 return vecCore::Mask<T>(x < T(i));
48}
49}
50#endif
51
52namespace ROOT {
53
54 namespace Fit {
55
56/**
57 namespace defining utility free functions using in Fit for evaluating the various fit method
58 functions (chi2, likelihood, etc..) given the data and the model function
59
60 @ingroup FitMain
61*/
62namespace FitUtil {
63
66
67 template <class T>
69
70 template <class T>
72
73 // internal class defining
74 template <class T>
76 public:
77 LikelihoodAux(T logv = {}, T w = {}, T w2 = {}) : logvalue(logv), weight(w), weight2(w2) {}
78
80 {
81 return LikelihoodAux<T>(logvalue + l.logvalue, weight + l.weight, weight2 + l.weight2);
82 }
83
85 {
86 logvalue += l.logvalue;
87 weight += l.weight;
88 weight2 += l.weight2;
89 return *this;
90 }
91
95 };
96
97 template <>
99 public:
100 LikelihoodAux(double logv = 0.0, double w = 0.0, double w2 = 0.0) : logvalue(logv), weight(w), weight2(w2){};
101
103 {
104 return LikelihoodAux<double>(logvalue + l.logvalue, weight + l.weight, weight2 + l.weight2);
105 }
106
108 {
109 logvalue += l.logvalue;
110 weight += l.weight;
111 weight2 += l.weight2;
112 return *this;
113 }
114
115 double logvalue;
116 double weight;
117 double weight2;
118 };
119
120 // internal class to evaluate the function or the integral
121 // and cached internal integration details
122 // if useIntegral is false no allocation is done
123 // and this is a dummy class
124 // class is templated on any parametric functor implementing operator()(x,p) and NDim()
125 // contains a constant pointer to the function
126
127 template <class ParamFunc = ROOT::Math::IParamMultiFunctionTempl<double>>
129
130 public:
131 IntegralEvaluator(const ParamFunc &func, const double *p, bool useIntegral = true,
133 : fDim(0), fParams(nullptr), fFunc(nullptr), fIg1Dim(nullptr), fIgNDim(nullptr), fFunc1Dim(nullptr), fFuncNDim(nullptr)
134 {
135 if (useIntegral) {
136 SetFunction(func, p, igType);
137 }
138 }
139
140 void SetFunction(const ParamFunc &func, const double *p = nullptr,
142 {
143 // set the integrand function and create required wrapper
144 // to perform integral in (x) of a generic f(x,p)
145 fParams = p;
146 fDim = func.NDim();
147 // copy the function object to be able to modify the parameters
148 // fFunc = dynamic_cast<ROOT::Math::IParamMultiFunction *>( func.Clone() );
149 fFunc = &func;
150 assert(fFunc != nullptr);
151 // set parameters in function
152 // fFunc->SetParameters(p);
153 if (fDim == 1) {
154 fFunc1Dim =
156 *this, &IntegralEvaluator::F1);
158 // fIg1Dim->SetFunction( static_cast<const ROOT::Math::IMultiGenFunction & >(*fFunc),false);
159 fIg1Dim->SetFunction(static_cast<const ROOT::Math::IGenFunction &>(*fFunc1Dim));
160 } else if (fDim > 1) {
161 fFuncNDim =
163 const>(*this, &IntegralEvaluator::FN, fDim);
166 } else
167 assert(fDim > 0);
168 }
169
170 void SetParameters(const double *p)
171 {
172 // copy just the pointer
173 fParams = p;
174 }
175
177 {
178 if (fIg1Dim)
179 delete fIg1Dim;
180 if (fIgNDim)
181 delete fIgNDim;
182 if (fFunc1Dim)
183 delete fFunc1Dim;
184 if (fFuncNDim)
185 delete fFuncNDim;
186 // if (fFunc) delete fFunc;
187 }
188
189 // evaluation of integrand function (one-dim)
190 double F1(double x) const
191 {
192 double xx = x;
193 return ExecFunc(fFunc, &xx, fParams);
194 }
195 // evaluation of integrand function (multi-dim)
196 double FN(const double *x) const { return ExecFunc(fFunc, x, fParams); }
197
198 double Integral(const double *x1, const double *x2)
199 {
200 // return unnormalized integral
201 return (fIg1Dim) ? fIg1Dim->Integral(*x1, *x2) : fIgNDim->Integral(x1, x2);
202 }
203
204 double operator()(const double *x1, const double *x2)
205 {
206 // return normalized integral, divided by bin volume (dx1*dx...*dxn)
207 if (fIg1Dim) {
208 double dV = *x2 - *x1;
209 return fIg1Dim->Integral(*x1, *x2) / dV;
210 } else if (fIgNDim) {
211 double dV = 1;
212 for (unsigned int i = 0; i < fDim; ++i)
213 dV *= (x2[i] - x1[i]);
214 return fIgNDim->Integral(x1, x2) / dV;
215 // std::cout << " do integral btw x " << x1[0] << " " << x2[0] << " y " << x1[1] << " "
216 // << x2[1] << " dV = " << dV << " result = " << result << std::endl; return result;
217 } else
218 assert(1.); // should never be here
219 return 0;
220 }
221
222 private:
223 template <class T>
224 inline double ExecFunc(T *f, const double *x, const double *p) const
225 {
226 return (*f)(x, p);
227 }
228
229#ifdef R__HAS_VECCORE
230 inline double ExecFunc(const IModelFunctionTempl<ROOT::Double_v> *f, const double *x, const double *p) const
231 {
232 // Figure out the size of the SIMD vectors.
233 constexpr static int vecSize = sizeof(ROOT::Double_v) / sizeof(double);
234 double xBuffer[vecSize];
236 for (unsigned int i = 0; i < fDim; ++i) {
237 // The Load() function reads multiple values from the pointed-to
238 // memory into xx. This is why we have to copy the input values from
239 // the x array into a zero-padded buffer to read from. Otherwise,
240 // Load() would access the x array out of bounds.
241 *xBuffer = x[i];
242 for(int j = 1; j < vecSize; ++j) {
243 xBuffer[j] = 0.0;
244 }
245 vecCore::Load<ROOT::Double_v>(xx[i], xBuffer);
246 }
247 auto res = (*f)(xx, p);
248 return vecCore::Get<ROOT::Double_v>(res, 0);
249 }
250#endif
251
252 // objects of this class are not meant to be copied / assigned
255
256 unsigned int fDim;
257 const double *fParams;
258 // ROOT::Math::IParamMultiFunction * fFunc; // copy of function in order to be able to change parameters
259 // const ParamFunc * fFunc; // reference to a generic parametric function
260 const ParamFunc *fFunc;
265 };
266
267 /** Chi2 Functions */
268
269 /**
270 evaluate the Chi2 given a model function and the data at the point x.
271 return also nPoints as the effective number of used points in the Chi2 evaluation
272 */
273 double EvaluateChi2(const IModelFunction &func, const BinData &data, const double *p, unsigned int &nPoints,
274 ::ROOT::EExecutionPolicy executionPolicy, unsigned nChunks = 0);
275
276 /**
277 evaluate the effective Chi2 given a model function and the data at the point x.
278 The effective chi2 uses the errors on the coordinates : W = 1/(sigma_y**2 + ( sigma_x_i * df/dx_i )**2 )
279 return also nPoints as the effective number of used points in the Chi2 evaluation
280 */
281 double EvaluateChi2Effective(const IModelFunction &func, const BinData &data, const double *x, unsigned int &nPoints);
282
283 /**
284 evaluate the Chi2 gradient given a model function and the data at the point p.
285 return also nPoints as the effective number of used points in the Chi2 evaluation
286 */
287 void EvaluateChi2Gradient(const IModelFunction &func, const BinData &data, const double *p, double *grad,
288 unsigned int &nPoints,
290 unsigned nChunks = 0);
291
292 /**
293 evaluate the LogL given a model function and the data at the point x.
294 return also nPoints as the effective number of used points in the LogL evaluation
295 */
296 double EvaluateLogL(const IModelFunction &func, const UnBinData &data, const double *p, int iWeight, bool extended,
297 unsigned int &nPoints, ::ROOT::EExecutionPolicy executionPolicy, unsigned nChunks = 0);
298
299 /**
300 evaluate the LogL gradient given a model function and the data at the point p.
301 return also nPoints as the effective number of used points in the LogL evaluation
302 */
303 void EvaluateLogLGradient(const IModelFunction &func, const UnBinData &data, const double *p, double *grad,
304 unsigned int &nPoints,
306 unsigned nChunks = 0);
307
308 // #ifdef R__HAS_VECCORE
309 // template <class NotCompileIfScalarBackend = std::enable_if<!(std::is_same<double, ROOT::Double_v>::value)>>
310 // void EvaluateLogLGradient(const IModelFunctionTempl<ROOT::Double_v> &, const UnBinData &, const double *, double
311 // *, unsigned int & ) {}
312 // #endif
313
314 /**
315 evaluate the Poisson LogL given a model function and the data at the point p.
316 return also nPoints as the effective number of used points in the LogL evaluation
317 By default is extended, pass extend to false if want to be not extended (MultiNomial)
318 */
319 double EvaluatePoissonLogL(const IModelFunction &func, const BinData &data, const double *p, int iWeight,
320 bool extended, unsigned int &nPoints, ::ROOT::EExecutionPolicy executionPolicy,
321 unsigned nChunks = 0);
322
323 /**
324 evaluate the Poisson LogL given a model function and the data at the point p.
325 return also nPoints as the effective number of used points in the LogL evaluation
326 */
327 void EvaluatePoissonLogLGradient(const IModelFunction &func, const BinData &data, const double *p, double *grad,
328 unsigned int &nPoints,
330 unsigned nChunks = 0);
331
332 // methods required by dedicate minimizer like Fumili
333
334 /**
335 evaluate the residual contribution to the Chi2 given a model function and the BinPoint data
336 and if the pointer g is not null evaluate also the gradient of the residual.
337 If the function provides parameter derivatives they are used otherwise a simple derivative calculation
338 is used
339 */
340 double EvaluateChi2Residual(const IModelFunction &func, const BinData &data, const double *p, unsigned int ipoint,
341 double *g = nullptr, double * h = nullptr, bool hasGrad = false, bool fullHessian = false);
342
343 /**
344 evaluate the pdf contribution to the LogL given a model function and the BinPoint data.
345 If the pointer g is not null evaluate also the gradient of the pdf.
346 If the function provides parameter derivatives they are used otherwise a simple derivative calculation
347 is used
348 */
349 double
350 EvaluatePdf(const IModelFunction &func, const UnBinData &data, const double *p, unsigned int ipoint, double *g = nullptr, double * h = nullptr, bool hasGrad = false, bool fullHessian = false);
351
352
353 /**
354 evaluate the pdf contribution to the Poisson LogL given a model function and the BinPoint data.
355 If the pointer g is not null evaluate also the gradient of the Poisson pdf.
356 If the function provides parameter derivatives they are used otherwise a simple derivative calculation
357 is used
358 */
359 double EvaluatePoissonBinPdf(const IModelFunction & func, const BinData & data, const double * x, unsigned int ipoint, double * g = nullptr, double * h = nullptr, bool hasGrad = false, bool fullHessian = false);
360
361 unsigned setAutomaticChunking(unsigned nEvents);
362
363 template<class T>
364 struct Evaluate {
365#ifdef R__HAS_VECCORE
366 static double EvalChi2(const IModelFunctionTempl<T> &func, const BinData &data, const double *p,
367 unsigned int &nPoints, ::ROOT::EExecutionPolicy executionPolicy, unsigned nChunks = 0)
368 {
369 // evaluate the chi2 given a vectorized function reference , the data and returns the value and also in nPoints
370 // the actual number of used points
371 // normal chi2 using only error on values (from fitting histogram)
372 // optionally the integral of function in the bin is used
373
374 //Info("EvalChi2","Using vectorized implementation %d",(int) data.Opt().fIntegral);
375
376 unsigned int n = data.Size();
377 nPoints = data.Size(); // npoints
378
379 // set parameters of the function to cache integral value
380#ifdef USE_PARAMCACHE
381 (const_cast<IModelFunctionTempl<T> &>(func)).SetParameters(p);
382#endif
383 // do not cache parameter values (it is not thread safe)
384 //func.SetParameters(p);
385
386
387 // get fit option and check case if using integral of bins
388 const DataOptions &fitOpt = data.Opt();
389 if (fitOpt.fBinVolume || fitOpt.fIntegral || fitOpt.fExpErrors)
390 Error("FitUtil::EvaluateChi2", "The vectorized implementation doesn't support Integrals, BinVolume or ExpErrors\n. Aborting operation.");
391
392 (const_cast<IModelFunctionTempl<T> &>(func)).SetParameters(p);
393
394 double maxResValue = std::numeric_limits<double>::max() / n;
395 std::vector<double> ones{1., 1., 1., 1.};
396 auto vecSize = vecCore::VectorSize<T>();
397
398 auto mapFunction = [&](unsigned int i) {
399 // in case of no error in y invError=1 is returned
400 T x1, y, invErrorVec;
401 vecCore::Load<T>(x1, data.GetCoordComponent(i * vecSize, 0));
402 vecCore::Load<T>(y, data.ValuePtr(i * vecSize));
403 const auto invError = data.ErrorPtr(i * vecSize);
404 auto invErrorptr = (invError != nullptr) ? invError : &ones.front();
405 vecCore::Load<T>(invErrorVec, invErrorptr);
406
407 const T *x;
408 std::vector<T> xc;
409 if(data.NDim() > 1) {
410 xc.resize(data.NDim());
411 xc[0] = x1;
412 for (unsigned int j = 1; j < data.NDim(); ++j)
413 vecCore::Load<T>(xc[j], data.GetCoordComponent(i * vecSize, j));
414 x = xc.data();
415 } else {
416 x = &x1;
417 }
418
419 T fval{};
420
421#ifdef USE_PARAMCACHE
422 fval = func(x);
423#else
424 fval = func(x, p);
425#endif
426
427 T tmp = (y - fval) * invErrorVec;
428 T chi2 = tmp * tmp;
429
430
431 // avoid infinity or nan in chi2 values due to wrong function values
432 auto m = vecCore::Mask_v<T>(chi2 > maxResValue);
433
434 vecCore::MaskedAssign<T>(chi2, m, maxResValue);
435
436 return chi2;
437 };
438
439 auto redFunction = [](const std::vector<T> &objs) {
440 return std::accumulate(objs.begin(), objs.end(), T{});
441 };
442
443#ifndef R__USE_IMT
444 (void)nChunks;
445
446 // If IMT is disabled, force the execution policy to the serial case
447 if (executionPolicy == ::ROOT::EExecutionPolicy::kMultiThread) {
448 Warning("FitUtil::EvaluateChi2", "Multithread execution policy requires IMT, which is disabled. Changing "
449 "to ::ROOT::EExecutionPolicy::kSequential.");
451 }
452#endif
453
454 T res{};
455 if (executionPolicy == ::ROOT::EExecutionPolicy::kSequential) {
457 res = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, data.Size()/vecSize), redFunction);
458#ifdef R__USE_IMT
459 } else if (executionPolicy == ::ROOT::EExecutionPolicy::kMultiThread) {
461 auto chunks = nChunks != 0 ? nChunks : setAutomaticChunking(data.Size() / vecSize);
462 res = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, data.Size() / vecSize), redFunction, chunks);
463#endif
464 } else {
465 Error("FitUtil::EvaluateChi2", "Execution policy unknown. Available choices:\n ::ROOT::EExecutionPolicy::kSequential (default)\n ::ROOT::EExecutionPolicy::kMultiThread (requires IMT)\n");
466 }
467
468 // Last SIMD vector of elements (if padding needed)
469 if (data.Size() % vecSize != 0)
470 vecCore::MaskedAssign(res, vecCore::Int2Mask<T>(data.Size() % vecSize),
471 res + mapFunction(data.Size() / vecSize));
472
473 return vecCore::ReduceAdd(res);
474 }
475
476 static double EvalLogL(const IModelFunctionTempl<T> &func, const UnBinData &data, const double *const p,
477 int iWeight, bool extended, unsigned int &nPoints,
478 ::ROOT::EExecutionPolicy executionPolicy, unsigned nChunks = 0)
479 {
480 // evaluate the LogLikelihood
481 unsigned int n = data.Size();
482 nPoints = data.Size(); // npoints
483
484 //unsigned int nRejected = 0;
485 bool normalizeFunc = false;
486
487 // set parameters of the function to cache integral value
488#ifdef USE_PARAMCACHE
489 (const_cast<IModelFunctionTempl<T> &>(func)).SetParameters(p);
490#endif
491
492#ifdef R__USE_IMT
493 // in case parameter needs to be propagated to user function use trick to set parameters by calling one time the function
494 // this will be done in sequential mode and parameters can be set in a thread safe manner
495 if (!normalizeFunc) {
496 if (data.NDim() == 1) {
497 T x;
498 vecCore::Load<T>(x, data.GetCoordComponent(0, 0));
499 func( &x, p);
500 }
501 else {
502 std::vector<T> x(data.NDim());
503 for (unsigned int j = 0; j < data.NDim(); ++j)
504 vecCore::Load<T>(x[j], data.GetCoordComponent(0, j));
505 func( x.data(), p);
506 }
507 }
508#endif
509
510 // this is needed if function must be normalized
511 double norm = 1.0;
512 if (normalizeFunc) {
513 // compute integral of the function
514 std::vector<double> xmin(data.NDim());
515 std::vector<double> xmax(data.NDim());
516 IntegralEvaluator<IModelFunctionTempl<T>> igEval(func, p, true);
517 // compute integral in the ranges where is defined
518 if (data.Range().Size() > 0) {
519 norm = 0;
520 for (unsigned int ir = 0; ir < data.Range().Size(); ++ir) {
521 data.Range().GetRange(&xmin[0], &xmax[0], ir);
522 norm += igEval.Integral(xmin.data(), xmax.data());
523 }
524 } else {
525 // use (-inf +inf)
526 data.Range().GetRange(&xmin[0], &xmax[0]);
527 // check if function is zero at +- inf
528 T xmin_v, xmax_v;
529 vecCore::Load<T>(xmin_v, xmin.data());
530 vecCore::Load<T>(xmax_v, xmax.data());
531 if (vecCore::ReduceAdd(func(&xmin_v, p)) != 0 || vecCore::ReduceAdd(func(&xmax_v, p)) != 0) {
532 MATH_ERROR_MSG("FitUtil::EvaluateLogLikelihood", "A range has not been set and the function is not zero at +/- inf");
533 return 0;
534 }
535 norm = igEval.Integral(&xmin[0], &xmax[0]);
536 }
537 }
538
539 // needed to compute effective global weight in case of extended likelihood
540
541 auto vecSize = vecCore::VectorSize<T>();
542 unsigned int numVectors = n / vecSize;
543
544 auto mapFunction = [ &, p](const unsigned i) {
545 T W{};
546 T W2{};
547 T fval{};
548
549 (void)p; /* avoid unused lambda capture warning if PARAMCACHE is disabled */
550
551 T x1;
552 vecCore::Load<T>(x1, data.GetCoordComponent(i * vecSize, 0));
553 const T *x = nullptr;
554 unsigned int ndim = data.NDim();
555 std::vector<T> xc;
556 if (ndim > 1) {
557 xc.resize(ndim);
558 xc[0] = x1;
559 for (unsigned int j = 1; j < ndim; ++j)
560 vecCore::Load<T>(xc[j], data.GetCoordComponent(i * vecSize, j));
561 x = xc.data();
562 } else {
563 x = &x1;
564 }
565
566#ifdef USE_PARAMCACHE
567 fval = func(x);
568#else
569 fval = func(x, p);
570#endif
571
572#ifdef DEBUG_FITUTIL
573 if (i < 5 || (i > numVectors-5) ) {
574 if (ndim == 1) std::cout << i << " x " << x[0] << " fval = " << fval;
575 else std::cout << i << " x " << x[0] << " y " << x[1] << " fval = " << fval;
576 }
577#endif
578
579 if (normalizeFunc) fval = fval * (1 / norm);
580
581 // function EvalLog protects against negative or too small values of fval
582 auto logval = ROOT::Math::Util::EvalLog(fval);
583 if (iWeight > 0) {
584 T weight{};
585 if (data.WeightsPtr(i) == nullptr)
586 weight = 1;
587 else
588 vecCore::Load<T>(weight, data.WeightsPtr(i*vecSize));
589 logval *= weight;
590 if (iWeight == 2) {
591 logval *= weight; // use square of weights in likelihood
592 if (!extended) {
593 // needed sum of weights and sum of weight square if likelkihood is extended
594 W = weight;
595 W2 = weight * weight;
596 }
597 }
598 }
599#ifdef DEBUG_FITUTIL
600 if (i < 5 || (i > numVectors-5)) {
601 std::cout << " " << fval << " logfval " << logval << std::endl;
602 }
603#endif
604
605 return LikelihoodAux<T>(logval, W, W2);
606 };
607
608 auto redFunction = [](const std::vector<LikelihoodAux<T>> &objs) {
609 return std::accumulate(objs.begin(), objs.end(), LikelihoodAux<T>(),
610 [](const LikelihoodAux<T> &l1, const LikelihoodAux<T> &l2) {
611 return l1 + l2;
612 });
613 };
614
615#ifndef R__USE_IMT
616 (void)nChunks;
617
618 // If IMT is disabled, force the execution policy to the serial case
619 if (executionPolicy == ::ROOT::EExecutionPolicy::kMultiThread) {
620 Warning("FitUtil::EvaluateLogL", "Multithread execution policy requires IMT, which is disabled. Changing "
621 "to ::ROOT::EExecutionPolicy::kSequential.");
623 }
624#endif
625
626 T logl_v{};
627 T sumW_v{};
628 T sumW2_v{};
630 if (executionPolicy == ::ROOT::EExecutionPolicy::kSequential) {
632 resArray = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, data.Size() / vecSize), redFunction);
633#ifdef R__USE_IMT
634 } else if (executionPolicy == ::ROOT::EExecutionPolicy::kMultiThread) {
636 auto chunks = nChunks != 0 ? nChunks : setAutomaticChunking( numVectors);
637 resArray = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, data.Size() / vecSize), redFunction, chunks);
638#endif
639 } else {
640 Error("FitUtil::EvaluateLogL", "Execution policy unknown. Available choices:\n ::ROOT::EExecutionPolicy::kSequential (default)\n ::ROOT::EExecutionPolicy::kMultiThread (requires IMT)\n");
641 }
642
643 logl_v = resArray.logvalue;
644 sumW_v = resArray.weight;
645 sumW2_v = resArray.weight2;
646
647 // Compute the contribution from the remaining points ( Last padded SIMD vector of elements )
648 unsigned int remainingPoints = n % vecSize;
649 if (remainingPoints > 0) {
650 auto remainingPointsContribution = mapFunction(numVectors);
651 // Add the contribution from the valid remaining points and store the result in the output variable
652 auto remainingMask = vecCore::Int2Mask<T>(remainingPoints);
653 vecCore::MaskedAssign(logl_v, remainingMask, logl_v + remainingPointsContribution.logvalue);
654 vecCore::MaskedAssign(sumW_v, remainingMask, sumW_v + remainingPointsContribution.weight);
655 vecCore::MaskedAssign(sumW2_v, remainingMask, sumW2_v + remainingPointsContribution.weight2);
656 }
657
658
659 //reduce vector type to double.
660 double logl = vecCore::ReduceAdd(logl_v);
661 double sumW = vecCore::ReduceAdd(sumW_v);
662 double sumW2 = vecCore::ReduceAdd(sumW2_v);
663
664 if (extended) {
665 // add Poisson extended term
666 double extendedTerm = 0; // extended term in likelihood
667 double nuTot = 0;
668 // nuTot is integral of function in the range
669 // if function has been normalized integral has been already computed
670 if (!normalizeFunc) {
671 IntegralEvaluator<IModelFunctionTempl<T>> igEval(func, p, true);
672 std::vector<double> xmin(data.NDim());
673 std::vector<double> xmax(data.NDim());
674
675 // compute integral in the ranges where is defined
676 if (data.Range().Size() > 0) {
677 nuTot = 0;
678 for (unsigned int ir = 0; ir < data.Range().Size(); ++ir) {
679 data.Range().GetRange(&xmin[0], &xmax[0], ir);
680 nuTot += igEval.Integral(xmin.data(), xmax.data());
681 }
682 } else {
683 // use (-inf +inf)
684 data.Range().GetRange(&xmin[0], &xmax[0]);
685 // check if function is zero at +- inf
686 T xmin_v, xmax_v;
687 vecCore::Load<T>(xmin_v, xmin.data());
688 vecCore::Load<T>(xmax_v, xmax.data());
689 if (vecCore::ReduceAdd(func(&xmin_v, p)) != 0 || vecCore::ReduceAdd(func(&xmax_v, p)) != 0) {
690 MATH_ERROR_MSG("FitUtil::EvaluateLogLikelihood", "A range has not been set and the function is not zero at +/- inf");
691 return 0;
692 }
693 nuTot = igEval.Integral(&xmin[0], &xmax[0]);
694 }
695
696 // force to be last parameter value
697 //nutot = p[func.NDim()-1];
698 if (iWeight != 2)
699 extendedTerm = - nuTot; // no need to add in this case n log(nu) since is already computed before
700 else {
701 // case use weight square in likelihood : compute total effective weight = sw2/sw
702 // ignore for the moment case when sumW is zero
703 extendedTerm = - (sumW2 / sumW) * nuTot;
704 }
705
706 } else {
707 nuTot = norm;
708 extendedTerm = - nuTot + double(n) * ROOT::Math::Util::EvalLog(nuTot);
709 // in case of weights need to use here sum of weights (to be done)
710 }
711
712 logl += extendedTerm;
713 }
714
715#ifdef DEBUG_FITUTIL
716 std::cout << "Evaluated log L for parameters (";
717 for (unsigned int ip = 0; ip < func.NPar(); ++ip)
718 std::cout << " " << p[ip];
719 std::cout << ") nll = " << -logl << std::endl;
720#endif
721
722 return -logl;
723
724 }
725
726 static double EvalPoissonLogL(const IModelFunctionTempl<T> &func, const BinData &data, const double *p,
727 int iWeight, bool extended, unsigned int,
728 ::ROOT::EExecutionPolicy executionPolicy, unsigned nChunks = 0)
729 {
730 // evaluate the Poisson Log Likelihood
731 // for binned likelihood fits
732 // this is Sum ( f(x_i) - y_i * log( f (x_i) ) )
733 // add as well constant term for saturated model to make it like a Chi2/2
734 // by default is extended. If extended is false the fit is not extended and
735 // the global poisson term is removed (i.e is a binomial fit)
736 // (remember that in this case one needs to have a function with a fixed normalization
737 // like in a non extended binned fit)
738 //
739 // if use Weight use a weighted dataset
740 // iWeight = 1 ==> logL = Sum( w f(x_i) )
741 // case of iWeight==1 is actually identical to weight==0
742 // iWeight = 2 ==> logL = Sum( w*w * f(x_i) )
743 //
744
745#ifdef USE_PARAMCACHE
746 (const_cast<IModelFunctionTempl<T> &>(func)).SetParameters(p);
747#endif
748 auto vecSize = vecCore::VectorSize<T>();
749 // get fit option and check case of using integral of bins
750 const DataOptions &fitOpt = data.Opt();
751 if (fitOpt.fExpErrors || fitOpt.fIntegral)
752 Error("FitUtil::EvaluateChi2",
753 "The vectorized implementation doesn't support Integrals or BinVolume\n. Aborting operation.");
754 bool useW2 = (iWeight == 2);
755
756 auto mapFunction = [&](unsigned int i) {
757 T y;
758 vecCore::Load<T>(y, data.ValuePtr(i * vecSize));
759 T fval{};
760
761 if (data.NDim() > 1) {
762 std::vector<T> x(data.NDim());
763 for (unsigned int j = 0; j < data.NDim(); ++j)
764 vecCore::Load<T>(x[j], data.GetCoordComponent(i * vecSize, j));
765#ifdef USE_PARAMCACHE
766 fval = func(x.data());
767#else
768 fval = func(x.data(), p);
769#endif
770 // one -dim case
771 } else {
772 T x;
773 vecCore::Load<T>(x, data.GetCoordComponent(i * vecSize, 0));
774#ifdef USE_PARAMCACHE
775 fval = func(&x);
776#else
777 fval = func(&x, p);
778#endif
779 }
780
781 // EvalLog protects against 0 values of fval but don't want to add in the -log sum
782 // negative values of fval
783 vecCore::MaskedAssign<T>(fval, fval < 0.0, 0.0);
784
785 T nloglike{}; // negative loglikelihood
786
787 if (useW2) {
788 // apply weight correction . Effective weight is error^2/ y
789 // and expected events in bins is fval/weight
790 // can apply correction only when y is not zero otherwise weight is undefined
791 // (in case of weighted likelihood I don't care about the constant term due to
792 // the saturated model)
793 assert (data.GetErrorType() != ROOT::Fit::BinData::ErrorType::kNoError);
794 T error = 0.0;
795 vecCore::Load<T>(error, data.ErrorPtr(i * vecSize));
796 // for empty bin use the average weight computed from the total data weight
797 auto m = vecCore::Mask_v<T>(y != 0.0);
798 auto weight = vecCore::Blend(m,(error * error) / y, T(data.SumOfError2()/ data.SumOfContent()) );
799 if (extended) {
800 nloglike = weight * ( fval - y);
801 }
802 vecCore::MaskedAssign<T>(nloglike, y != 0, nloglike + weight * y *( ROOT::Math::Util::EvalLog(y) - ROOT::Math::Util::EvalLog(fval)) );
803
804 } else {
805 // standard case no weights or iWeight=1
806 // this is needed for Poisson likelihood (which are extended and not for multinomial)
807 // the formula below include constant term due to likelihood of saturated model (f(x) = y)
808 // (same formula as in Baker-Cousins paper, page 439 except a factor of 2
809 if (extended) nloglike = fval - y;
810
811 vecCore::MaskedAssign<T>(
812 nloglike, y > 0, nloglike + y * (ROOT::Math::Util::EvalLog(y) - ROOT::Math::Util::EvalLog(fval)));
813 }
814
815 return nloglike;
816 };
817
818#ifdef R__USE_IMT
819 auto redFunction = [](const std::vector<T> &objs) { return std::accumulate(objs.begin(), objs.end(), T{}); };
820#else
821 (void)nChunks;
822
823 // If IMT is disabled, force the execution policy to the serial case
824 if (executionPolicy == ::ROOT::EExecutionPolicy::kMultiThread) {
825 Warning("FitUtil::Evaluate<T>::EvalPoissonLogL",
826 "Multithread execution policy requires IMT, which is disabled. Changing "
827 "to ::ROOT::EExecutionPolicy::kSequential.");
829 }
830#endif
831
832 T res{};
833 if (executionPolicy == ::ROOT::EExecutionPolicy::kSequential) {
834 for (unsigned int i = 0; i < (data.Size() / vecSize); i++) {
835 res += mapFunction(i);
836 }
837#ifdef R__USE_IMT
838 } else if (executionPolicy == ::ROOT::EExecutionPolicy::kMultiThread) {
840 auto chunks = nChunks != 0 ? nChunks : setAutomaticChunking(data.Size() / vecSize);
841 res = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, data.Size() / vecSize), redFunction, chunks);
842#endif
843 } else {
844 Error(
845 "FitUtil::Evaluate<T>::EvalPoissonLogL",
846 "Execution policy unknown. Available choices:\n ::ROOT::EExecutionPolicy::kSequential (default)\n ::ROOT::EExecutionPolicy::kMultiThread (requires IMT)\n");
847 }
848
849 // Last padded SIMD vector of elements
850 if (data.Size() % vecSize != 0)
851 vecCore::MaskedAssign(res, vecCore::Int2Mask<T>(data.Size() % vecSize),
852 res + mapFunction(data.Size() / vecSize));
853
854 return vecCore::ReduceAdd(res);
855 }
856
857 static double EvalChi2Effective(const IModelFunctionTempl<T> &, const BinData &, const double *, unsigned int &)
858 {
859 Error("FitUtil::Evaluate<T>::EvalChi2Effective", "The vectorized evaluation of the Chi2 with coordinate errors is still not supported");
860 return -1.;
861 }
862
863 // Compute a mask to filter out infinite numbers and NaN values.
864 // The argument rval is updated so infinite numbers and NaN values are replaced by
865 // maximum finite values (preserving the original sign).
866 static vecCore::Mask<T> CheckInfNaNValues(T &rval)
867 {
868 auto mask = rval > -vecCore::NumericLimits<T>::Max() && rval < vecCore::NumericLimits<T>::Max();
869
870 // Case +inf or nan
871 vecCore::MaskedAssign(rval, !mask, +vecCore::NumericLimits<T>::Max());
872
873 // Case -inf
874 vecCore::MaskedAssign(rval, !mask && rval < 0, -vecCore::NumericLimits<T>::Max());
875
876 return mask;
877 }
878
879 static void EvalChi2Gradient(const IModelFunctionTempl<T> &f, const BinData &data, const double *p, double *grad,
880 unsigned int &nPoints,
882 unsigned nChunks = 0)
883 {
884 // evaluate the gradient of the chi2 function
885 // this function is used when the model function knows how to calculate the derivative and we can
886 // avoid that the minimizer re-computes them
887 //
888 // case of chi2 effective (errors on coordinate) is not supported
889
890 if (data.HaveCoordErrors()) {
891 MATH_ERROR_MSG("FitUtil::EvaluateChi2Gradient",
892 "Error on the coordinates are not used in calculating Chi2 gradient");
893 return; // it will assert otherwise later in GetPoint
894 }
895
896 const IGradModelFunctionTempl<T> *fg = dynamic_cast<const IGradModelFunctionTempl<T> *>(&f);
897 assert(fg != nullptr); // must be called by a gradient function
898
899 const IGradModelFunctionTempl<T> &func = *fg;
900
901 const DataOptions &fitOpt = data.Opt();
902 if (fitOpt.fBinVolume || fitOpt.fIntegral || fitOpt.fExpErrors)
903 Error("FitUtil::EvaluateChi2Gradient", "The vectorized implementation doesn't support Integrals,"
904 "BinVolume or ExpErrors\n. Aborting operation.");
905
906 unsigned int npar = func.NPar();
907 auto vecSize = vecCore::VectorSize<T>();
908 unsigned initialNPoints = data.Size();
909 unsigned numVectors = initialNPoints / vecSize;
910
911 // numVectors + 1 because of the padded data (call to mapFunction with i = numVectors after the main loop)
912 std::vector<vecCore::Mask<T>> validPointsMasks(numVectors + 1);
913
914 auto mapFunction = [&](const unsigned int i) {
915 // set all vector values to zero
916 std::vector<T> gradFunc(npar);
917 std::vector<T> pointContributionVec(npar);
918
919 T x1, y, invError;
920
921 vecCore::Load<T>(x1, data.GetCoordComponent(i * vecSize, 0));
922 vecCore::Load<T>(y, data.ValuePtr(i * vecSize));
923 const auto invErrorPtr = data.ErrorPtr(i * vecSize);
924
925 if (invErrorPtr == nullptr)
926 invError = 1;
927 else
928 vecCore::Load<T>(invError, invErrorPtr);
929
930 // TODO: Check error options and invert if needed
931
932 T fval = 0;
933
934 const T *x = nullptr;
935
936 unsigned int ndim = data.NDim();
937 // need to declare vector outside if statement
938 // otherwise pointer will be invalid
939 std::vector<T> xc;
940 if (ndim > 1) {
941 xc.resize(ndim);
942 xc[0] = x1;
943 for (unsigned int j = 1; j < ndim; ++j)
944 vecCore::Load<T>(xc[j], data.GetCoordComponent(i * vecSize, j));
945 x = xc.data();
946 } else {
947 x = &x1;
948 }
949
950 fval = func(x, p);
951 func.ParameterGradient(x, p, &gradFunc[0]);
952
953 validPointsMasks[i] = CheckInfNaNValues(fval);
954 if (vecCore::MaskEmpty(validPointsMasks[i])) {
955 // Return a zero contribution to all partial derivatives on behalf of the current points
956 return pointContributionVec;
957 }
958
959 // loop on the parameters
960 for (unsigned int ipar = 0; ipar < npar; ++ipar) {
961 // avoid singularity in the function (infinity and nan ) in the chi2 sum
962 // eventually add possibility of excluding some points (like singularity)
963 validPointsMasks[i] = CheckInfNaNValues(gradFunc[ipar]);
964
965 if (vecCore::MaskEmpty(validPointsMasks[i])) {
966 break; // exit loop on parameters
967 }
968
969 // calculate derivative point contribution (only for valid points)
970 vecCore::MaskedAssign(pointContributionVec[ipar], validPointsMasks[i],
971 -2.0 * (y - fval) * invError * invError * gradFunc[ipar]);
972 }
973
974 return pointContributionVec;
975 };
976
977 // Reduce the set of vectors by summing its equally-indexed components
978 auto redFunction = [&](const std::vector<std::vector<T>> &partialResults) {
979 std::vector<T> result(npar);
980
981 for (auto const &pointContributionVec : partialResults) {
982 for (unsigned int parameterIndex = 0; parameterIndex < npar; parameterIndex++)
983 result[parameterIndex] += pointContributionVec[parameterIndex];
984 }
985
986 return result;
987 };
988
989 std::vector<T> gVec(npar);
990 std::vector<double> g(npar);
991
992#ifndef R__USE_IMT
993 // to fix compiler warning
994 (void)nChunks;
995
996 // If IMT is disabled, force the execution policy to the serial case
997 if (executionPolicy == ::ROOT::EExecutionPolicy::kMultiThread) {
998 Warning("FitUtil::EvaluateChi2Gradient",
999 "Multithread execution policy requires IMT, which is disabled. Changing "
1000 "to ::ROOT::EExecutionPolicy::kSequential.");
1001 executionPolicy = ::ROOT::EExecutionPolicy::kSequential;
1002 }
1003#endif
1004
1005 if (executionPolicy == ::ROOT::EExecutionPolicy::kSequential) {
1007 gVec = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, numVectors), redFunction);
1008 }
1009#ifdef R__USE_IMT
1010 else if (executionPolicy == ::ROOT::EExecutionPolicy::kMultiThread) {
1012 auto chunks = nChunks != 0 ? nChunks : setAutomaticChunking(numVectors);
1013 gVec = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, numVectors), redFunction, chunks);
1014 }
1015#endif
1016 else {
1017 Error(
1018 "FitUtil::EvaluateChi2Gradient",
1019 "Execution policy unknown. Available choices:\n 0: Serial (default)\n 1: MultiThread (requires IMT)\n");
1020 }
1021
1022 // Compute the contribution from the remaining points
1023 unsigned int remainingPoints = initialNPoints % vecSize;
1024 if (remainingPoints > 0) {
1025 auto remainingPointsContribution = mapFunction(numVectors);
1026 // Add the contribution from the valid remaining points and store the result in the output variable
1027 auto remainingMask = vecCore::Int2Mask<T>(remainingPoints);
1028 for (unsigned int param = 0; param < npar; param++) {
1029 vecCore::MaskedAssign(gVec[param], remainingMask, gVec[param] + remainingPointsContribution[param]);
1030 }
1031 }
1032 // reduce final gradient result from T to double
1033 for (unsigned int param = 0; param < npar; param++) {
1034 grad[param] = vecCore::ReduceAdd(gVec[param]);
1035 }
1036
1037 // correct the number of points
1038 nPoints = initialNPoints;
1039
1040 if (std::any_of(validPointsMasks.begin(), validPointsMasks.end(),
1041 [](vecCore::Mask<T> validPoints) { return !vecCore::MaskFull(validPoints); })) {
1042 unsigned nRejected = 0;
1043
1044 for (const auto &mask : validPointsMasks) {
1045 for (unsigned int i = 0; i < vecSize; i++) {
1046 nRejected += !vecCore::Get(mask, i);
1047 }
1048 }
1049
1050 assert(nRejected <= initialNPoints);
1051 nPoints = initialNPoints - nRejected;
1052
1053 if (nPoints < npar) {
1054 MATH_ERROR_MSG("FitUtil::EvaluateChi2Gradient",
1055 "Too many points rejected for overflow in gradient calculation");
1056 }
1057 }
1058 }
1059
1060 static double EvalChi2Residual(const IModelFunctionTempl<T> &, const BinData &, const double *, unsigned int, double *, double *, bool, bool)
1061 {
1062 Error("FitUtil::Evaluate<T>::EvalChi2Residual", "The vectorized evaluation of the Chi2 with the ith residual is still not supported");
1063 return -1.;
1064 }
1065
1066 /// evaluate the pdf (Poisson) contribution to the logl (return actually log of pdf)
1067 /// and its gradient
1068 static double EvalPoissonBinPdf(const IModelFunctionTempl<T> &, const BinData &, const double *, unsigned int , double * , double * , bool, bool) {
1069 Error("FitUtil::Evaluate<T>::EvaluatePoissonBinPdf", "The vectorized evaluation of the BinnedLikelihood fit evaluated point by point is still not supported");
1070 return -1.;
1071 }
1072
1073 static double EvalPdf(const IModelFunctionTempl<T> &, const UnBinData &, const double *, unsigned int , double * , double * , bool, bool) {
1074 Error("FitUtil::Evaluate<T>::EvalPdf", "The vectorized evaluation of the LogLikelihood fit evaluated point by point is still not supported");
1075 return -1.;
1076 }
1077
1078 //template <class NotCompileIfScalarBackend = std::enable_if<!(std::is_same<double, ROOT::Double_v>::value)>>
1079 // static double EvalPdf(const IModelFunctionTempl<ROOT::Double_v> &func, const UnBinData &data, const double *p, unsigned int i, double *, double *, bool, bool) {
1080 // // evaluate the pdf contribution to the generic logl function in case of bin data
1081 // // return actually the log of the pdf and its derivatives
1082 // // func.SetParameters(p);
1083 // const auto x = vecCore::FromPtr<ROOT::Double_v>(data.GetCoordComponent(i, 0));
1084 // auto fval = func(&x, p);
1085 // auto logPdf = ROOT::Math::Util::EvalLog(fval);
1086 // return vecCore::Get<ROOT::Double_v>(logPdf, 0);
1087
1088 static void
1089 EvalPoissonLogLGradient(const IModelFunctionTempl<T> &f, const BinData &data, const double *p, double *grad,
1090 unsigned int &,
1092 unsigned nChunks = 0)
1093 {
1094 // evaluate the gradient of the Poisson log likelihood function
1095
1096 const IGradModelFunctionTempl<T> *fg = dynamic_cast<const IGradModelFunctionTempl<T> *>(&f);
1097 assert(fg != nullptr); // must be called by a grad function
1098
1099 const IGradModelFunctionTempl<T> &func = *fg;
1100
1101 (const_cast<IGradModelFunctionTempl<T> &>(func)).SetParameters(p);
1102
1103
1104 const DataOptions &fitOpt = data.Opt();
1105 if (fitOpt.fBinVolume || fitOpt.fIntegral || fitOpt.fExpErrors)
1106 Error("FitUtil::EvaluatePoissonLogLGradient", "The vectorized implementation doesn't support Integrals,"
1107 "BinVolume or ExpErrors\n. Aborting operation.");
1108
1109 unsigned int npar = func.NPar();
1110 auto vecSize = vecCore::VectorSize<T>();
1111 unsigned initialNPoints = data.Size();
1112 unsigned numVectors = initialNPoints / vecSize;
1113
1114 auto mapFunction = [&](const unsigned int i) {
1115 // set all vector values to zero
1116 std::vector<T> gradFunc(npar);
1117 std::vector<T> pointContributionVec(npar);
1118
1119 T x1, y;
1120
1121 vecCore::Load<T>(x1, data.GetCoordComponent(i * vecSize, 0));
1122 vecCore::Load<T>(y, data.ValuePtr(i * vecSize));
1123
1124 T fval = 0;
1125
1126 const T *x = nullptr;
1127
1128 unsigned ndim = data.NDim();
1129 std::vector<T> xc;
1130 if (ndim > 1) {
1131 xc.resize(ndim);
1132 xc[0] = x1;
1133 for (unsigned int j = 1; j < ndim; ++j)
1134 vecCore::Load<T>(xc[j], data.GetCoordComponent(i * vecSize, j));
1135 x = xc.data();
1136 } else {
1137 x = &x1;
1138 }
1139
1140 fval = func(x, p);
1141 func.ParameterGradient(x, p, &gradFunc[0]);
1142
1143 // correct the gradient
1144 for (unsigned int ipar = 0; ipar < npar; ++ipar) {
1145 vecCore::Mask<T> positiveValuesMask = fval > 0;
1146
1147 // df/dp * (1. - y/f )
1148 vecCore::MaskedAssign(pointContributionVec[ipar], positiveValuesMask, gradFunc[ipar] * (1. - y / fval));
1149
1150 vecCore::Mask<T> validNegativeValuesMask = !positiveValuesMask && gradFunc[ipar] != 0;
1151
1152 if (!vecCore::MaskEmpty(validNegativeValuesMask)) {
1153 const T kdmax1 = vecCore::math::Sqrt(vecCore::NumericLimits<T>::Max());
1154 const T kdmax2 = vecCore::NumericLimits<T>::Max() / (4 * initialNPoints);
1155 T gg = kdmax1 * gradFunc[ipar];
1156 pointContributionVec[ipar] = -vecCore::Blend(gg > 0, vecCore::math::Min(gg, kdmax2), vecCore::math::Max(gg, -kdmax2));
1157 }
1158 }
1159
1160#ifdef DEBUG_FITUTIL
1161 {
1163 if (i < 5 || (i > data.Size()-5) ) {
1164 if (data.NDim() > 1) std::cout << i << " x " << x[0] << " y " << x[1];
1165 else std::cout << i << " x " << x[0];
1166 std::cout << " func " << fval << " gradient ";
1167 for (unsigned int ii = 0; ii < npar; ++ii) std::cout << " " << pointContributionVec[ii];
1168 std::cout << "\n";
1169 }
1170 }
1171#endif
1172
1173 return pointContributionVec;
1174 };
1175
1176 // Vertically reduce the set of vectors by summing its equally-indexed components
1177 auto redFunction = [&](const std::vector<std::vector<T>> &partialResults) {
1178 std::vector<T> result(npar);
1179
1180 for (auto const &pointContributionVec : partialResults) {
1181 for (unsigned int parameterIndex = 0; parameterIndex < npar; parameterIndex++)
1182 result[parameterIndex] += pointContributionVec[parameterIndex];
1183 }
1184
1185 return result;
1186 };
1187
1188 std::vector<T> gVec(npar);
1189
1190#ifndef R__USE_IMT
1191 // to fix compiler warning
1192 (void)nChunks;
1193
1194 // If IMT is disabled, force the execution policy to the serial case
1195 if (executionPolicy == ::ROOT::EExecutionPolicy::kMultiThread) {
1196 Warning("FitUtil::EvaluatePoissonLogLGradient",
1197 "Multithread execution policy requires IMT, which is disabled. Changing "
1198 "to ::ROOT::EExecutionPolicy::kSequential.");
1199 executionPolicy = ::ROOT::EExecutionPolicy::kSequential;
1200 }
1201#endif
1202
1203 if (executionPolicy == ::ROOT::EExecutionPolicy::kSequential) {
1205 gVec = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, numVectors), redFunction);
1206 }
1207#ifdef R__USE_IMT
1208 else if (executionPolicy == ::ROOT::EExecutionPolicy::kMultiThread) {
1210 auto chunks = nChunks != 0 ? nChunks : setAutomaticChunking(numVectors);
1211 gVec = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, numVectors), redFunction, chunks);
1212 }
1213#endif
1214 else {
1215 Error("FitUtil::EvaluatePoissonLogLGradient", "Execution policy unknown. Available choices:\n "
1216 "::ROOT::EExecutionPolicy::kSequential (default)\n "
1217 "::ROOT::EExecutionPolicy::kMultiThread (requires IMT)\n");
1218 }
1219
1220
1221 // Compute the contribution from the remaining points
1222 unsigned int remainingPoints = initialNPoints % vecSize;
1223 if (remainingPoints > 0) {
1224 auto remainingPointsContribution = mapFunction(numVectors);
1225 // Add the contribution from the valid remaining points and store the result in the output variable
1226 auto remainingMask = vecCore::Int2Mask<T>(remainingPoints);
1227 for (unsigned int param = 0; param < npar; param++) {
1228 vecCore::MaskedAssign(gVec[param], remainingMask, gVec[param] + remainingPointsContribution[param]);
1229 }
1230 }
1231 // reduce final gradient result from T to double
1232 for (unsigned int param = 0; param < npar; param++) {
1233 grad[param] = vecCore::ReduceAdd(gVec[param]);
1234 }
1235
1236#ifdef DEBUG_FITUTIL
1237 std::cout << "***** Final gradient : ";
1238 for (unsigned int ii = 0; ii< npar; ++ii) std::cout << grad[ii] << " ";
1239 std::cout << "\n";
1240#endif
1241
1242 }
1243
1244 static void EvalLogLGradient(const IModelFunctionTempl<T> &f, const UnBinData &data, const double *p,
1245 double *grad, unsigned int &,
1247 unsigned nChunks = 0)
1248 {
1249 // evaluate the gradient of the log likelihood function
1250
1251 const IGradModelFunctionTempl<T> *fg = dynamic_cast<const IGradModelFunctionTempl<T> *>(&f);
1252 assert(fg != nullptr); // must be called by a grad function
1253
1254 const IGradModelFunctionTempl<T> &func = *fg;
1255
1256
1257 unsigned int npar = func.NPar();
1258 auto vecSize = vecCore::VectorSize<T>();
1259 unsigned initialNPoints = data.Size();
1260 unsigned numVectors = initialNPoints / vecSize;
1261
1262#ifdef DEBUG_FITUTIL
1263 std::cout << "\n===> Evaluate Gradient for parameters ";
1264 for (unsigned int ip = 0; ip < npar; ++ip)
1265 std::cout << " " << p[ip];
1266 std::cout << "\n";
1267#endif
1268
1269 (const_cast<IGradModelFunctionTempl<T> &>(func)).SetParameters(p);
1270
1271 const T kdmax1 = vecCore::math::Sqrt(vecCore::NumericLimits<T>::Max());
1272 const T kdmax2 = vecCore::NumericLimits<T>::Max() / (4 * initialNPoints);
1273
1274 auto mapFunction = [&](const unsigned int i) {
1275 std::vector<T> gradFunc(npar);
1276 std::vector<T> pointContributionVec(npar);
1277
1278 T x1;
1279 vecCore::Load<T>(x1, data.GetCoordComponent(i * vecSize, 0));
1280
1281 const T *x = nullptr;
1282
1283 unsigned int ndim = data.NDim();
1284 std::vector<T> xc(ndim);
1285 if (ndim > 1) {
1286 xc.resize(ndim);
1287 xc[0] = x1;
1288 for (unsigned int j = 1; j < ndim; ++j)
1289 vecCore::Load<T>(xc[j], data.GetCoordComponent(i * vecSize, j));
1290 x = xc.data();
1291 } else {
1292 x = &x1;
1293 }
1294
1295
1296 T fval = func(x, p);
1297 func.ParameterGradient(x, p, &gradFunc[0]);
1298
1299#ifdef DEBUG_FITUTIL
1300 if (i < 5 || (i > numVectors-5) ) {
1301 if (ndim > 1) std::cout << i << " x " << x[0] << " y " << x[1] << " gradient " << gradFunc[0] << " " << gradFunc[1] << " " << gradFunc[3] << std::endl;
1302 else std::cout << i << " x " << x[0] << " gradient " << gradFunc[0] << " " << gradFunc[1] << " " << gradFunc[3] << std::endl;
1303 }
1304#endif
1305
1306 vecCore::Mask<T> positiveValues = fval > 0;
1307
1308 for (unsigned int kpar = 0; kpar < npar; ++kpar) {
1309 if (!vecCore::MaskEmpty(positiveValues))
1310 vecCore::MaskedAssign<T>(pointContributionVec[kpar], positiveValues, -1. / fval * gradFunc[kpar]);
1311
1312 vecCore::Mask<T> nonZeroGradientValues = !positiveValues && gradFunc[kpar] != 0;
1313 if (!vecCore::MaskEmpty(nonZeroGradientValues)) {
1314 T gg = kdmax1 * gradFunc[kpar];
1315 pointContributionVec[kpar] =
1316 vecCore::Blend(nonZeroGradientValues && gg > 0, -vecCore::math::Min(gg, kdmax2),
1317 -vecCore::math::Max(gg, -kdmax2));
1318 }
1319 // if func derivative is zero term is also zero so do not add in g[kpar]
1320 }
1321
1322 return pointContributionVec;
1323 };
1324
1325 // Vertically reduce the set of vectors by summing its equally-indexed components
1326 auto redFunction = [&](const std::vector<std::vector<T>> &pointContributions) {
1327 std::vector<T> result(npar);
1328
1329 for (auto const &pointContributionVec : pointContributions) {
1330 for (unsigned int parameterIndex = 0; parameterIndex < npar; parameterIndex++)
1331 result[parameterIndex] += pointContributionVec[parameterIndex];
1332 }
1333
1334 return result;
1335 };
1336
1337 std::vector<T> gVec(npar);
1338 std::vector<double> g(npar);
1339
1340#ifndef R__USE_IMT
1341 // to fix compiler warning
1342 (void)nChunks;
1343
1344 // If IMT is disabled, force the execution policy to the serial case
1345 if (executionPolicy == ::ROOT::EExecutionPolicy::kMultiThread) {
1346 Warning("FitUtil::EvaluateLogLGradient",
1347 "Multithread execution policy requires IMT, which is disabled. Changing "
1348 "to ::ROOT::EExecutionPolicy::kSequential.");
1349 executionPolicy = ::ROOT::EExecutionPolicy::kSequential;
1350 }
1351#endif
1352
1353 if (executionPolicy == ::ROOT::EExecutionPolicy::kSequential) {
1355 gVec = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, numVectors), redFunction);
1356 }
1357#ifdef R__USE_IMT
1358 else if (executionPolicy == ::ROOT::EExecutionPolicy::kMultiThread) {
1360 auto chunks = nChunks != 0 ? nChunks : setAutomaticChunking(numVectors);
1361 gVec = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, numVectors), redFunction, chunks);
1362 }
1363#endif
1364 else {
1365 Error("FitUtil::EvaluateLogLGradient", "Execution policy unknown. Available choices:\n "
1366 "::ROOT::EExecutionPolicy::kSequential (default)\n "
1367 "::ROOT::EExecutionPolicy::kMultiThread (requires IMT)\n");
1368 }
1369
1370 // Compute the contribution from the remaining points
1371 unsigned int remainingPoints = initialNPoints % vecSize;
1372 if (remainingPoints > 0) {
1373 auto remainingPointsContribution = mapFunction(numVectors);
1374 // Add the contribution from the valid remaining points and store the result in the output variable
1375 auto remainingMask = vecCore::Int2Mask<T>(initialNPoints % vecSize);
1376 for (unsigned int param = 0; param < npar; param++) {
1377 vecCore::MaskedAssign(gVec[param], remainingMask, gVec[param] + remainingPointsContribution[param]);
1378 }
1379 }
1380 // reduce final gradient result from T to double
1381 for (unsigned int param = 0; param < npar; param++) {
1382 grad[param] = vecCore::ReduceAdd(gVec[param]);
1383 }
1384
1385#ifdef DEBUG_FITUTIL
1386 std::cout << "Final gradient ";
1387 for (unsigned int param = 0; param < npar; param++) {
1388 std::cout << " " << grad[param];
1389 }
1390 std::cout << "\n";
1391#endif
1392 }
1393 };
1394
1395 template<>
1396 struct Evaluate<double>{
1397#endif
1398
1399 static double EvalChi2(const IModelFunction &func, const BinData &data, const double *p, unsigned int &nPoints,
1400 ::ROOT::EExecutionPolicy executionPolicy, unsigned nChunks = 0)
1401 {
1402 // evaluate the chi2 given a function reference, the data and returns the value and also in nPoints
1403 // the actual number of used points
1404 // normal chi2 using only error on values (from fitting histogram)
1405 // optionally the integral of function in the bin is used
1406
1407
1408 //Info("EvalChi2","Using non-vectorized implementation %d",(int) data.Opt().fIntegral);
1409
1410 return FitUtil::EvaluateChi2(func, data, p, nPoints, executionPolicy, nChunks);
1411 }
1412
1413 static double EvalLogL(const IModelFunctionTempl<double> &func, const UnBinData &data, const double *p,
1414 int iWeight, bool extended, unsigned int &nPoints,
1415 ::ROOT::EExecutionPolicy executionPolicy, unsigned nChunks = 0)
1416 {
1417 return FitUtil::EvaluateLogL(func, data, p, iWeight, extended, nPoints, executionPolicy, nChunks);
1418 }
1419
1420 static double EvalPoissonLogL(const IModelFunctionTempl<double> &func, const BinData &data, const double *p,
1421 int iWeight, bool extended, unsigned int &nPoints,
1422 ::ROOT::EExecutionPolicy executionPolicy, unsigned nChunks = 0)
1423 {
1424 return FitUtil::EvaluatePoissonLogL(func, data, p, iWeight, extended, nPoints, executionPolicy, nChunks);
1425 }
1426
1427 static double EvalChi2Effective(const IModelFunctionTempl<double> &func, const BinData & data, const double * p, unsigned int &nPoints)
1428 {
1429 return FitUtil::EvaluateChi2Effective(func, data, p, nPoints);
1430 }
1431 static void EvalChi2Gradient(const IModelFunctionTempl<double> &func, const BinData &data, const double *p,
1432 double *g, unsigned int &nPoints,
1434 unsigned nChunks = 0)
1435 {
1436 FitUtil::EvaluateChi2Gradient(func, data, p, g, nPoints, executionPolicy, nChunks);
1437 }
1438
1439 static double EvalChi2Residual(const IModelFunctionTempl<double> &func, const BinData & data, const double * p, unsigned int i, double *g, double * h,
1440 bool hasGrad, bool fullHessian)
1441 {
1442 return FitUtil::EvaluateChi2Residual(func, data, p, i, g, h, hasGrad, fullHessian);
1443 }
1444
1445 /// evaluate the pdf (Poisson) contribution to the logl (return actually log of pdf)
1446 /// and its gradient
1447 static double EvalPoissonBinPdf(const IModelFunctionTempl<double> &func, const BinData & data, const double *p, unsigned int i, double *g, double * h, bool hasGrad, bool fullHessian) {
1448 return FitUtil::EvaluatePoissonBinPdf(func, data, p, i, g, h, hasGrad, fullHessian);
1449 }
1450
1451 static double EvalPdf(const IModelFunctionTempl<double> &func, const UnBinData & data, const double *p, unsigned int i, double *g, double * h, bool hasGrad, bool fullHessian) {
1452 return FitUtil::EvaluatePdf(func, data, p, i, g, h, hasGrad, fullHessian);
1453 }
1454
1455 static void
1456 EvalPoissonLogLGradient(const IModelFunctionTempl<double> &func, const BinData &data, const double *p, double *g,
1457 unsigned int &nPoints,
1459 unsigned nChunks = 0)
1460 {
1461 FitUtil::EvaluatePoissonLogLGradient(func, data, p, g, nPoints, executionPolicy, nChunks);
1462 }
1463
1464 static void EvalLogLGradient(const IModelFunctionTempl<double> &func, const UnBinData &data, const double *p,
1465 double *g, unsigned int &nPoints,
1467 unsigned nChunks = 0)
1468 {
1469 FitUtil::EvaluateLogLGradient(func, data, p, g, nPoints, executionPolicy, nChunks);
1470 }
1471 };
1472
1473} // end namespace FitUtil
1474
1475 } // end namespace Fit
1476
1477} // end namespace ROOT
1478
1479#if defined (R__HAS_VECCORE) && defined(R__HAS_VC)
1480//Fixes alignment for structures of SIMD structures
1482#endif
1483
1484#endif /* ROOT_Fit_FitUtil */
#define MATH_ERROR_MSG(loc, str)
Definition Error.h:83
#define f(i)
Definition RSha256.hxx:104
#define g(i)
Definition RSha256.hxx:105
#define h(i)
Definition RSha256.hxx:106
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition TError.cxx:185
void Warning(const char *location, const char *msgfmt,...)
Use this function in warning situations.
Definition TError.cxx:229
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t mask
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
Option_t Option_t TPoint TPoint const char x2
Option_t Option_t TPoint TPoint const char x1
float xmin
float xmax
R__EXTERN TVirtualMutex * gROOTMutex
Definition TROOT.h:63
#define R__LOCKGUARD(mutex)
Class describing the binned data sets : vectors of x coordinates, y values and optionally error on y ...
Definition BinData.h:52
double Integral(const double *x1, const double *x2)
Definition FitUtil.h:198
ROOT::Math::IGenFunction * fFunc1Dim
Definition FitUtil.h:263
ROOT::Math::IntegratorMultiDim * fIgNDim
Definition FitUtil.h:262
void SetFunction(const ParamFunc &func, const double *p=nullptr, ROOT::Math::IntegrationOneDim::Type igType=ROOT::Math::IntegrationOneDim::kDEFAULT)
Definition FitUtil.h:140
IntegralEvaluator(const IntegralEvaluator &rhs)=delete
IntegralEvaluator(const ParamFunc &func, const double *p, bool useIntegral=true, ROOT::Math::IntegrationOneDim::Type igType=ROOT::Math::IntegrationOneDim::kDEFAULT)
Definition FitUtil.h:131
IntegralEvaluator & operator=(const IntegralEvaluator &rhs)=delete
double F1(double x) const
Definition FitUtil.h:190
ROOT::Math::IntegratorOneDim * fIg1Dim
Definition FitUtil.h:261
double FN(const double *x) const
Definition FitUtil.h:196
ROOT::Math::IMultiGenFunction * fFuncNDim
Definition FitUtil.h:264
void SetParameters(const double *p)
Definition FitUtil.h:170
double ExecFunc(T *f, const double *x, const double *p) const
Definition FitUtil.h:224
double operator()(const double *x1, const double *x2)
Definition FitUtil.h:204
LikelihoodAux(double logv=0.0, double w=0.0, double w2=0.0)
Definition FitUtil.h:100
LikelihoodAux & operator+=(const LikelihoodAux &l)
Definition FitUtil.h:107
LikelihoodAux operator+(const LikelihoodAux &l) const
Definition FitUtil.h:102
LikelihoodAux operator+(const LikelihoodAux &l) const
Definition FitUtil.h:79
LikelihoodAux(T logv={}, T w={}, T w2={})
Definition FitUtil.h:77
LikelihoodAux & operator+=(const LikelihoodAux &l)
Definition FitUtil.h:84
Class describing the un-binned data sets (just x coordinates values) of any dimensions.
Definition UnBinData.h:46
Documentation for the abstract class IBaseFunctionMultiDim.
Definition IFunction.h:61
Interface (abstract class) for generic functions objects of one-dimension Provides a method to evalua...
Definition IFunction.h:112
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 ...
User class for performing multidimensional integration.
double Integral(const double *xmin, const double *xmax)
evaluate the integral with the previously given function between xmin[] and xmax[]
void SetFunction(Function &f, unsigned int dim)
set integration function using a generic function implementing the operator()(double *x) The dimensio...
User Class for performing numerical integration of a function in one dimension.
Definition Integrator.h:98
void SetFunction(Function &f)
method to set the a generic integration function
Definition Integrator.h:492
double Integral(Function &f, double a, double b)
evaluate the Integral of a function f over the defined interval (a,b)
Definition Integrator.h:499
Template class to wrap any member function of a class taking a double and returning a double in a 1D ...
auto MapReduce(F func, unsigned nTimes, R redfunc) -> InvokeResult_t< F >
Execute a function without arguments several times (Map) and accumulate the results into a single val...
A pseudo container class which is a generator of indices.
Definition TSeq.hxx:67
This class provides a simple interface to execute the same task multiple times in parallel threads,...
auto MapReduce(F func, unsigned nTimes, R redfunc) -> InvokeResult_t< F >
Execute a function nTimes in parallel (Map) and accumulate the results into a single value (Reduce).
Type
enumeration specifying the integration types.
@ 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
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
ROOT::Math::IParamMultiFunction IModelFunction
Definition FitUtil.h:64
double EvaluatePoissonBinPdf(const IModelFunction &func, const BinData &data, const double *x, unsigned int ipoint, double *g=nullptr, double *h=nullptr, bool hasGrad=false, bool fullHessian=false)
evaluate the pdf contribution to the Poisson LogL given a model function and the BinPoint data.
Definition FitUtil.cxx:1297
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 EvaluateChi2Residual(const IModelFunction &func, const BinData &data, const double *p, unsigned int ipoint, double *g=nullptr, double *h=nullptr, bool hasGrad=false, bool fullHessian=false)
evaluate the residual contribution to the Chi2 given a model function and the BinPoint data and if th...
Definition FitUtil.cxx:545
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:424
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 EvaluatePdf(const IModelFunction &func, const UnBinData &data, const double *p, unsigned int ipoint, double *g=nullptr, double *h=nullptr, bool hasGrad=false, bool fullHessian=false)
evaluate the pdf contribution to the LogL given a model function and the BinPoint data.
Definition FitUtil.cxx:891
unsigned setAutomaticChunking(unsigned nEvents)
Definition FitUtil.cxx:1841
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
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
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
Double_t Double_v
Definition Types.h:55
DataOptions : simple structure holding the options on how the data are filled.
Definition DataOptions.h:28
bool fIntegral
use integral of bin content instead of bin center (default is false)
Definition DataOptions.h:47
bool fExpErrors
use expected errors from the function and not from the data
Definition DataOptions.h:53
bool fBinVolume
normalize data by the bin volume (it is used in the Poisson likelihood fits)
Definition DataOptions.h:48
static double EvalLogL(const IModelFunctionTempl< double > &func, const UnBinData &data, const double *p, int iWeight, bool extended, unsigned int &nPoints, ::ROOT::EExecutionPolicy executionPolicy, unsigned nChunks=0)
Definition FitUtil.h:1413
static double EvalChi2Effective(const IModelFunctionTempl< double > &func, const BinData &data, const double *p, unsigned int &nPoints)
Definition FitUtil.h:1427
static void EvalLogLGradient(const IModelFunctionTempl< double > &func, const UnBinData &data, const double *p, double *g, unsigned int &nPoints, ::ROOT::EExecutionPolicy executionPolicy=::ROOT::EExecutionPolicy::kSequential, unsigned nChunks=0)
Definition FitUtil.h:1464
static void EvalChi2Gradient(const IModelFunctionTempl< double > &func, const BinData &data, const double *p, double *g, unsigned int &nPoints, ::ROOT::EExecutionPolicy executionPolicy=::ROOT::EExecutionPolicy::kSequential, unsigned nChunks=0)
Definition FitUtil.h:1431
static double EvalPdf(const IModelFunctionTempl< double > &func, const UnBinData &data, const double *p, unsigned int i, double *g, double *h, bool hasGrad, bool fullHessian)
Definition FitUtil.h:1451
static double EvalChi2(const IModelFunction &func, const BinData &data, const double *p, unsigned int &nPoints, ::ROOT::EExecutionPolicy executionPolicy, unsigned nChunks=0)
Definition FitUtil.h:1399
static double EvalPoissonBinPdf(const IModelFunctionTempl< double > &func, const BinData &data, const double *p, unsigned int i, double *g, double *h, bool hasGrad, bool fullHessian)
evaluate the pdf (Poisson) contribution to the logl (return actually log of pdf) and its gradient
Definition FitUtil.h:1447
static void EvalPoissonLogLGradient(const IModelFunctionTempl< double > &func, const BinData &data, const double *p, double *g, unsigned int &nPoints, ::ROOT::EExecutionPolicy executionPolicy=::ROOT::EExecutionPolicy::kSequential, unsigned nChunks=0)
Definition FitUtil.h:1456
static double EvalPoissonLogL(const IModelFunctionTempl< double > &func, const BinData &data, const double *p, int iWeight, bool extended, unsigned int &nPoints, ::ROOT::EExecutionPolicy executionPolicy, unsigned nChunks=0)
Definition FitUtil.h:1420
static double EvalChi2Residual(const IModelFunctionTempl< double > &func, const BinData &data, const double *p, unsigned int i, double *g, double *h, bool hasGrad, bool fullHessian)
Definition FitUtil.h:1439
TMarker m
Definition textangle.C:8
TLine l
Definition textangle.C:4