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