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