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