Logo ROOT   6.14/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 
16 #include "Math/IParamFunctionfwd.h"
17 #include "Math/IParamFunction.h"
18 
19 #ifdef R__USE_IMT
20 #include "ROOT/TThreadExecutor.hxx"
21 #endif
23 
24 #include "Fit/BinData.h"
25 #include "Fit/UnBinData.h"
26 #include "Fit/FitExecutionPolicy.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
40 namespace vecCore {
41 template <class T>
42 vecCore::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 
52 namespace 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 */
62 namespace FitUtil {
63 
66 
67  template <class T>
69 
70  template <class T>
72 
73  // internal class defining
74  template <class T>
75  class LikelihoodAux {
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);
155  fIg1Dim = new ROOT::Math::IntegratorOneDim();
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 =
160  new ROOT::Math::WrappedMemMultiFunction<IntegralEvaluator, double (IntegralEvaluator::*)(const double *)
161  const>(*this, &IntegralEvaluator::FN, fDim);
162  fIgNDim = new ROOT::Math::IntegratorMultiDim();
163  fIgNDim->SetFunction(*fFuncNDim);
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) {
231  ROOT::Double_v xx;
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
249  IntegralEvaluator &operator=(const IntegralEvaluator &rhs);
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 
385  nPoints = 0; // count the effective non-zero points
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  nPoints++;
434 
435  T tmp = (y - fval) * invErrorVec;
436  T chi2 = tmp * tmp;
437 
438 
439  // avoid inifinity or nan in chi2 values due to wrong function values
440  auto m = vecCore::Mask_v<T>(chi2 > maxResValue);
441 
442  vecCore::MaskedAssign<T>(chi2, m, maxResValue);
443 
444  return chi2;
445  };
446 
447  auto redFunction = [](const std::vector<T> &objs) {
448  return std::accumulate(objs.begin(), objs.end(), T{});
449  };
450 
451 #ifndef R__USE_IMT
452  (void)nChunks;
453 
454  // If IMT is disabled, force the execution policy to the serial case
455  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
456  Warning("FitUtil::EvaluateChi2", "Multithread execution policy requires IMT, which is disabled. Changing "
457  "to ROOT::Fit::ExecutionPolicy::kSerial.");
458  executionPolicy = ROOT::Fit::ExecutionPolicy::kSerial;
459  }
460 #endif
461 
462  T res{};
463  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kSerial) {
465  res = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, data.Size()/vecSize), redFunction);
466 #ifdef R__USE_IMT
467  } else if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
468  auto chunks = nChunks != 0 ? nChunks : setAutomaticChunking(data.Size() / vecSize);
470  res = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, data.Size() / vecSize), redFunction, chunks);
471 #endif
472  } else {
473  Error("FitUtil::EvaluateChi2", "Execution policy unknown. Avalaible choices:\n ROOT::Fit::ExecutionPolicy::kSerial (default)\n ROOT::Fit::ExecutionPolicy::kMultithread (requires IMT)\n");
474  }
475  nPoints = n;
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::Fit::ExecutionPolicy executionPolicy, unsigned nChunks = 0)
488  {
489  // evaluate the LogLikelihood
490  unsigned int n = data.Size();
491 
492  //unsigned int nRejected = 0;
493  bool normalizeFunc = false;
494 
495  // set parameters of the function to cache integral value
496 #ifdef USE_PARAMCACHE
497  (const_cast<IModelFunctionTempl<T> &>(func)).SetParameters(p);
498 #endif
499 
500 #ifdef R__USE_IMT
501  // in case parameter needs to be propagated to user function use trick to set parameters by calling one time the function
502  // this will be done in sequential mode and parameters can be set in a thread safe manner
503  if (!normalizeFunc) {
504  if (data.NDim() == 1) {
505  T x;
506  vecCore::Load<T>(x, data.GetCoordComponent(0, 0));
507  func( &x, p);
508  }
509  else {
510  std::vector<T> x(data.NDim());
511  for (unsigned int j = 0; j < data.NDim(); ++j)
512  vecCore::Load<T>(x[j], data.GetCoordComponent(0, j));
513  func( x.data(), p);
514  }
515  }
516 #endif
517 
518  // this is needed if function must be normalized
519  double norm = 1.0;
520  if (normalizeFunc) {
521  // compute integral of the function
522  std::vector<double> xmin(data.NDim());
523  std::vector<double> xmax(data.NDim());
524  IntegralEvaluator<IModelFunctionTempl<T>> igEval(func, p, true);
525  // compute integral in the ranges where is defined
526  if (data.Range().Size() > 0) {
527  norm = 0;
528  for (unsigned int ir = 0; ir < data.Range().Size(); ++ir) {
529  data.Range().GetRange(&xmin[0], &xmax[0], ir);
530  norm += igEval.Integral(xmin.data(), xmax.data());
531  }
532  } else {
533  // use (-inf +inf)
534  data.Range().GetRange(&xmin[0], &xmax[0]);
535  // check if funcition is zero at +- inf
536  T xmin_v, xmax_v;
537  vecCore::Load<T>(xmin_v, xmin.data());
538  vecCore::Load<T>(xmax_v, xmax.data());
539  if (vecCore::ReduceAdd(func(&xmin_v, p)) != 0 || vecCore::ReduceAdd(func(&xmax_v, p)) != 0) {
540  MATH_ERROR_MSG("FitUtil::EvaluateLogLikelihood", "A range has not been set and the function is not zero at +/- inf");
541  return 0;
542  }
543  norm = igEval.Integral(&xmin[0], &xmax[0]);
544  }
545  }
546 
547  // needed to compute effective global weight in case of extended likelihood
548 
549  auto vecSize = vecCore::VectorSize<T>();
550  unsigned int numVectors = n / vecSize;
551 
552  auto mapFunction = [ &, p](const unsigned i) {
553  T W{};
554  T W2{};
555  T fval{};
556 
557  (void)p; /* avoid unused lambda capture warning if PARAMCACHE is disabled */
558 
559  T x1;
560  vecCore::Load<T>(x1, data.GetCoordComponent(i * vecSize, 0));
561  const T *x = nullptr;
562  unsigned int ndim = data.NDim();
563  std::vector<T> xc;
564  if (ndim > 1) {
565  xc.resize(ndim);
566  xc[0] = x1;
567  for (unsigned int j = 1; j < ndim; ++j)
568  vecCore::Load<T>(xc[j], data.GetCoordComponent(i * vecSize, j));
569  x = xc.data();
570  } else {
571  x = &x1;
572  }
573 
574 #ifdef USE_PARAMCACHE
575  fval = func(x);
576 #else
577  fval = func(x, p);
578 #endif
579 
580 #ifdef DEBUG_FITUTIL
581  if (i < 5 || (i > numVectors-5) ) {
582  if (ndim == 1) std::cout << i << " x " << x[0] << " fval = " << fval;
583  else std::cout << i << " x " << x[0] << " y " << x[1] << " fval = " << fval;
584  }
585 #endif
586 
587  if (normalizeFunc) fval = fval * (1 / norm);
588 
589  // function EvalLog protects against negative or too small values of fval
590  auto logval = ROOT::Math::Util::EvalLog(fval);
591  if (iWeight > 0) {
592  T weight{};
593  if (data.WeightsPtr(i) == nullptr)
594  weight = 1;
595  else
596  vecCore::Load<T>(weight, data.WeightsPtr(i*vecSize));
597  logval *= weight;
598  if (iWeight == 2) {
599  logval *= weight; // use square of weights in likelihood
600  if (!extended) {
601  // needed sum of weights and sum of weight square if likelkihood is extended
602  W = weight;
603  W2 = weight * weight;
604  }
605  }
606  }
607 #ifdef DEBUG_FITUTIL
608  if (i < 5 || (i > numVectors-5) ) {
609  std::cout << " " << fval << " logfval " << logval << std::endl;
610  }
611 #endif
612 
613  nPoints++;
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::Fit::ExecutionPolicy::kMultithread) {
629  Warning("FitUtil::EvaluateLogL", "Multithread execution policy requires IMT, which is disabled. Changing "
630  "to ROOT::Fit::ExecutionPolicy::kSerial.");
631  executionPolicy = ROOT::Fit::ExecutionPolicy::kSerial;
632  }
633 #endif
634 
635  T logl_v{};
636  T sumW_v{};
637  T sumW2_v{};
639  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kSerial) {
641  resArray = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, data.Size() / vecSize), redFunction);
642 #ifdef R__USE_IMT
643  } else if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
644  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::Fit::ExecutionPolicy::kSerial (default)\n ROOT::Fit::ExecutionPolicy::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 = 0.;
670  double sumW = 0.;
671  double sumW2 = 0;;
672 
673  for (unsigned vIt = 0; vIt < vecSize; vIt++) {
674  logl += logl_v[vIt];
675  sumW += sumW_v[vIt];
676  sumW2 += sumW2_v[vIt];
677  }
678 
679  if (extended) {
680  // add Poisson extended term
681  double extendedTerm = 0; // extended term in likelihood
682  double nuTot = 0;
683  // nuTot is integral of function in the range
684  // if function has been normalized integral has been already computed
685  if (!normalizeFunc) {
686  IntegralEvaluator<IModelFunctionTempl<T>> igEval(func, p, true);
687  std::vector<double> xmin(data.NDim());
688  std::vector<double> xmax(data.NDim());
689 
690  // compute integral in the ranges where is defined
691  if (data.Range().Size() > 0) {
692  nuTot = 0;
693  for (unsigned int ir = 0; ir < data.Range().Size(); ++ir) {
694  data.Range().GetRange(&xmin[0], &xmax[0], ir);
695  nuTot += igEval.Integral(xmin.data(), xmax.data());
696  }
697  } else {
698  // use (-inf +inf)
699  data.Range().GetRange(&xmin[0], &xmax[0]);
700  // check if funcition is zero at +- inf
701  T xmin_v, xmax_v;
702  vecCore::Load<T>(xmin_v, xmin.data());
703  vecCore::Load<T>(xmax_v, xmax.data());
704  if (vecCore::ReduceAdd(func(&xmin_v, p)) != 0 || vecCore::ReduceAdd(func(&xmax_v, p)) != 0) {
705  MATH_ERROR_MSG("FitUtil::EvaluateLogLikelihood", "A range has not been set and the function is not zero at +/- inf");
706  return 0;
707  }
708  nuTot = igEval.Integral(&xmin[0], &xmax[0]);
709  }
710 
711  // force to be last parameter value
712  //nutot = p[func.NDim()-1];
713  if (iWeight != 2)
714  extendedTerm = - nuTot; // no need to add in this case n log(nu) since is already computed before
715  else {
716  // case use weight square in likelihood : compute total effective weight = sw2/sw
717  // ignore for the moment case when sumW is zero
718  extendedTerm = - (sumW2 / sumW) * nuTot;
719  }
720 
721  } else {
722  nuTot = norm;
723  extendedTerm = - nuTot + double(n) * ROOT::Math::Util::EvalLog(nuTot);
724  // in case of weights need to use here sum of weights (to be done)
725  }
726 
727  logl += extendedTerm;
728  }
729 
730 #ifdef DEBUG_FITUTIL
731  std::cout << "Evaluated log L for parameters (";
732  for (unsigned int ip = 0; ip < func.NPar(); ++ip)
733  std::cout << " " << p[ip];
734  std::cout << ") nll = " << -logl << std::endl;
735 #endif
736 
737  // reset the number of fitting data points
738  // nPoints = n;
739  // std::cout<<", n: "<<nPoints<<std::endl;
740  nPoints = 0;
741  return -logl;
742 
743  }
744 
745  static double EvalPoissonLogL(const IModelFunctionTempl<T> &func, const BinData &data, const double *p,
746  int iWeight, bool extended, unsigned int,
747  ROOT::Fit::ExecutionPolicy executionPolicy, unsigned nChunks = 0)
748  {
749  // evaluate the Poisson Log Likelihood
750  // for binned likelihood fits
751  // this is Sum ( f(x_i) - y_i * log( f (x_i) ) )
752  // add as well constant term for saturated model to make it like a Chi2/2
753  // by default is etended. If extended is false the fit is not extended and
754  // the global poisson term is removed (i.e is a binomial fit)
755  // (remember that in this case one needs to have a function with a fixed normalization
756  // like in a non extended binned fit)
757  //
758  // if use Weight use a weighted dataset
759  // iWeight = 1 ==> logL = Sum( w f(x_i) )
760  // case of iWeight==1 is actually identical to weight==0
761  // iWeight = 2 ==> logL = Sum( w*w * f(x_i) )
762  //
763 
764 #ifdef USE_PARAMCACHE
765  (const_cast<IModelFunctionTempl<T> &>(func)).SetParameters(p);
766 #endif
767  auto vecSize = vecCore::VectorSize<T>();
768  // get fit option and check case of using integral of bins
769  const DataOptions &fitOpt = data.Opt();
770  if (fitOpt.fExpErrors || fitOpt.fIntegral)
771  Error("FitUtil::EvaluateChi2",
772  "The vectorized implementation doesn't support Integrals or BinVolume\n. Aborting operation.");
773  bool useW2 = (iWeight == 2);
774 
775  auto mapFunction = [&](unsigned int i) {
776  T y;
777  vecCore::Load<T>(y, data.ValuePtr(i * vecSize));
778  T fval{};
779 
780  if (data.NDim() > 1) {
781  std::vector<T> x(data.NDim());
782  for (unsigned int j = 0; j < data.NDim(); ++j)
783  vecCore::Load<T>(x[j], data.GetCoordComponent(i * vecSize, j));
784 #ifdef USE_PARAMCACHE
785  fval = func(x.data());
786 #else
787  fval = func(x.data(), p);
788 #endif
789  // one -dim case
790  } else {
791  T x;
792  vecCore::Load<T>(x, data.GetCoordComponent(i * vecSize, 0));
793 #ifdef USE_PARAMCACHE
794  fval = func(&x);
795 #else
796  fval = func(&x, p);
797 #endif
798  }
799 
800  // EvalLog protects against 0 values of fval but don't want to add in the -log sum
801  // negative values of fval
802  vecCore::MaskedAssign<T>(fval, fval < 0.0, 0.0);
803 
804  T nloglike{}; // negative loglikelihood
805 
806  if (useW2) {
807  // apply weight correction . Effective weight is error^2/ y
808  // and expected events in bins is fval/weight
809  // can apply correction only when y is not zero otherwise weight is undefined
810  // (in case of weighted likelihood I don't care about the constant term due to
811  // the saturated model)
812  assert (data.GetErrorType() != ROOT::Fit::BinData::ErrorType::kNoError);
813  T error = 0.0;
814  vecCore::Load<T>(error, data.ErrorPtr(i * vecSize));
815  // for empty bin use the average weight computed from the total data weight
816  auto m = vecCore::Mask_v<T>(y != 0.0);
817  auto weight = vecCore::Blend(m,(error * error) / y, T(data.SumOfError2()/ data.SumOfContent()) );
818  if (extended) {
819  nloglike = weight * ( fval - y);
820  }
821  vecCore::MaskedAssign<T>(nloglike, y != 0, nloglike + weight * y *( ROOT::Math::Util::EvalLog(y) - ROOT::Math::Util::EvalLog(fval)) );
822 
823  } else {
824  // standard case no weights or iWeight=1
825  // this is needed for Poisson likelihood (which are extened and not for multinomial)
826  // the formula below include constant term due to likelihood of saturated model (f(x) = y)
827  // (same formula as in Baker-Cousins paper, page 439 except a factor of 2
828  if (extended) nloglike = fval - y;
829 
830  vecCore::MaskedAssign<T>(
831  nloglike, y > 0, nloglike + y * (ROOT::Math::Util::EvalLog(y) - ROOT::Math::Util::EvalLog(fval)));
832  }
833 
834  return nloglike;
835  };
836 
837 #ifdef R__USE_IMT
838  auto redFunction = [](const std::vector<T> &objs) { return std::accumulate(objs.begin(), objs.end(), T{}); };
839 #else
840  (void)nChunks;
841 
842  // If IMT is disabled, force the execution policy to the serial case
843  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
844  Warning("FitUtil::Evaluate<T>::EvalPoissonLogL",
845  "Multithread execution policy requires IMT, which is disabled. Changing "
846  "to ROOT::Fit::ExecutionPolicy::kSerial.");
847  executionPolicy = ROOT::Fit::ExecutionPolicy::kSerial;
848  }
849 #endif
850 
851  T res{};
852  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kSerial) {
853  for (unsigned int i = 0; i < (data.Size() / vecSize); i++) {
854  res += mapFunction(i);
855  }
856 #ifdef R__USE_IMT
857  } else if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
858  auto chunks = nChunks != 0 ? nChunks : setAutomaticChunking(data.Size() / vecSize);
860  res = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, data.Size() / vecSize), redFunction, chunks);
861 #endif
862  } else {
863  Error(
864  "FitUtil::Evaluate<T>::EvalPoissonLogL",
865  "Execution policy unknown. Avalaible choices:\n ROOT::Fit::ExecutionPolicy::kSerial (default)\n ROOT::Fit::ExecutionPolicy::kMultithread (requires IMT)\n");
866  }
867 
868  // Last padded SIMD vector of elements
869  if (data.Size() % vecSize != 0)
870  vecCore::MaskedAssign(res, vecCore::Int2Mask<T>(data.Size() % vecSize),
871  res + mapFunction(data.Size() / vecSize));
872 
873  return vecCore::ReduceAdd(res);
874  }
875 
876  static double EvalChi2Effective(const IModelFunctionTempl<T> &, const BinData &, const double *, unsigned int &)
877  {
878  Error("FitUtil::Evaluate<T>::EvalChi2Effective", "The vectorized evaluation of the Chi2 with coordinate errors is still not supported");
879  return -1.;
880  }
881 
882  // Compute a mask to filter out infinite numbers and NaN values.
883  // The argument rval is updated so infinite numbers and NaN values are replaced by
884  // maximum finite values (preserving the original sign).
885  static vecCore::Mask<T> CheckInfNaNValues(T &rval)
886  {
888 
889  // Case +inf or nan
890  vecCore::MaskedAssign(rval, !mask, +vecCore::NumericLimits<T>::Max());
891 
892  // Case -inf
893  vecCore::MaskedAssign(rval, !mask && rval < 0, -vecCore::NumericLimits<T>::Max());
894 
895  return mask;
896  }
897 
898  static void EvalChi2Gradient(const IModelFunctionTempl<T> &f, const BinData &data, const double *p, double *grad,
899  unsigned int &nPoints,
901  unsigned nChunks = 0)
902  {
903  // evaluate the gradient of the chi2 function
904  // this function is used when the model function knows how to calculate the derivative and we can
905  // avoid that the minimizer re-computes them
906  //
907  // case of chi2 effective (errors on coordinate) is not supported
908 
909  if (data.HaveCoordErrors()) {
910  MATH_ERROR_MSG("FitUtil::EvaluateChi2Gradient",
911  "Error on the coordinates are not used in calculating Chi2 gradient");
912  return; // it will assert otherwise later in GetPoint
913  }
914 
915  const IGradModelFunctionTempl<T> *fg = dynamic_cast<const IGradModelFunctionTempl<T> *>(&f);
916  assert(fg != nullptr); // must be called by a gradient function
917 
918  const IGradModelFunctionTempl<T> &func = *fg;
919 
920  const DataOptions &fitOpt = data.Opt();
921  if (fitOpt.fBinVolume || fitOpt.fIntegral || fitOpt.fExpErrors)
922  Error("FitUtil::EvaluateChi2Gradient", "The vectorized implementation doesn't support Integrals,"
923  "BinVolume or ExpErrors\n. Aborting operation.");
924 
925  unsigned int npar = func.NPar();
926  auto vecSize = vecCore::VectorSize<T>();
927  unsigned initialNPoints = data.Size();
928  unsigned numVectors = initialNPoints / vecSize;
929 
930  // numVectors + 1 because of the padded data (call to mapFunction with i = numVectors after the main loop)
931  std::vector<vecCore::Mask<T>> validPointsMasks(numVectors + 1);
932 
933  auto mapFunction = [&](const unsigned int i) {
934  // set all vector values to zero
935  std::vector<T> gradFunc(npar);
936  std::vector<T> pointContributionVec(npar);
937 
938  T x1, y, invError;
939 
940  vecCore::Load<T>(x1, data.GetCoordComponent(i * vecSize, 0));
941  vecCore::Load<T>(y, data.ValuePtr(i * vecSize));
942  const auto invErrorPtr = data.ErrorPtr(i * vecSize);
943 
944  if (invErrorPtr == nullptr)
945  invError = 1;
946  else
947  vecCore::Load<T>(invError, invErrorPtr);
948 
949  // TODO: Check error options and invert if needed
950 
951  T fval = 0;
952 
953  const T *x = nullptr;
954 
955  unsigned int ndim = data.NDim();
956  // need to declare vector outside if statement
957  // otherwise pointer will be invalid
958  std::vector<T> xc;
959  if (ndim > 1) {
960  xc.resize(ndim);
961  xc[0] = x1;
962  for (unsigned int j = 1; j < ndim; ++j)
963  vecCore::Load<T>(xc[j], data.GetCoordComponent(i * vecSize, j));
964  x = xc.data();
965  } else {
966  x = &x1;
967  }
968 
969  fval = func(x, p);
970  func.ParameterGradient(x, p, &gradFunc[0]);
971 
972  validPointsMasks[i] = CheckInfNaNValues(fval);
973  if (vecCore::MaskEmpty(validPointsMasks[i])) {
974  // Return a zero contribution to all partial derivatives on behalf of the current points
975  return pointContributionVec;
976  }
977 
978  // loop on the parameters
979  for (unsigned int ipar = 0; ipar < npar; ++ipar) {
980  // avoid singularity in the function (infinity and nan ) in the chi2 sum
981  // eventually add possibility of excluding some points (like singularity)
982  validPointsMasks[i] = CheckInfNaNValues(gradFunc[ipar]);
983 
984  if (vecCore::MaskEmpty(validPointsMasks[i])) {
985  break; // exit loop on parameters
986  }
987 
988  // calculate derivative point contribution (only for valid points)
989  vecCore::MaskedAssign(pointContributionVec[ipar], validPointsMasks[i],
990  -2.0 * (y - fval) * invError * invError * gradFunc[ipar]);
991  }
992 
993  return pointContributionVec;
994  };
995 
996  // Reduce the set of vectors by summing its equally-indexed components
997  auto redFunction = [&](const std::vector<std::vector<T>> &partialResults) {
998  std::vector<T> result(npar);
999 
1000  for (auto const &pointContributionVec : partialResults) {
1001  for (unsigned int parameterIndex = 0; parameterIndex < npar; parameterIndex++)
1002  result[parameterIndex] += pointContributionVec[parameterIndex];
1003  }
1004 
1005  return result;
1006  };
1007 
1008  std::vector<T> gVec(npar);
1009  std::vector<double> g(npar);
1010 
1011 #ifndef R__USE_IMT
1012  // to fix compiler warning
1013  (void)nChunks;
1014 
1015  // If IMT is disabled, force the execution policy to the serial case
1016  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
1017  Warning("FitUtil::EvaluateChi2Gradient",
1018  "Multithread execution policy requires IMT, which is disabled. Changing "
1019  "to ROOT::Fit::ExecutionPolicy::kSerial.");
1020  executionPolicy = ROOT::Fit::ExecutionPolicy::kSerial;
1021  }
1022 #endif
1023 
1024  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kSerial) {
1026  gVec = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, numVectors), redFunction);
1027  }
1028 #ifdef R__USE_IMT
1029  else if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
1030  auto chunks = nChunks != 0 ? nChunks : setAutomaticChunking(numVectors);
1031  ROOT::TThreadExecutor pool;
1032  gVec = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, numVectors), redFunction, chunks);
1033  }
1034 #endif
1035  else {
1036  Error(
1037  "FitUtil::EvaluateChi2Gradient",
1038  "Execution policy unknown. Avalaible choices:\n 0: Serial (default)\n 1: MultiThread (requires IMT)\n");
1039  }
1040 
1041  // Compute the contribution from the remaining points
1042  unsigned int remainingPoints = initialNPoints % vecSize;
1043  if (remainingPoints > 0) {
1044  auto remainingPointsContribution = mapFunction(numVectors);
1045  // Add the contribution from the valid remaining points and store the result in the output variable
1046  auto remainingMask = vecCore::Int2Mask<T>(remainingPoints);
1047  for (unsigned int param = 0; param < npar; param++) {
1048  vecCore::MaskedAssign(gVec[param], remainingMask, gVec[param] + remainingPointsContribution[param]);
1049  }
1050  }
1051  // reduce final gradient result from T to double
1052  for (unsigned int param = 0; param < npar; param++) {
1053  grad[param] = vecCore::ReduceAdd(gVec[param]);
1054  }
1055 
1056  // correct the number of points
1057  nPoints = initialNPoints;
1058 
1059  if (std::any_of(validPointsMasks.begin(), validPointsMasks.end(),
1060  [](vecCore::Mask<T> validPoints) { return !vecCore::MaskFull(validPoints); })) {
1061  unsigned nRejected = 0;
1062 
1063  for (const auto &mask : validPointsMasks) {
1064  for (unsigned int i = 0; i < vecSize; i++) {
1065  nRejected += !vecCore::Get(mask, i);
1066  }
1067  }
1068 
1069  assert(nRejected <= initialNPoints);
1070  nPoints = initialNPoints - nRejected;
1071 
1072  if (nPoints < npar) {
1073  MATH_ERROR_MSG("FitUtil::EvaluateChi2Gradient",
1074  "Too many points rejected for overflow in gradient calculation");
1075  }
1076  }
1077  }
1078 
1079  static double EvalChi2Residual(const IModelFunctionTempl<T> &, const BinData &, const double *, unsigned int, double *)
1080  {
1081  Error("FitUtil::Evaluate<T>::EvalChi2Residual", "The vectorized evaluation of the Chi2 with the ith residual is still not supported");
1082  return -1.;
1083  }
1084 
1085  /// evaluate the pdf (Poisson) contribution to the logl (return actually log of pdf)
1086  /// and its gradient
1087  static double EvalPoissonBinPdf(const IModelFunctionTempl<T> &, const BinData &, const double *, unsigned int , double * ) {
1088  Error("FitUtil::Evaluate<T>::EvaluatePoissonBinPdf", "The vectorized evaluation of the BinnedLikelihood fit evaluated point by point is still not supported");
1089  return -1.;
1090  }
1091 
1092  static void
1093  EvalPoissonLogLGradient(const IModelFunctionTempl<T> &f, const BinData &data, const double *p, double *grad,
1094  unsigned int &,
1096  unsigned nChunks = 0)
1097  {
1098  // evaluate the gradient of the Poisson log likelihood function
1099 
1100  const IGradModelFunctionTempl<T> *fg = dynamic_cast<const IGradModelFunctionTempl<T> *>(&f);
1101  assert(fg != nullptr); // must be called by a grad function
1102 
1103  const IGradModelFunctionTempl<T> &func = *fg;
1104 
1105  (const_cast<IGradModelFunctionTempl<T> &>(func)).SetParameters(p);
1106 
1107 
1108  const DataOptions &fitOpt = data.Opt();
1109  if (fitOpt.fBinVolume || fitOpt.fIntegral || fitOpt.fExpErrors)
1110  Error("FitUtil::EvaluatePoissonLogLGradient", "The vectorized implementation doesn't support Integrals,"
1111  "BinVolume or ExpErrors\n. Aborting operation.");
1112 
1113  unsigned int npar = func.NPar();
1114  auto vecSize = vecCore::VectorSize<T>();
1115  unsigned initialNPoints = data.Size();
1116  unsigned numVectors = initialNPoints / vecSize;
1117 
1118  auto mapFunction = [&](const unsigned int i) {
1119  // set all vector values to zero
1120  std::vector<T> gradFunc(npar);
1121  std::vector<T> pointContributionVec(npar);
1122 
1123  T x1, y;
1124 
1125  vecCore::Load<T>(x1, data.GetCoordComponent(i * vecSize, 0));
1126  vecCore::Load<T>(y, data.ValuePtr(i * vecSize));
1127 
1128  T fval = 0;
1129 
1130  const T *x = nullptr;
1131 
1132  unsigned ndim = data.NDim();
1133  std::vector<T> xc;
1134  if (ndim > 1) {
1135  xc.resize(ndim);
1136  xc[0] = x1;
1137  for (unsigned int j = 1; j < ndim; ++j)
1138  vecCore::Load<T>(xc[j], data.GetCoordComponent(i * vecSize, j));
1139  x = xc.data();
1140  } else {
1141  x = &x1;
1142  }
1143 
1144  fval = func(x, p);
1145  func.ParameterGradient(x, p, &gradFunc[0]);
1146 
1147  // correct the gradient
1148  for (unsigned int ipar = 0; ipar < npar; ++ipar) {
1149  vecCore::Mask<T> positiveValuesMask = fval > 0;
1150 
1151  // df/dp * (1. - y/f )
1152  vecCore::MaskedAssign(pointContributionVec[ipar], positiveValuesMask, gradFunc[ipar] * (1. - y / fval));
1153 
1154  vecCore::Mask<T> validNegativeValuesMask = !positiveValuesMask && gradFunc[ipar] != 0;
1155 
1156  if (!vecCore::MaskEmpty(validNegativeValuesMask)) {
1158  const T kdmax2 = vecCore::NumericLimits<T>::Max() / (4 * initialNPoints);
1159  T gg = kdmax1 * gradFunc[ipar];
1160  pointContributionVec[ipar] = -vecCore::Blend(gg > 0, vecCore::math::Min(gg, kdmax2), vecCore::math::Max(gg, -kdmax2));
1161  }
1162  }
1163 
1164 #ifdef DEBUG_FITUTIL
1165  {
1167  if (i < 5 || (i > data.Size()-5) ) {
1168  if (data.NDim() > 1) std::cout << i << " x " << x[0] << " y " << x[1];
1169  else std::cout << i << " x " << x[0];
1170  std::cout << " func " << fval << " gradient ";
1171  for (unsigned int ii = 0; ii < npar; ++ii) std::cout << " " << pointContributionVec[ii];
1172  std::cout << "\n";
1173  }
1174  }
1175 #endif
1176 
1177  return pointContributionVec;
1178  };
1179 
1180  // Vertically reduce the set of vectors by summing its equally-indexed components
1181  auto redFunction = [&](const std::vector<std::vector<T>> &partialResults) {
1182  std::vector<T> result(npar);
1183 
1184  for (auto const &pointContributionVec : partialResults) {
1185  for (unsigned int parameterIndex = 0; parameterIndex < npar; parameterIndex++)
1186  result[parameterIndex] += pointContributionVec[parameterIndex];
1187  }
1188 
1189  return result;
1190  };
1191 
1192  std::vector<T> gVec(npar);
1193 
1194 #ifndef R__USE_IMT
1195  // to fix compiler warning
1196  (void)nChunks;
1197 
1198  // If IMT is disabled, force the execution policy to the serial case
1199  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
1200  Warning("FitUtil::EvaluatePoissonLogLGradient",
1201  "Multithread execution policy requires IMT, which is disabled. Changing "
1202  "to ROOT::Fit::ExecutionPolicy::kSerial.");
1203  executionPolicy = ROOT::Fit::ExecutionPolicy::kSerial;
1204  }
1205 #endif
1206 
1207  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kSerial) {
1209  gVec = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, numVectors), redFunction);
1210  }
1211 #ifdef R__USE_IMT
1212  else if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
1213  auto chunks = nChunks != 0 ? nChunks : setAutomaticChunking(numVectors);
1214  ROOT::TThreadExecutor pool;
1215  gVec = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, numVectors), redFunction, chunks);
1216  }
1217 #endif
1218  else {
1219  Error("FitUtil::EvaluatePoissonLogLGradient", "Execution policy unknown. Avalaible choices:\n "
1220  "ROOT::Fit::ExecutionPolicy::kSerial (default)\n "
1221  "ROOT::Fit::ExecutionPolicy::kMultithread (requires IMT)\n");
1222  }
1223 
1224 
1225  // Compute the contribution from the remaining points
1226  unsigned int remainingPoints = initialNPoints % vecSize;
1227  if (remainingPoints > 0) {
1228  auto remainingPointsContribution = mapFunction(numVectors);
1229  // Add the contribution from the valid remaining points and store the result in the output variable
1230  auto remainingMask = vecCore::Int2Mask<T>(remainingPoints);
1231  for (unsigned int param = 0; param < npar; param++) {
1232  vecCore::MaskedAssign(gVec[param], remainingMask, gVec[param] + remainingPointsContribution[param]);
1233  }
1234  }
1235  // reduce final gradient result from T to double
1236  for (unsigned int param = 0; param < npar; param++) {
1237  grad[param] = vecCore::ReduceAdd(gVec[param]);
1238  }
1239 
1240 #ifdef DEBUG_FITUTIL
1241  std::cout << "***** Final gradient : ";
1242  for (unsigned int ii = 0; ii< npar; ++ii) std::cout << grad[ii] << " ";
1243  std::cout << "\n";
1244 #endif
1245 
1246  }
1247 
1248  static void EvalLogLGradient(const IModelFunctionTempl<T> &f, const UnBinData &data, const double *p,
1249  double *grad, unsigned int &,
1251  unsigned nChunks = 0)
1252  {
1253  // evaluate the gradient of the log likelihood function
1254 
1255  const IGradModelFunctionTempl<T> *fg = dynamic_cast<const IGradModelFunctionTempl<T> *>(&f);
1256  assert(fg != nullptr); // must be called by a grad function
1257 
1258  const IGradModelFunctionTempl<T> &func = *fg;
1259 
1260 
1261  unsigned int npar = func.NPar();
1262  auto vecSize = vecCore::VectorSize<T>();
1263  unsigned initialNPoints = data.Size();
1264  unsigned numVectors = initialNPoints / vecSize;
1265 
1266 #ifdef DEBUG_FITUTIL
1267  std::cout << "\n===> Evaluate Gradient for parameters ";
1268  for (unsigned int ip = 0; ip < npar; ++ip)
1269  std::cout << " " << p[ip];
1270  std::cout << "\n";
1271 #endif
1272 
1273  (const_cast<IGradModelFunctionTempl<T> &>(func)).SetParameters(p);
1274 
1276  const T kdmax2 = vecCore::NumericLimits<T>::Max() / (4 * initialNPoints);
1277 
1278  auto mapFunction = [&](const unsigned int i) {
1279  std::vector<T> gradFunc(npar);
1280  std::vector<T> pointContributionVec(npar);
1281 
1282  T x1;
1283  vecCore::Load<T>(x1, data.GetCoordComponent(i * vecSize, 0));
1284 
1285  const T *x = nullptr;
1286 
1287  unsigned int ndim = data.NDim();
1288  std::vector<T> xc(ndim);
1289  if (ndim > 1) {
1290  xc.resize(ndim);
1291  xc[0] = x1;
1292  for (unsigned int j = 1; j < ndim; ++j)
1293  vecCore::Load<T>(xc[j], data.GetCoordComponent(i * vecSize, j));
1294  x = xc.data();
1295  } else {
1296  x = &x1;
1297  }
1298 
1299 
1300  T fval = func(x, p);
1301  func.ParameterGradient(x, p, &gradFunc[0]);
1302 
1303 #ifdef DEBUG_FITUTIL
1304  if (i < 5 || (i > numVectors-5) ) {
1305  if (ndim > 1) std::cout << i << " x " << x[0] << " y " << x[1] << " gradient " << gradFunc[0] << " " << gradFunc[1] << " " << gradFunc[3] << std::endl;
1306  else std::cout << i << " x " << x[0] << " gradient " << gradFunc[0] << " " << gradFunc[1] << " " << gradFunc[3] << std::endl;
1307  }
1308 #endif
1309 
1310  vecCore::Mask<T> positiveValues = fval > 0;
1311 
1312  for (unsigned int kpar = 0; kpar < npar; ++kpar) {
1313  if (!vecCore::MaskEmpty(positiveValues))
1314  vecCore::MaskedAssign<T>(pointContributionVec[kpar], positiveValues, -1. / fval * gradFunc[kpar]);
1315 
1316  vecCore::Mask<T> nonZeroGradientValues = !positiveValues && gradFunc[kpar] != 0;
1317  if (!vecCore::MaskEmpty(nonZeroGradientValues)) {
1318  T gg = kdmax1 * gradFunc[kpar];
1319  pointContributionVec[kpar] =
1320  vecCore::Blend(nonZeroGradientValues && gg > 0, -vecCore::math::Min(gg, kdmax2),
1321  -vecCore::math::Max(gg, -kdmax2));
1322  }
1323  // if func derivative is zero term is also zero so do not add in g[kpar]
1324  }
1325 
1326  return pointContributionVec;
1327  };
1328 
1329  // Vertically reduce the set of vectors by summing its equally-indexed components
1330  auto redFunction = [&](const std::vector<std::vector<T>> &pointContributions) {
1331  std::vector<T> result(npar);
1332 
1333  for (auto const &pointContributionVec : pointContributions) {
1334  for (unsigned int parameterIndex = 0; parameterIndex < npar; parameterIndex++)
1335  result[parameterIndex] += pointContributionVec[parameterIndex];
1336  }
1337 
1338  return result;
1339  };
1340 
1341  std::vector<T> gVec(npar);
1342  std::vector<double> g(npar);
1343 
1344 #ifndef R__USE_IMT
1345  // to fix compiler warning
1346  (void)nChunks;
1347 
1348  // If IMT is disabled, force the execution policy to the serial case
1349  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
1350  Warning("FitUtil::EvaluateLogLGradient",
1351  "Multithread execution policy requires IMT, which is disabled. Changing "
1352  "to ROOT::Fit::ExecutionPolicy::kSerial.");
1353  executionPolicy = ROOT::Fit::ExecutionPolicy::kSerial;
1354  }
1355 #endif
1356 
1357  if (executionPolicy == ROOT::Fit::ExecutionPolicy::kSerial) {
1359  gVec = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, numVectors), redFunction);
1360  }
1361 #ifdef R__USE_IMT
1362  else if (executionPolicy == ROOT::Fit::ExecutionPolicy::kMultithread) {
1363  auto chunks = nChunks != 0 ? nChunks : setAutomaticChunking(numVectors);
1364  ROOT::TThreadExecutor pool;
1365  gVec = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, numVectors), redFunction, chunks);
1366  }
1367 #endif
1368  else {
1369  Error("FitUtil::EvaluateLogLGradient", "Execution policy unknown. Avalaible choices:\n "
1370  "ROOT::Fit::ExecutionPolicy::kSerial (default)\n "
1371  "ROOT::Fit::ExecutionPolicy::kMultithread (requires IMT)\n");
1372  }
1373 
1374  // Compute the contribution from the remaining points
1375  unsigned int remainingPoints = initialNPoints % vecSize;
1376  if (remainingPoints > 0) {
1377  auto remainingPointsContribution = mapFunction(numVectors);
1378  // Add the contribution from the valid remaining points and store the result in the output variable
1379  auto remainingMask = vecCore::Int2Mask<T>(initialNPoints % vecSize);
1380  for (unsigned int param = 0; param < npar; param++) {
1381  vecCore::MaskedAssign(gVec[param], remainingMask, gVec[param] + remainingPointsContribution[param]);
1382  }
1383  }
1384  // reduce final gradient result from T to double
1385  for (unsigned int param = 0; param < npar; param++) {
1386  grad[param] = vecCore::ReduceAdd(gVec[param]);
1387  }
1388 
1389 #ifdef DEBUG_FITUTIL
1390  std::cout << "Final gradient ";
1391  for (unsigned int param = 0; param < npar; param++) {
1392  std::cout << " " << grad[param];
1393  }
1394  std::cout << "\n";
1395 #endif
1396  }
1397  };
1398 
1399  template<>
1400  struct Evaluate<double>{
1401 #endif
1402 
1403  static double EvalChi2(const IModelFunction &func, const BinData &data, const double *p, unsigned int &nPoints,
1404  ::ROOT::Fit::ExecutionPolicy executionPolicy, unsigned nChunks = 0)
1405  {
1406  // evaluate the chi2 given a function reference, the data and returns the value and also in nPoints
1407  // the actual number of used points
1408  // normal chi2 using only error on values (from fitting histogram)
1409  // optionally the integral of function in the bin is used
1410 
1411 
1412  //Info("EvalChi2","Using non-vecorized implementation %d",(int) data.Opt().fIntegral);
1413 
1414  return FitUtil::EvaluateChi2(func, data, p, nPoints, executionPolicy, nChunks);
1415  }
1416 
1417  static double EvalLogL(const IModelFunctionTempl<double> &func, const UnBinData &data, const double *p,
1418  int iWeight, bool extended, unsigned int &nPoints,
1419  ::ROOT::Fit::ExecutionPolicy executionPolicy, unsigned nChunks = 0)
1420  {
1421  return FitUtil::EvaluateLogL(func, data, p, iWeight, extended, nPoints, executionPolicy, nChunks);
1422  }
1423 
1424  static double EvalPoissonLogL(const IModelFunctionTempl<double> &func, const BinData &data, const double *p,
1425  int iWeight, bool extended, unsigned int &nPoints,
1426  ::ROOT::Fit::ExecutionPolicy executionPolicy, unsigned nChunks = 0)
1427  {
1428  return FitUtil::EvaluatePoissonLogL(func, data, p, iWeight, extended, nPoints, executionPolicy, nChunks);
1429  }
1430 
1431  static double EvalChi2Effective(const IModelFunctionTempl<double> &func, const BinData & data, const double * p, unsigned int &nPoints)
1432  {
1433  return FitUtil::EvaluateChi2Effective(func, data, p, nPoints);
1434  }
1435  static void EvalChi2Gradient(const IModelFunctionTempl<double> &func, const BinData &data, const double *p,
1436  double *g, unsigned int &nPoints,
1438  unsigned nChunks = 0)
1439  {
1440  FitUtil::EvaluateChi2Gradient(func, data, p, g, nPoints, executionPolicy, nChunks);
1441  }
1442  static double EvalChi2Residual(const IModelFunctionTempl<double> &func, const BinData & data, const double * p, unsigned int i, double *g = 0)
1443  {
1444  return FitUtil::EvaluateChi2Residual(func, data, p, i, g);
1445  }
1446 
1447  /// evaluate the pdf (Poisson) contribution to the logl (return actually log of pdf)
1448  /// and its gradient
1449  static double EvalPoissonBinPdf(const IModelFunctionTempl<double> &func, const BinData & data, const double *p, unsigned int i, double *g ) {
1450  return FitUtil::EvaluatePoissonBinPdf(func, data, p, i, g);
1451  }
1452 
1453  static void
1454  EvalPoissonLogLGradient(const IModelFunctionTempl<double> &func, const BinData &data, const double *p, double *g,
1455  unsigned int &nPoints,
1457  unsigned nChunks = 0)
1458  {
1459  FitUtil::EvaluatePoissonLogLGradient(func, data, p, g, nPoints, executionPolicy, nChunks);
1460  }
1461 
1462  static void EvalLogLGradient(const IModelFunctionTempl<double> &func, const UnBinData &data, const double *p,
1463  double *g, unsigned int &nPoints,
1465  unsigned nChunks = 0)
1466  {
1467  FitUtil::EvaluateLogLGradient(func, data, p, g, nPoints, executionPolicy, nChunks);
1468  }
1469  };
1470 
1471 } // end namespace FitUtil
1472 
1473  } // end namespace Fit
1474 
1475 } // end namespace ROOT
1476 
1477 #if defined (R__HAS_VECCORE) && defined(R__HAS_VC)
1478 //Fixes alignment for structures of SIMD structures
1480 #endif
1481 
1482 #endif /* ROOT_Fit_FitUtil */
LikelihoodAux & operator+=(const LikelihoodAux &l)
Definition: FitUtil.h:84
void SetFunction(const ParamFunc &func, const double *p=0)
Definition: FitUtil.h:139
static double EvalChi2Residual(const IModelFunctionTempl< double > &func, const BinData &data, const double *p, unsigned int i, double *g=0)
Definition: FitUtil.h:1442
ROOT::Math::IntegratorOneDim * fIg1Dim
Definition: FitUtil.h:256
LikelihoodAux operator+(const LikelihoodAux &l) const
Definition: FitUtil.h:79
float xmin
Definition: THbookFile.cxx:93
Interface (abstract class) for generic functions objects of one-dimension Provides a method to evalua...
Definition: IFunction.h:135
void SetParameters(const double *p)
Definition: FitUtil.h:168
unsigned int NDim() const
return coordinate data dimension
Definition: UnBinData.h:280
double ExecFunc(T *f, const double *x, const double *p) const
Definition: FitUtil.h:222
auto * m
Definition: textangle.C:8
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
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
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:1126
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:1403
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
ROOT::Math::IMultiGenFunction * fFuncNDim
Definition: FitUtil.h:259
ErrorType GetErrorType() const
retrieve the errortype
Definition: BinData.h:545
auto MapReduce(F func, unsigned nTimes, R redfunc) -> typename std::result_of< F()>::type
double T(double x)
Definition: ChebyshevPol.h:34
#define g(i)
Definition: RSha256.hxx:105
const double * ValuePtr(unsigned int ipoint) const
return a pointer to the value for the given fit point
Definition: BinData.h:229
unsigned setAutomaticChunking(unsigned nEvents)
Definition: FitUtil.cxx:1756
Class describing the unbinned data sets (just x coordinates values) of any dimensions.
Definition: UnBinData.h:42
#define f(i)
Definition: RSha256.hxx:104
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:168
Interface (abstract class) for parametric gradient multi-dimensional functions providing in addition ...
R__EXTERN TVirtualMutex * gROOTMutex
Definition: TROOT.h:57
ROOT::Math::IParamMultiGradFunction IGradModelFunction
Definition: FitUtil.h:65
ROOT::Math::IntegratorMultiDim * fIgNDim
Definition: FitUtil.h:257
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:1449
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
LikelihoodAux(double logv=0.0, double w=0.0, double w2=0.0)
Definition: FitUtil.h:100
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:1454
static const double x2[5]
Double_t x[n]
Definition: legend1.C:17
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
const double * GetCoordComponent(unsigned int ipoint, unsigned int icoord) const
returns a single coordinate component of a point.
Definition: FitData.h:229
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:1572
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:1462
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:1424
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:1371
double FN(const double *x) const
Definition: FitUtil.h:194
#define F1(x, y, z)
Definition: TMD5.cxx:267
unsigned int Size() const
return number of fit points
Definition: FitData.h:303
#define MATH_ERROR_MSG(loc, str)
Definition: Error.h:50
void Error(const char *location, const char *msgfmt,...)
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 Double_v
Definition: Types.h:50
Documentation for the abstract class IBaseFunctionMultiDim.
Definition: IFunction.h:62
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
This class provides a simple interface to execute the same task multiple times in parallel...
virtual unsigned int NPar() const =0
Return the number of Parameters.
Double_t Sqrt(Double_t x)
DataOptions : simple structure holding the options on how the data are filled.
Definition: DataOptions.h:28
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
LikelihoodAux & operator+=(const LikelihoodAux &l)
Definition: FitUtil.h:107
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
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
IntegralEvaluator(const ParamFunc &func, const double *p, bool useIntegral=true)
Definition: FitUtil.h:131
User Class for performing numerical integration of a function in one dimension.
Definition: Integrator.h:94
LikelihoodAux operator+(const LikelihoodAux &l) const
Definition: FitUtil.h:102
const DataOptions & Opt() const
access to options
Definition: FitData.h:319
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 ...
Class describing the binned data sets : vectors of x coordinates, y values and optionally error on y ...
Definition: BinData.h:53
float xmax
Definition: THbookFile.cxx:93
bool HaveCoordErrors() const
flag to control if data provides error on the coordinates
Definition: BinData.h:131
void Warning(const char *location, const char *msgfmt,...)
ROOT::Math::IParamMultiFunction IModelFunction
Definition: FitUtil.h:64
ROOT::Math::IGenFunction * fFunc1Dim
Definition: FitUtil.h:258
IParamFunction interface (abstract class) describing multi-dimensional parameteric functions It is a ...
virtual void SetParameters(const double *p)=0
Set the parameter values.
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
static const double x1[5]
A pseudo container class which is a generator of indices.
Definition: TSeq.hxx:66
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
double F1(double x) const
Definition: FitUtil.h:188
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
static double EvalChi2Effective(const IModelFunctionTempl< double > &func, const BinData &data, const double *p, unsigned int &nPoints)
Definition: FitUtil.h:1431
Double_t y[n]
Definition: legend1.C:17
#define R__LOCKGUARD(mutex)
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.
Binding & operator=(OUT(*fun)(void))
double SumOfContent() const
compute the total sum of the data content (sum of weights in case of weighted data set) ...
Definition: BinData.h:554
LikelihoodAux(T logv={}, T w={}, T w2={})
Definition: FitUtil.h:77
typedef void((*Func_t)())
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.
User class for performing multidimensional integration.
auto * l
Definition: textangle.C:4
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:200
unsigned int NDim() const
return coordinate data dimension
Definition: FitData.h:311
const DataRange & Range() const
access to range
Definition: FitData.h:331
double Integral(const double *x1, const double *x2)
Definition: FitUtil.h:196
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:1435
const Int_t n
Definition: legend1.C:16
Template class to wrap any member function of a class taking a double and returning a double in a 1D ...
double operator()(const double *x1, const double *x2)
Definition: FitUtil.h:202
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:1417
const double * WeightsPtr(unsigned int ipoint) const
Definition: UnBinData.h:265