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