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