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