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