ROOT   6.10/09 Reference Guide
FitUtil.cxx
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 // Implementation file for class FitUtil
12
13 #include "Fit/FitUtil.h"
14
15 #include "Fit/BinData.h"
16 #include "Fit/UnBinData.h"
17
18 #include "Math/IFunctionfwd.h"
19 #include "Math/IParamFunction.h"
20 #include "Math/IParamFunctionfwd.h"
21 #include "Math/Integrator.h"
23 #include "Math/WrappedFunction.h"
26
27 #include "Math/Error.h"
28 #include "Math/Util.h" // for safe log(x)
29
30 #include <limits>
31 #include <cmath>
32 #include <cassert>
33 #include <algorithm>
34 //#include <memory>
35
36 //#define DEBUG
37 #ifdef DEBUG
38 #define NSAMPLE 10
39 #include <iostream>
40 #endif
41
42 // using parameter cache is not thread safe but needed for normalizing the functions
43 #define USE_PARAMCACHE
44
45 // need to implement integral option
46
47 namespace ROOT {
48
49  namespace Fit {
50
51  namespace FitUtil {
52
53  // internal class to evaluate the function or the integral
54  // and cached internal integration details
55  // if useIntegral is false no allocation is done
56  // and this is a dummy class
57  // class is templated on any parametric functor implementing operator()(x,p) and NDim()
58  // contains a constant pointer to the function
59  template <class ParamFunc = ROOT::Math::IParamMultiFunction>
60  class IntegralEvaluator {
61
62  public:
63
64  IntegralEvaluator(const ParamFunc & func, const double * p, bool useIntegral = true) :
65  fDim(0),
66  fParams(0),
67  fFunc(0),
68  fIg1Dim(0),
69  fIgNDim(0),
70  fFunc1Dim(0),
71  fFuncNDim(0)
72  {
73  if (useIntegral) {
74  SetFunction(func, p);
75  }
76  }
77
78  void SetFunction(const ParamFunc & func, const double * p = 0) {
79  // set the integrand function and create required wrapper
80  // to perform integral in (x) of a generic f(x,p)
81  fParams = p;
82  fDim = func.NDim();
83  // copy the function object to be able to modify the parameters
84  //fFunc = dynamic_cast<ROOT::Math::IParamMultiFunction *>( func.Clone() );
85  fFunc = &func;
86  assert(fFunc != 0);
87  // set parameters in function
88  //fFunc->SetParameters(p);
89  if (fDim == 1) {
91  fIg1Dim = new ROOT::Math::IntegratorOneDim();
92  //fIg1Dim->SetFunction( static_cast<const ROOT::Math::IMultiGenFunction & >(*fFunc),false);
93  fIg1Dim->SetFunction( static_cast<const ROOT::Math::IGenFunction &>(*fFunc1Dim) );
94  }
95  else if (fDim > 1) {
97  fIgNDim = new ROOT::Math::IntegratorMultiDim();
98  fIgNDim->SetFunction(*fFuncNDim);
99  }
100  else
101  assert(fDim > 0);
102  }
103
104  void SetParameters(const double *p) {
105  // copy just the pointer
106  fParams = p;
107  }
108
109  ~IntegralEvaluator() {
110  if (fIg1Dim) delete fIg1Dim;
111  if (fIgNDim) delete fIgNDim;
112  if (fFunc1Dim) delete fFunc1Dim;
113  if (fFuncNDim) delete fFuncNDim;
114  //if (fFunc) delete fFunc;
115  }
116
117  // evaluation of integrand function (one-dim)
118  double F1 (double x) const {
119  double xx[1]; xx[0] = x;
120  return (*fFunc)( xx, fParams);
121  }
122  // evaluation of integrand function (multi-dim)
123  double FN(const double * x) const {
124  return (*fFunc)( x, fParams);
125  }
126
127  double Integral(const double *x1, const double * x2) {
128  // return unormalized integral
129  return (fIg1Dim) ? fIg1Dim->Integral( *x1, *x2) : fIgNDim->Integral( x1, x2);
130  }
131  double operator()(const double *x1, const double * x2) {
132  // return normalized integral, divided by bin volume (dx1*dx...*dxn)
133  if (fIg1Dim) {
134  double dV = *x2 - *x1;
135  return fIg1Dim->Integral( *x1, *x2)/dV;
136  }
137  else if (fIgNDim) {
138  double dV = 1;
139  for (unsigned int i = 0; i < fDim; ++i)
140  dV *= ( x2[i] - x1[i] );
141  return fIgNDim->Integral( x1, x2)/dV;
142 // std::cout << " do integral btw x " << x1[0] << " " << x2[0] << " y " << x1[1] << " " << x2[1] << " dV = " << dV << " result = " << result << std::endl;
143 // return result;
144  }
145  else
146  assert(1.); // should never be here
147  return 0;
148  }
149
150  private:
151
152  // objects of this class are not meant to be copied / assigned
153  IntegralEvaluator(const IntegralEvaluator& rhs);
154  IntegralEvaluator& operator=(const IntegralEvaluator& rhs);
155
156  unsigned int fDim;
157  const double * fParams;
158  //ROOT::Math::IParamMultiFunction * fFunc; // copy of function in order to be able to change parameters
159  const ParamFunc * fFunc; // reference to a generic parametric function
162  ROOT::Math::IGenFunction * fFunc1Dim;
163  ROOT::Math::IMultiGenFunction * fFuncNDim;
164  };
165
166
167  // derivative with respect of the parameter to be integrated
169  struct ParamDerivFunc {
170  ParamDerivFunc(const GradFunc & f) : fFunc(f), fIpar(0) {}
171  void SetDerivComponent(unsigned int ipar) { fIpar = ipar; }
172  double operator() (const double *x, const double *p) const {
173  return fFunc.ParameterDerivative( x, p, fIpar );
174  }
175  unsigned int NDim() const { return fFunc.NDim(); }
177  unsigned int fIpar;
178  };
179
180  // simple gradient calculator using the 2 points rule
181
183
184  public:
185  // construct from function and gradient dimension gdim
186  // gdim = npar for parameter gradient
187  // gdim = ndim for coordinate gradients
188  // construct (the param values will be passed later)
189  // one can choose between 2 points rule (1 extra evaluation) istrat=1
190  // or two point rule (2 extra evaluation)
191  // (found 2 points rule does not work correctly - minuit2FitBench fails)
192  SimpleGradientCalculator(int gdim, const IModelFunction & func,double eps = 2.E-8, int istrat = 1) :
193  fEps(eps),
194  fPrecision(1.E-8 ), // sqrt(epsilon)
195  fStrategy(istrat),
196  fN(gdim ),
197  fFunc(func),
198  fVec(std::vector<double>(gdim) ) // this can be probably optimized
199  {}
200
201  // internal method to calculate single partial derivative
202  // assume cached vector fVec is already set
203  double DoParameterDerivative(const double *x, const double *p, double f0, int k) const {
204  double p0 = p[k];
205  double h = std::max( fEps* std::abs(p0), 8.0*fPrecision*(std::abs(p0) + fPrecision) );
206  fVec[k] += h;
207  double deriv = 0;
208  // t.b.d : treat case of infinities
209  //if (fval > - std::numeric_limits<double>::max() && fval < std::numeric_limits<double>::max() )
210  double f1 = fFunc(x, &fVec.front() );
211  if (fStrategy > 1) {
212  fVec[k] = p0 - h;
213  double f2 = fFunc(x, &fVec.front() );
214  deriv = 0.5 * ( f2 - f1 )/h;
215  }
216  else
217  deriv = ( f1 - f0 )/h;
218
219  fVec[k] = p[k]; // restore original p value
220  return deriv;
221  }
222  // number of dimension in x (needed when calculating the integrals)
223  unsigned int NDim() const {
224  return fFunc.NDim();
225  }
226  // number of parameters (needed for grad ccalculation)
227  unsigned int NPar() const {
228  return fFunc.NPar();
229  }
230
231  double ParameterDerivative(const double *x, const double *p, int ipar) const {
232  // fVec are the cached parameter values
233  std::copy(p, p+fN, fVec.begin());
234  double f0 = fFunc(x, p);
235  return DoParameterDerivative(x,p,f0,ipar);
236  }
237
238  // calculate all gradient at point (x,p) knnowing already value f0 (we gain a function eval.)
239  void ParameterGradient(const double * x, const double * p, double f0, double * g) {
240  // fVec are the cached parameter values
241  std::copy(p, p+fN, fVec.begin());
242  for (unsigned int k = 0; k < fN; ++k) {
243  g[k] = DoParameterDerivative(x,p,f0,k);
244  }
245  }
246
247  // calculate gradient w.r coordinate values
248  void Gradient(const double * x, const double * p, double f0, double * g) {
249  // fVec are the cached coordinate values
250  std::copy(x, x+fN, fVec.begin());
251  for (unsigned int k = 0; k < fN; ++k) {
252  double x0 = x[k];
253  double h = std::max( fEps* std::abs(x0), 8.0*fPrecision*(std::abs(x0) + fPrecision) );
254  fVec[k] += h;
255  // t.b.d : treat case of infinities
256  //if (fval > - std::numeric_limits<double>::max() && fval < std::numeric_limits<double>::max() )
257  double f1 = fFunc( &fVec.front(), p );
258  if (fStrategy > 1) {
259  fVec[k] = x0 - h;
260  double f2 = fFunc( &fVec.front(), p );
261  g[k] = 0.5 * ( f2 - f1 )/h;
262  }
263  else
264  g[k] = ( f1 - f0 )/h;
265
266  fVec[k] = x[k]; // restore original x value
267  }
268  }
269
270  private:
271
272  double fEps;
273  double fPrecision;
274  int fStrategy; // strategy in calculation ( =1 use 2 point rule( 1 extra func) , = 2 use r point rule)
275  unsigned int fN; // gradient dimension
276  const IModelFunction & fFunc;
277  mutable std::vector<double> fVec; // cached coordinates (or parameter values in case of gradientpar)
278  };
279
280
281  // function to avoid infinities or nan
282  double CorrectValue(double rval) {
283  // avoid infinities or nan in rval
284  if (rval > - std::numeric_limits<double>::max() && rval < std::numeric_limits<double>::max() )
285  return rval;
286  else if (rval < 0)
287  // case -inf
288  return -std::numeric_limits<double>::max();
289  else
290  // case + inf or nan
291  return + std::numeric_limits<double>::max();
292  }
293  bool CheckValue(double & rval) {
294  if (rval > - std::numeric_limits<double>::max() && rval < std::numeric_limits<double>::max() )
295  return true;
296  else if (rval < 0) {
297  // case -inf
298  rval = -std::numeric_limits<double>::max();
299  return false;
300  }
301  else {
302  // case + inf or nan
303  rval = + std::numeric_limits<double>::max();
304  return false;
305  }
306  }
307
308
309  // calculation of the integral of the gradient functions
310  // for a function providing derivative w.r.t parameters
311  // x1 and x2 defines the integration interval , p the parameters
312  template <class GFunc>
313  void CalculateGradientIntegral(const GFunc & gfunc,
314  const double *x1, const double * x2, const double * p, double *g) {
315
316  // needs to calculate the integral for each partial derivative
317  ParamDerivFunc<GFunc> pfunc( gfunc);
318  IntegralEvaluator<ParamDerivFunc<GFunc> > igDerEval( pfunc, p, true);
319  // loop on the parameters
320  unsigned int npar = gfunc.NPar();
321  for (unsigned int k = 0; k < npar; ++k ) {
322  pfunc.SetDerivComponent(k);
323  g[k] = igDerEval( x1, x2 );
324  }
325  }
326
327
328
329  } // end namespace FitUtil
330
331
332
333 //___________________________________________________________________________________________________________________________
334 // for chi2 functions
335 //___________________________________________________________________________________________________________________________
336
337 double FitUtil::EvaluateChi2(const IModelFunction & func, const BinData & data, const double * p, unsigned int & nPoints) {
338  // evaluate the chi2 given a function reference , the data and returns the value and also in nPoints
339  // the actual number of used points
340  // normal chi2 using only error on values (from fitting histogram)
341  // optionally the integral of function in the bin is used
342
343  unsigned int n = data.Size();
344
345  double chi2 = 0;
346  nPoints = 0; // count the effective non-zero points
347  // set parameters of the function to cache integral value
348 #ifdef USE_PARAMCACHE
349  (const_cast<IModelFunction &>(func)).SetParameters(p);
350 #endif
351  // do not cache parameter values (it is not thread safe)
352  //func.SetParameters(p);
353
354
355  // get fit option and check case if using integral of bins
356  const DataOptions & fitOpt = data.Opt();
357  bool useBinIntegral = fitOpt.fIntegral && data.HasBinEdges();
358  bool useBinVolume = (fitOpt.fBinVolume && data.HasBinEdges());
359  bool useExpErrors = (fitOpt.fExpErrors);
360
361 #ifdef DEBUG
362  std::cout << "\n\nFit data size = " << n << std::endl;
363  std::cout << "evaluate chi2 using function " << &func << " " << p << std::endl;
364  std::cout << "use empty bins " << fitOpt.fUseEmpty << std::endl;
365  std::cout << "use integral " << fitOpt.fIntegral << std::endl;
366  std::cout << "use all error=1 " << fitOpt.fErrors1 << std::endl;
367 #endif
368
369 #ifdef USE_PARAMCACHE
370  IntegralEvaluator<> igEval( func, 0, useBinIntegral);
371 #else
372  IntegralEvaluator<> igEval( func, p, useBinIntegral);
373 #endif
374  double maxResValue = std::numeric_limits<double>::max() /n;
375  double wrefVolume = 1.0;
376  std::vector<double> xc;
377  if (useBinVolume) {
378  if (fitOpt.fNormBinVolume) wrefVolume /= data.RefVolume();
379  xc.resize(data.NDim() );
380  }
381
382  (const_cast<IModelFunction &>(func)).SetParameters(p);
383  for (unsigned int i = 0; i < n; ++ i) {
384
385  double y = 0, invError = 1.;
386
387  // in case of no error in y invError=1 is returned
388  const double * x1 = data.GetPoint(i,y, invError);
389
390  double fval = 0;
391
392  double binVolume = 1.0;
393  if (useBinVolume) {
394  unsigned int ndim = data.NDim();
395  const double * x2 = data.BinUpEdge(i);
396  for (unsigned int j = 0; j < ndim; ++j) {
397  binVolume *= std::abs( x2[j]-x1[j] );
398  xc[j] = 0.5*(x2[j]+ x1[j]);
399  }
400  // normalize the bin volume using a reference value
401  binVolume *= wrefVolume;
402  }
403
404  const double * x = (useBinVolume) ? &xc.front() : x1;
405
406  if (!useBinIntegral) {
407 #ifdef USE_PARAMCACHE
408  fval = func ( x );
409 #else
410  fval = func ( x, p );
411 #endif
412  }
413  else {
414  // calculate integral normalized by bin volume
415  // need to set function and parameters here in case loop is parallelized
416  fval = igEval( x1, data.BinUpEdge(i)) ;
417  }
418  // normalize result if requested according to bin volume
419  if (useBinVolume) fval *= binVolume;
420
421  // expected errors
422  if (useExpErrors) {
423  // we need first to check if a weight factor needs to be applied
424  // weight = sumw2/sumw = error**2/content
425  double invWeight = y * invError * invError;
426  // if (invError == 0) invWeight = (data.SumOfError2() > 0) ? data.SumOfContent()/ data.SumOfError2() : 1.0;
427  // compute expected error as f(x) / weight
428  double invError2 = (fval > 0) ? invWeight / fval : 0.0;
429  invError = std::sqrt(invError2);
430  }
431
432 //#define DEBUG
433 #ifdef DEBUG
434  std::cout << x[0] << " " << y << " " << 1./invError << " params : ";
435  for (unsigned int ipar = 0; ipar < func.NPar(); ++ipar)
436  std::cout << p[ipar] << "\t";
437  std::cout << "\tfval = " << fval << " bin volume " << binVolume << " ref " << wrefVolume << std::endl;
438 #endif
439 //#undef DEBUG
440
441
442  if (invError > 0) {
443  nPoints++;
444
445  double tmp = ( y -fval )* invError;
446  double resval = tmp * tmp;
447
448
449  // avoid inifinity or nan in chi2 values due to wrong function values
450  if ( resval < maxResValue )
451  chi2 += resval;
452  else {
453  //nRejected++;
454  chi2 += maxResValue;
455  }
456  }
457
458
459  }
460  nPoints=n;
461
462 #ifdef DEBUG
463  std::cout << "chi2 = " << chi2 << " n = " << nPoints /*<< " rejected = " << nRejected */ << std::endl;
464 #endif
465
466
467  return chi2;
468 }
469
470
471 //___________________________________________________________________________________________________________________________
472
473 double FitUtil::EvaluateChi2Effective(const IModelFunction & func, const BinData & data, const double * p, unsigned int & nPoints) {
474  // evaluate the chi2 given a function reference , the data and returns the value and also in nPoints
475  // the actual number of used points
476  // method using the error in the coordinates
477  // integral of bin does not make sense in this case
478
479  unsigned int n = data.Size();
480
481 #ifdef DEBUG
482  std::cout << "\n\nFit data size = " << n << std::endl;
483  std::cout << "evaluate effective chi2 using function " << &func << " " << p << std::endl;
484 #endif
485
486  assert(data.HaveCoordErrors() || data.HaveAsymErrors());
487
488  double chi2 = 0;
489  //int nRejected = 0;
490
491
492  //func.SetParameters(p);
493
494  unsigned int ndim = func.NDim();
495
496  // use Richardson derivator
498
499  double maxResValue = std::numeric_limits<double>::max() /n;
500
501
502
503  for (unsigned int i = 0; i < n; ++ i) {
504
505
506  double y = 0;
507  const double * x = data.GetPoint(i,y);
508
509  double fval = func( x, p );
510
511  double delta_y_func = y - fval;
512
513
514  double ey = 0;
515  const double * ex = 0;
516  if (!data.HaveAsymErrors() )
517  ex = data.GetPointError(i, ey);
518  else {
519  double eylow, eyhigh = 0;
520  ex = data.GetPointError(i, eylow, eyhigh);
521  if ( delta_y_func < 0)
522  ey = eyhigh; // function is higher than points
523  else
524  ey = eylow;
525  }
526  double e2 = ey * ey;
527  // before calculating the gradient check that all error in x are not zero
528  unsigned int j = 0;
529  while ( j < ndim && ex[j] == 0.) { j++; }
530  // if j is less ndim some elements are not zero
531  if (j < ndim) {
532  // need an adapter from a multi-dim function to a one-dimensional
534  // select optimal step size (use 10--2 by default as was done in TF1:
535  double kEps = 0.01;
536  double kPrecision = 1.E-8;
537  for (unsigned int icoord = 0; icoord < ndim; ++icoord) {
538  // calculate derivative for each coordinate
539  if (ex[icoord] > 0) {
541  f1D.SetCoord(icoord);
542  // optimal spep size (take ex[] as scale for the points and 1% of it
543  double x0= x[icoord];
544  double h = std::max( kEps* std::abs(ex[icoord]), 8.0*kPrecision*(std::abs(x0) + kPrecision) );
545  double deriv = derivator.Derivative1(f1D, x[icoord], h);
546  double edx = ex[icoord] * deriv;
547  e2 += edx * edx;
548 #ifdef DEBUG
549  std::cout << "error for coord " << icoord << " = " << ex[icoord] << " deriv " << deriv << std::endl;
550 #endif
551  }
552  }
553  }
554  double w2 = (e2 > 0) ? 1.0/e2 : 0;
555  double resval = w2 * ( y - fval ) * ( y - fval);
556
557 #ifdef DEBUG
558  std::cout << x[0] << " " << y << " ex " << ex[0] << " ey " << ey << " params : ";
559  for (unsigned int ipar = 0; ipar < func.NPar(); ++ipar)
560  std::cout << p[ipar] << "\t";
561  std::cout << "\tfval = " << fval << "\tresval = " << resval << std::endl;
562 #endif
563
564  // avoid (infinity and nan ) in the chi2 sum
565  // eventually add possibility of excluding some points (like singularity)
566  if ( resval < maxResValue )
567  chi2 += resval;
568  else
569  chi2 += maxResValue;
570  //nRejected++;
571
572  }
573
574  // reset the number of fitting data points
575  nPoints = n; // no points are rejected
576  //if (nRejected != 0) nPoints = n - nRejected;
577
578 #ifdef DEBUG
579  std::cout << "chi2 = " << chi2 << " n = " << nPoints << std::endl;
580 #endif
581
582  return chi2;
583
584 }
585
586
587 ////////////////////////////////////////////////////////////////////////////////
588 /// evaluate the chi2 contribution (residual term) only for data with no coord-errors
589 /// This function is used in the specialized least square algorithms like FUMILI or L.M.
590 /// if we have error on the coordinates the method is not yet implemented
591 /// integral option is also not yet implemented
592 /// one can use in that case normal chi2 method
593
594 double FitUtil::EvaluateChi2Residual(const IModelFunction & func, const BinData & data, const double * p, unsigned int i, double * g) {
595  if (data.GetErrorType() == BinData::kCoordError && data.Opt().fCoordErrors ) {
596  MATH_ERROR_MSG("FitUtil::EvaluateChi2Residual","Error on the coordinates are not used in calculating Chi2 residual");
597  return 0; // it will assert otherwise later in GetPoint
598  }
599
600
601  //func.SetParameters(p);
602
603  double y, invError = 0;
604
605  const DataOptions & fitOpt = data.Opt();
606  bool useBinIntegral = fitOpt.fIntegral && data.HasBinEdges();
607  bool useBinVolume = (fitOpt.fBinVolume && data.HasBinEdges());
608  bool useExpErrors = (fitOpt.fExpErrors);
609
610  const double * x1 = data.GetPoint(i,y, invError);
611
612  IntegralEvaluator<> igEval( func, p, useBinIntegral);
613  double fval = 0;
614  unsigned int ndim = data.NDim();
615  double binVolume = 1.0;
616  const double * x2 = 0;
617  if (useBinVolume || useBinIntegral) x2 = data.BinUpEdge(i);
618
619  double * xc = 0;
620
621  if (useBinVolume) {
622  xc = new double[ndim];
623  for (unsigned int j = 0; j < ndim; ++j) {
624  binVolume *= std::abs( x2[j]-x1[j] );
625  xc[j] = 0.5*(x2[j]+ x1[j]);
626  }
627  // normalize the bin volume using a reference value
628  binVolume /= data.RefVolume();
629  }
630
631  const double * x = (useBinVolume) ? xc : x1;
632
633  if (!useBinIntegral) {
634  fval = func ( x, p );
635  }
636  else {
637  // calculate integral (normalized by bin volume)
638  // need to set function and parameters here in case loop is parallelized
639  fval = igEval( x1, x2) ;
640  }
641  // normalize result if requested according to bin volume
642  if (useBinVolume) fval *= binVolume;
643
644  // expected errors
645  if (useExpErrors) {
646  // we need first to check if a weight factor needs to be applied
647  // weight = sumw2/sumw = error**2/content
648  double invWeight = y * invError * invError;
649  // if (invError == 0) invWeight = (data.SumOfError2() > 0) ? data.SumOfContent()/ data.SumOfError2() : 1.0;
650  // compute expected error as f(x) / weight
651  double invError2 = (fval > 0) ? invWeight / fval : 0.0;
652  invError = std::sqrt(invError2);
653  }
654
655
656  double resval = ( y -fval )* invError;
657
658  // avoid infinities or nan in resval
659  resval = CorrectValue(resval);
660
662  if (g != 0) {
663
664  unsigned int npar = func.NPar();
666
667  if (gfunc != 0) {
669  if (!useBinIntegral ) {
670  gfunc->ParameterGradient( x , p, g );
671  }
672  else {
673  // needs to calculate the integral for each partial derivative
674  CalculateGradientIntegral( *gfunc, x1, x2, p, g);
675  }
676  }
677  else {
679  if (!useBinIntegral )
681  else {
682  // needs to calculate the integral for each partial derivative
683  CalculateGradientIntegral( gc, x1, x2, p, g);
684  }
685  }
686  // mutiply by - 1 * weight
687  for (unsigned int k = 0; k < npar; ++k) {
688  g[k] *= - invError;
689  if (useBinVolume) g[k] *= binVolume;
690  }
691  }
692
693  if (useBinVolume) delete [] xc;
694
695  return resval;
696
697 }
698
699 void FitUtil::EvaluateChi2Gradient(const IModelFunction & f, const BinData & data, const double * p, double * grad, unsigned int & nPoints) {
700  // evaluate the gradient of the chi2 function
701  // this function is used when the model function knows how to calculate the derivative and we can
702  // avoid that the minimizer re-computes them
703  //
704  // case of chi2 effective (errors on coordinate) is not supported
705
706  if ( data.HaveCoordErrors() ) {
707  MATH_ERROR_MSG("FitUtil::EvaluateChi2Residual","Error on the coordinates are not used in calculating Chi2 gradient"); return; // it will assert otherwise later in GetPoint
708  }
709
710  unsigned int nRejected = 0;
711
713  assert (fg != 0); // must be called by a gradient function
714
715  const IGradModelFunction & func = *fg;
716  unsigned int n = data.Size();
717
718
719 #ifdef DEBUG
720  std::cout << "\n\nFit data size = " << n << std::endl;
721  std::cout << "evaluate chi2 using function gradient " << &func << " " << p << std::endl;
722 #endif
723
724  const DataOptions & fitOpt = data.Opt();
725  bool useBinIntegral = fitOpt.fIntegral && data.HasBinEdges();
726  bool useBinVolume = (fitOpt.fBinVolume && data.HasBinEdges());
727
728  double wrefVolume = 1.0;
729  std::vector<double> xc;
730  if (useBinVolume) {
731  if (fitOpt.fNormBinVolume) wrefVolume /= data.RefVolume();
732  xc.resize(data.NDim() );
733  }
734
735  IntegralEvaluator<> igEval( func, p, useBinIntegral);
736
737  //int nRejected = 0;
738  // set values of parameters
739
740  unsigned int npar = func.NPar();
741  // assert (npar == NDim() ); // npar MUST be Chi2 dimension
743  // set all vector values to zero
744  std::vector<double> g( npar);
745
746  for (unsigned int i = 0; i < n; ++ i) {
747
748
749  double y, invError = 0;
750  const double * x1 = data.GetPoint(i,y, invError);
751
752  double fval = 0;
753  const double * x2 = 0;
754
755  double binVolume = 1;
756  if (useBinVolume) {
757  unsigned int ndim = data.NDim();
758  x2 = data.BinUpEdge(i);
759  for (unsigned int j = 0; j < ndim; ++j) {
760  binVolume *= std::abs( x2[j]-x1[j] );
761  xc[j] = 0.5*(x2[j]+ x1[j]);
762  }
763  // normalize the bin volume using a reference value
764  binVolume *= wrefVolume;
765  }
766
767  const double * x = (useBinVolume) ? &xc.front() : x1;
768
769  if (!useBinIntegral ) {
770  fval = func ( x, p );
772  }
773  else {
774  x2 = data.BinUpEdge(i);
775  // calculate normalized integral and gradient (divided by bin volume)
776  // need to set function and parameters here in case loop is parallelized
777  fval = igEval( x1, x2 ) ;
779  }
780  if (useBinVolume) fval *= binVolume;
781
782 #ifdef DEBUG
783  std::cout << x[0] << " " << y << " " << 1./invError << " params : ";
784  for (unsigned int ipar = 0; ipar < npar; ++ipar)
785  std::cout << p[ipar] << "\t";
786  std::cout << "\tfval = " << fval << std::endl;
787 #endif
788  if ( !CheckValue(fval) ) {
789  nRejected++;
790  continue;
791  }
792
793  // loop on the parameters
794  unsigned int ipar = 0;
795  for ( ; ipar < npar ; ++ipar) {
796
797  // correct gradient for bin volumes
798  if (useBinVolume) gradFunc[ipar] *= binVolume;
799
800  // avoid singularity in the function (infinity and nan ) in the chi2 sum
801  // eventually add possibility of excluding some points (like singularity)
803  if ( !CheckValue(dfval) ) {
804  break; // exit loop on parameters
805  }
806
807  // calculate derivative point contribution
808  double tmp = - 2.0 * ( y -fval )* invError * invError * gradFunc[ipar];
809  g[ipar] += tmp;
810
811  }
812
813  if ( ipar < npar ) {
814  // case loop was broken for an overflow in the gradient calculation
815  nRejected++;
816  continue;
817  }
818
819
820  }
821
822  // correct the number of points
823  nPoints = n;
824  if (nRejected != 0) {
825  assert(nRejected <= n);
826  nPoints = n - nRejected;
827  if (nPoints < npar) MATH_ERROR_MSG("FitUtil::EvaluateChi2Gradient","Error - too many points rejected for overflow in gradient calculation");
828  }
829
830  // copy result
832
833 }
834
835 //______________________________________________________________________________________________________
836 //
837 // Log Likelihood functions
838 //_______________________________________________________________________________________________________
839
840 // utility function used by the likelihoods
841
842 // for LogLikelihood functions
843
844 double FitUtil::EvaluatePdf(const IModelFunction & func, const UnBinData & data, const double * p, unsigned int i, double * g) {
845  // evaluate the pdf contribution to the generic logl function in case of bin data
846  // return actually the log of the pdf and its derivatives
847
848
849  //func.SetParameters(p);
850
851
852  const double * x = data.Coords(i);
853  double fval = func ( x, p );
854  double logPdf = ROOT::Math::Util::EvalLog(fval);
855  //return
856  if (g == 0) return logPdf;
857
859
861  if (gfunc != 0) {
863  gfunc->ParameterGradient( x , p, g );
864  }
865  else {
866  // estimate gradieant numerically with simple 2 point rule
867  // should probably calculate gradient of log(pdf) is more stable numerically
869  gc.ParameterGradient(x, p, fval, g );
870  }
871  // divide gradient by function value since returning the logs
872  for (unsigned int ipar = 0; ipar < func.NPar(); ++ipar) {
873  g[ipar] /= fval; // this should be checked against infinities
874  }
875
876 #ifdef DEBUG
877  std::cout << x[i] << "\t";
878  std::cout << "\tpar = [ " << func.NPar() << " ] = ";
879  for (unsigned int ipar = 0; ipar < func.NPar(); ++ipar)
880  std::cout << p[ipar] << "\t";
881  std::cout << "\tfval = " << fval;
882  std::cout << "\tgrad = [ ";
883  for (unsigned int ipar = 0; ipar < func.NPar(); ++ipar)
884  std::cout << g[ipar] << "\t";
885  std::cout << " ] " << std::endl;
886 #endif
887
888
889  return logPdf;
890 }
891
892 double FitUtil::EvaluateLogL(const IModelFunction & func, const UnBinData & data, const double * p,
893  int iWeight, bool extended, unsigned int &nPoints) {
894  // evaluate the LogLikelihood
895
896  unsigned int n = data.Size();
897
898 #ifdef DEBUG
899  std::cout << "\n\nFit data size = " << n << std::endl;
900  std::cout << "func pointer is " << typeid(func).name() << std::endl;
901 #endif
902
903  double logl = 0;
904  //unsigned int nRejected = 0;
905
906  // set parameters of the function to cache integral value
907 #ifdef USE_PARAMCACHE
908  (const_cast<IModelFunction &>(func)).SetParameters(p);
909 #endif
910
911  // this is needed if function must be normalized
912  bool normalizeFunc = false;
913  double norm = 1.0;
914  if (normalizeFunc) {
915  // compute integral of the function
916  std::vector<double> xmin(data.NDim());
917  std::vector<double> xmax(data.NDim());
918  IntegralEvaluator<> igEval( func, p, true);
919  // compute integral in the ranges where is defined
920  if (data.Range().Size() > 0 ) {
921  norm = 0;
922  for (unsigned int ir = 0; ir < data.Range().Size(); ++ir) {
923  data.Range().GetRange(&xmin[0],&xmax[0],ir);
924  norm += igEval.Integral(xmin.data(),xmax.data());
925  }
926  } else {
927  // use (-inf +inf)
928  data.Range().GetRange(&xmin[0],&xmax[0]);
929  // check if funcition is zero at +- inf
930  if (func(xmin.data()) != 0 || func(xmax.data()) != 0) {
931  MATH_ERROR_MSG("FitUtil::EvaluateLogLikelihood","A range has not been set and the function is not zero at +/- inf");
932  return 0;
933  }
934  norm = igEval.Integral(&xmin[0],&xmax[0]);
935  }
936  }
937
938  // needed to compue effective global weight in case of extended likelihood
939  double sumW = 0;
940  double sumW2 = 0;
941
942  for (unsigned int i = 0; i < n; ++ i) {
943  const double * x = data.Coords(i);
944 #ifdef USE_PARAMCACHE
945  double fval = func ( x );
946 #else
947  double fval = func ( x, p );
948 #endif
949  if (normalizeFunc) fval = fval / norm;
950
951 #ifdef DEBUG
952  if (i == 0) {
953  std::cout << "x [ " << data.NDim() << " ] = ";
954  for (unsigned int j = 0; j < data.NDim(); ++j)
955  std::cout << x[j] << "\t";
956  std::cout << "\tpar = [ " << func.NPar() << " ] = ";
957  for (unsigned int ipar = 0; ipar < func.NPar(); ++ipar)
958  std::cout << p[ipar] << "\t";
959  std::cout << "\tfval = " << fval << std::endl;
960  } else {
961  std::cout << ".";
962  }
963 #endif
964  // function EvalLog protects against negative or too small values of fval
965  double logval = ROOT::Math::Util::EvalLog( fval);
966  if (iWeight > 0) {
967  double weight = data.Weight(i);
968  logval *= weight;
969  if (iWeight ==2) {
970  logval *= weight; // use square of weights in likelihood
971  if (extended) {
972  // needed sum of weights and sum of weight square if likelkihood is extended
973  sumW += weight;
974  sumW2 += weight*weight;
975  }
976  }
977  }
978  logl += logval;
979  }
980
981 #ifdef DEBUG
982  std::cout << std::endl;
983 #endif
984
985  if (extended) {
986  // add Poisson extended term
987  double extendedTerm = 0; // extended term in likelihood
988  double nuTot = 0;
989  // nuTot is integral of function in the range
990  // if function has been normalized integral has been already computed
991  if (!normalizeFunc) {
992  IntegralEvaluator<> igEval( func, p, true);
993  std::vector<double> xmin(data.NDim());
994  std::vector<double> xmax(data.NDim());
995
996  // compute integral in the ranges where is defined
997  if (data.Range().Size() > 0 ) {
998  nuTot = 0;
999  for (unsigned int ir = 0; ir < data.Range().Size(); ++ir) {
1000  data.Range().GetRange(&xmin[0],&xmax[0],ir);
1001  nuTot += igEval.Integral(xmin.data(),xmax.data());
1002 #ifdef DEBUG
1003  std::cout << "After Integral bewteen " << xmin[0] << " and " << xmax[0] << " nexp is " << nuTot << std::endl;
1004 #endif
1005  }
1006  } else {
1007  // use (-inf +inf)
1008  data.Range().GetRange(&xmin[0],&xmax[0]);
1009  // check if funcition is zero at +- inf
1010  if (func(xmin.data()) != 0 || func(xmax.data()) != 0) {
1011  MATH_ERROR_MSG("FitUtil::EvaluateLogLikelihood","A range has not been set and the function is not zero at +/- inf");
1012  return 0;
1013  }
1014  nuTot = igEval.Integral(&xmin[0],&xmax[0]);
1015 #ifdef DEBUG
1016  std::cout << "Range not existig - integral bewteen " << xmin[0] << " and " << xmax[0] << " is " << nuTot << std::endl;
1017 #endif
1018  }
1019
1020  // force to be last parameter value
1021  //nutot = p[func.NDim()-1];
1022  if (iWeight != 2)
1023  extendedTerm = - nuTot; // no need to add in this case n log(nu) since is already computed before
1024  else {
1025  // case use weight square in likelihood : compute total effective weight = sw2/sw
1026  // ignore for the moment case when sumW is zero
1027  extendedTerm = - (sumW2 / sumW) * nuTot;
1028  }
1029
1030  }
1031  else {
1032  nuTot = norm;
1033  extendedTerm = - nuTot + double(n) * ROOT::Math::Util::EvalLog( nuTot);
1034  // in case of weights need to use here sum of weights (to be done)
1035  }
1036  logl += extendedTerm;
1037
1038 #ifdef DEBUG
1039  // for (unsigned int ipar = 0; ipar < func.NPar(); ++ipar)
1040  // std::cout << p[ipar] << "\t";
1041  // std::cout << std::endl;
1042  std::cout << "fit is extended n = " << n << " nutot " << nuTot << " extended LL term = " << extendedTerm << " logl = " << logl
1043  << std::endl;
1044 #endif
1045  }
1046
1047  // reset the number of fitting data points
1048  nPoints = n;
1049 // if (nRejected != 0) {
1050 // assert(nRejected <= n);
1051 // nPoints = n - nRejected;
1052 // if ( nPoints < func.NPar() )
1053 // MATH_ERROR_MSG("FitUtil::EvaluateLogL","Error too many points rejected because of bad pdf values");
1054
1055 // }
1056 #ifdef DEBUG
1057  std::cout << "Logl = " << logl << " np = " << nPoints << std::endl;
1058 #endif
1059
1060  return -logl;
1061 }
1062
1063 void FitUtil::EvaluateLogLGradient(const IModelFunction & f, const UnBinData & data, const double * p, double * grad, unsigned int & ) {
1064  // evaluate the gradient of the log likelihood function
1065
1067  assert (fg != 0); // must be called by a grad function
1068  const IGradModelFunction & func = *fg;
1069
1070  unsigned int n = data.Size();
1071  //int nRejected = 0;
1072
1073  unsigned int npar = func.NPar();
1075  std::vector<double> g( npar);
1076
1077  for (unsigned int i = 0; i < n; ++ i) {
1078  const double * x = data.Coords(i);
1079  double fval = func ( x , p);
1081  for (unsigned int kpar = 0; kpar < npar; ++ kpar) {
1082  if (fval > 0)
1083  g[kpar] -= 1./fval * gradFunc[ kpar ];
1084  else if (gradFunc [ kpar] != 0) {
1085  const double kdmax1 = std::sqrt( std::numeric_limits<double>::max() );
1086  const double kdmax2 = std::numeric_limits<double>::max() / (4*n);
1087  double gg = kdmax1 * gradFunc[ kpar ];
1088  if ( gg > 0) gg = std::min( gg, kdmax2);
1089  else gg = std::max(gg, - kdmax2);
1090  g[kpar] -= gg;
1091  }
1092  // if func derivative is zero term is also zero so do not add in g[kpar]
1093  }
1094
1095  // copy result
1097  }
1098 }
1099 //_________________________________________________________________________________________________
1100 // for binned log likelihood functions
1101 ////////////////////////////////////////////////////////////////////////////////
1102 /// evaluate the pdf (Poisson) contribution to the logl (return actually log of pdf)
1104
1105 double FitUtil::EvaluatePoissonBinPdf(const IModelFunction & func, const BinData & data, const double * p, unsigned int i, double * g ) {
1106  double y = 0;
1107  const double * x1 = data.GetPoint(i,y);
1108
1109  const DataOptions & fitOpt = data.Opt();
1110  bool useBinIntegral = fitOpt.fIntegral && data.HasBinEdges();
1111  bool useBinVolume = (fitOpt.fBinVolume && data.HasBinEdges());
1112
1113  IntegralEvaluator<> igEval( func, p, useBinIntegral);
1114  double fval = 0;
1115  const double * x2 = 0;
1116  // calculate the bin volume
1117  double binVolume = 1;
1118  std::vector<double> xc;
1119  if (useBinVolume) {
1120  unsigned int ndim = data.NDim();
1121  xc.resize(ndim);
1122  x2 = data.BinUpEdge(i);
1123  for (unsigned int j = 0; j < ndim; ++j) {
1124  binVolume *= std::abs( x2[j]-x1[j] );
1125  xc[j] = 0.5*(x2[j]+ x1[j]);
1126  }
1127  // normalize the bin volume using a reference value
1128  binVolume /= data.RefVolume();
1129  }
1130
1131  const double * x = (useBinVolume) ? &xc.front() : x1;
1132
1133  if (!useBinIntegral ) {
1134  fval = func ( x, p );
1135  }
1136  else {
1137  // calculate integral normalized (divided by bin volume)
1138  x2 = data.BinUpEdge(i);
1139  fval = igEval( x1, x2 ) ;
1140  }
1141  if (useBinVolume) fval *= binVolume;
1142
1143  // logPdf for Poisson: ignore constant term depending on N
1144  fval = std::max(fval, 0.0); // avoid negative or too small values
1145  double logPdf = - fval;
1146  if (y > 0.0) {
1147  // include also constants due to saturate model (see Baker-Cousins paper)
1148  logPdf += y * ROOT::Math::Util::EvalLog( fval / y) + y;
1149  }
1150  // need to return the pdf contribution (not the log)
1151
1152  //double pdfval = std::exp(logPdf);
1153
1154  //if (g == 0) return pdfval;
1155  if (g == 0) return logPdf;
1156
1157  unsigned int npar = func.NPar();
1159
1161  if (gfunc != 0) {
1163  if (!useBinIntegral )
1164  gfunc->ParameterGradient( x , p, g );
1165  else {
1166  // needs to calculate the integral for each partial derivative
1167  CalculateGradientIntegral( *gfunc, x1, x2, p, g);
1168  }
1169
1170  }
1171  else {
1173  if (!useBinIntegral )
1175  else {
1176  // needs to calculate the integral for each partial derivative
1177  CalculateGradientIntegral( gc, x1, x2, p, g);
1178  }
1179  }
1180  // correct g[] do be derivative of poisson term
1181  for (unsigned int k = 0; k < npar; ++k) {
1182  // apply bin volume correction
1183  if (useBinVolume) g[k] *= binVolume;
1184
1185  // correct for Poisson term
1186  if ( fval > 0)
1187  g[k] *= ( y/fval - 1.) ;//* pdfval;
1188  else if (y > 0) {
1189  const double kdmax1 = std::sqrt( std::numeric_limits<double>::max() );
1190  g[k] *= kdmax1;
1191  }
1192  else // y == 0 cannot have negative y
1193  g[k] *= -1;
1194  }
1195
1196
1197 #ifdef DEBUG
1198  std::cout << "x = " << x[0] << " logPdf = " << logPdf << " grad";
1199  for (unsigned int ipar = 0; ipar < npar; ++ipar)
1200  std::cout << g[ipar] << "\t";
1201  std::cout << std::endl;
1202 #endif
1203
1204 // return pdfval;
1205  return logPdf;
1206 }
1207
1209  const double * p, int iWeight, bool extended, unsigned int & nPoints ) {
1210  // evaluate the Poisson Log Likelihood
1211  // for binned likelihood fits
1212  // this is Sum ( f(x_i) - y_i * log( f (x_i) ) )
1213  // add as well constant term for saturated model to make it like a Chi2/2
1214  // by default is etended. If extended is false the fit is not extended and
1215  // the global poisson term is removed (i.e is a binomial fit)
1216  // (remember that in this case one needs to have a function with a fixed normalization
1217  // like in a non extended binned fit)
1218  //
1219  // if use Weight use a weighted dataset
1220  // iWeight = 1 ==> logL = Sum( w f(x_i) )
1221  // case of iWeight==1 is actually identical to weight==0
1222  // iWeight = 2 ==> logL = Sum( w*w * f(x_i) )
1223  //
1224  // nPoints returns the points where bin content is not zero
1225
1226
1227  unsigned int n = data.Size();
1228
1229 #ifdef USE_PARAMCACHE
1230  (const_cast<IModelFunction &>(func)).SetParameters(p);
1231 #endif
1232
1233  double nloglike = 0; // negative loglikelihood
1234  nPoints = 0; // npoints
1235
1236
1237  // get fit option and check case of using integral of bins
1238  const DataOptions & fitOpt = data.Opt();
1239  bool useBinIntegral = fitOpt.fIntegral && data.HasBinEdges();
1240  bool useBinVolume = (fitOpt.fBinVolume && data.HasBinEdges());
1241  bool useW2 = (iWeight == 2);
1242
1243  // normalize if needed by a reference volume value
1244  double wrefVolume = 1.0;
1245  std::vector<double> xc;
1246  if (useBinVolume) {
1247  if (fitOpt.fNormBinVolume) wrefVolume /= data.RefVolume();
1248  xc.resize(data.NDim() );
1249  }
1250
1251 #ifdef DEBUG
1252  std::cout << "Evaluate PoissonLogL for params = [ ";
1253  for (unsigned int j=0; j < func.NPar(); ++j) std::cout << p[j] << " , ";
1254  std::cout << "] - data size = " << n << " useBinIntegral " << useBinIntegral << " useBinVolume "
1255  << useBinVolume << " useW2 " << useW2 << " wrefVolume = " << wrefVolume << std::endl;
1256 #endif
1257
1258 #ifdef USE_PARAMCACHE
1259  IntegralEvaluator<> igEval( func, 0, useBinIntegral);
1260 #else
1261  IntegralEvaluator<> igEval( func, p, useBinIntegral);
1262 #endif
1263  // double nuTot = 0; // total number of expected events (needed for non-extended fits)
1264  // double wTot = 0; // sum of all weights
1265  // double w2Tot = 0; // sum of weight squared (these are needed for useW2)
1266
1267
1268  for (unsigned int i = 0; i < n; ++ i) {
1269  const double * x1 = data.Coords(i);
1270  double y = data.Value(i);
1271
1272  double fval = 0;
1273  double binVolume = 1.0;
1274
1275  if (useBinVolume) {
1276  unsigned int ndim = data.NDim();
1277  const double * x2 = data.BinUpEdge(i);
1278  for (unsigned int j = 0; j < ndim; ++j) {
1279  binVolume *= std::abs( x2[j]-x1[j] );
1280  xc[j] = 0.5*(x2[j]+ x1[j]);
1281  }
1282  // normalize the bin volume using a reference value
1283  binVolume *= wrefVolume;
1284  }
1285
1286  const double * x = (useBinVolume) ? &xc.front() : x1;
1287
1288  if (!useBinIntegral) {
1289 #ifdef USE_PARAMCACHE
1290  fval = func ( x );
1291 #else
1292  fval = func ( x, p );
1293 #endif
1294  }
1295  else {
1296  // calculate integral (normalized by bin volume)
1297  // need to set function and parameters here in case loop is parallelized
1298  fval = igEval( x1, data.BinUpEdge(i)) ;
1299  }
1300  if (useBinVolume) fval *= binVolume;
1301
1302
1303
1304 #ifdef DEBUG
1305  int NSAMPLE = 100;
1306  if (i%NSAMPLE == 0) {
1307  std::cout << "evt " << i << " x1 = [ ";
1308  for (unsigned int j=0; j < func.NDim(); ++j) std::cout << x[j] << " , ";
1309  std::cout << "] ";
1310  if (fitOpt.fIntegral) {
1311  std::cout << "x2 = [ ";
1312  for (unsigned int j=0; j < func.NDim(); ++j) std::cout << data.BinUpEdge(i)[j] << " , ";
1313  std::cout << "] ";
1314  }
1315  std::cout << " y = " << y << " fval = " << fval << std::endl;
1316  }
1317 #endif
1318
1319
1320  // EvalLog protects against 0 values of fval but don't want to add in the -log sum
1321  // negative values of fval
1322  fval = std::max(fval, 0.0);
1323
1324
1325  double tmp = 0;
1326  if (useW2) {
1327  // apply weight correction . Effective weight is error^2/ y
1328  // and expected events in bins is fval/weight
1329  // can apply correction only when y is not zero otherwise weight is undefined
1330  // (in case of weighted likelihood I don't care about the constant term due to
1331  // the saturated model)
1332  if (y != 0) {
1333  double error = data.Error(i);
1334  double weight = (error*error)/y; // this is the bin effective weight
1335  if (extended) {
1336  tmp = fval * weight;
1337  // wTot += weight;
1338  // w2Tot += weight*weight;
1339  }
1340  tmp -= weight * y * ROOT::Math::Util::EvalLog( fval);
1341  }
1342
1343  // need to compute total weight and weight-square
1344  // if (extended ) {
1345  // nuTot += fval;
1346  // }
1347
1348  }
1349  else {
1350  // standard case no weights or iWeight=1
1351  // this is needed for Poisson likelihood (which are extened and not for multinomial)
1352  // the formula below include constant term due to likelihood of saturated model (f(x) = y)
1353  // (same formula as in Baker-Cousins paper, page 439 except a factor of 2
1354  if (extended) tmp = fval -y ;
1355  if (y > 0) {
1356  tmp += y * (ROOT::Math::Util::EvalLog( y) - ROOT::Math::Util::EvalLog(fval));
1357  nPoints++;
1358  }
1359  }
1360
1361
1362  nloglike += tmp;
1363  }
1364
1365  // if (notExtended) {
1366  // // not extended : remove from the Likelihood the global Poisson term
1367  // if (!useW2)
1368  // nloglike -= nuTot - yTot * ROOT::Math::Util::EvalLog( nuTot);
1369  // else {
1370  // // this needs to be checked
1371  // nloglike -= wTot* nuTot - wTot* yTot * ROOT::Math::Util::EvalLog( nuTot);
1372  // }
1373
1374  // }
1375
1376  // if (extended && useW2) {
1377  // // effective total weight is total sum of weight square / sum of weights
1378  // //nloglike += (w2Tot/wTot) * nuTot;
1379  // }
1380
1381
1382 #ifdef DEBUG
1383  std::cout << "Loglikelihood = " << nloglike << std::endl;
1384 #endif
1385
1386  return nloglike;
1387 }
1388
1389 void FitUtil::EvaluatePoissonLogLGradient(const IModelFunction & f, const BinData & data, const double * p, double * grad ) {
1390  // evaluate the gradient of the Poisson log likelihood function
1391
1393  assert (fg != 0); // must be called by a grad function
1394  const IGradModelFunction & func = *fg;
1395
1396  unsigned int n = data.Size();
1397
1398  const DataOptions & fitOpt = data.Opt();
1399  bool useBinIntegral = fitOpt.fIntegral && data.HasBinEdges();
1400  bool useBinVolume = (fitOpt.fBinVolume && data.HasBinEdges());
1401
1402  double wrefVolume = 1.0;
1403  std::vector<double> xc;
1404  if (useBinVolume) {
1405  if (fitOpt.fNormBinVolume) wrefVolume /= data.RefVolume();
1406  xc.resize(data.NDim() );
1407  }
1408
1409  IntegralEvaluator<> igEval( func, p, useBinIntegral);
1410
1411  unsigned int npar = func.NPar();
1413  std::vector<double> g( npar);
1414
1415  for (unsigned int i = 0; i < n; ++ i) {
1416  const double * x1 = data.Coords(i);
1417  double y = data.Value(i);
1418  double fval = 0;
1419  const double * x2 = 0;
1420
1421  double binVolume = 1.0;
1422  if (useBinVolume) {
1423  x2 = data.BinUpEdge(i);
1424  unsigned int ndim = data.NDim();
1425  for (unsigned int j = 0; j < ndim; ++j) {
1426  binVolume *= std::abs( x2[j]-x1[j] );
1427  xc[j] = 0.5*(x2[j]+ x1[j]);
1428  }
1429  // normalize the bin volume using a reference value
1430  binVolume *= wrefVolume;
1431  }
1432
1433  const double * x = (useBinVolume) ? &xc.front() : x1;
1434
1435  if (!useBinIntegral) {
1436  fval = func ( x, p );
1438  }
1439  else {
1440  // calculate integral (normalized by bin volume)
1441  // need to set function and parameters here in case loop is parallelized
1442  x2 = data.BinUpEdge(i);
1443  fval = igEval( x1, x2) ;
1445  }
1446  if (useBinVolume) fval *= binVolume;
1447
1449  for (unsigned int kpar = 0; kpar < npar; ++ kpar) {
1450
1451  // correct gradient for bin volumes
1452  if (useBinVolume) gradFunc[kpar] *= binVolume;
1453
1454  // df/dp * (1. - y/f )
1455  if (fval > 0)
1456  g[kpar] += gradFunc[ kpar ] * ( 1. - y/fval );
1457  else if (gradFunc [ kpar] != 0) {
1458  const double kdmax1 = std::sqrt( std::numeric_limits<double>::max() );
1459  const double kdmax2 = std::numeric_limits<double>::max() / (4*n);
1460  double gg = kdmax1 * gradFunc[ kpar ];
1461  if ( gg > 0) gg = std::min( gg, kdmax2);
1462  else gg = std::max(gg, - kdmax2);
1463  g[kpar] -= gg;
1464  }
1465  }
1466
1467  // copy result
1469  }
1470 }
1471
1472 }
1473
1474 } // end namespace ROOT
1475
bool CheckValue(double &rval)
Definition: FitUtil.cxx:293
float xmin
Definition: THbookFile.cxx:93
Interface (abstract class) for generic functions objects of one-dimension Provides a method to evalua...
Definition: IFunction.h:134
bool HaveAsymErrors() const
flag to control if data provides asymmetric errors on the value
Definition: BinData.h:143
unsigned int NDim() const
return coordinate data dimension
Definition: UnBinData.h:280
bool HasBinEdges() const
query if the data store the bin edges instead of the center
Definition: BinData.h:505
void EvaluateChi2Gradient(const IModelFunction &func, const BinData &data, const double *x, double *grad, unsigned int &nPoints)
evaluate the Chi2 gradient given a model function and the data at the point x.
Definition: FitUtil.cxx:699
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
Evaluate the all the derivatives (gradient vector) of the function with respect to the parameters at ...
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:1105
void CalculateGradientIntegral(const GFunc &gfunc, const double *x1, const double *x2, const double *p, double *g)
Definition: FitUtil.cxx:313
ErrorType GetErrorType() const
retrieve the errortype
Definition: BinData.h:522
TH1 * h
Definition: legend2.C:5
Class describing the unbinned data sets (just x coordinates values) of any dimensions.
Definition: UnBinData.h:42
Interface (abstract class) for parametric gradient multi-dimensional functions providing in addition ...
TRObject operator()(const T1 &t1) const
const double * Coords(unsigned int ipoint) const
return a pointer to the coordinates data for the given fit point
Definition: FitData.h:245
const double * GetPoint(unsigned int ipoint, double &value) const
retrieve at the same time a pointer to the coordinate data and the fit value More efficient than call...
Definition: BinData.h:347
double sqrt(double)
static const double x2[5]
void EvaluatePoissonLogLGradient(const IModelFunction &func, const BinData &data, const double *x, double *grad)
evaluate the Poisson LogL given a model function and the data at the point x.
Definition: FitUtil.cxx:1389
void SetParameters(TFitEditor::FuncParams_t &pars, TF1 *func)
Restore the parameters from pars into the function.
Definition: TFitEditor.cxx:287
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
OneDimMultiFunctionAdapter class to wrap a multidimensional function in one dimensional one...
double RefVolume() const
retrieve the reference volume used to normalize the data when the option bin volume is set ...
Definition: BinData.h:512
double Derivative1(double x)
Returns the first derivative of the function at point x, computed by Richardson&#39;s extrapolation metho...
void EvaluateLogLGradient(const IModelFunction &func, const UnBinData &data, const double *x, double *grad, unsigned int &nPoints)
evaluate the LogL gradient given a model function and the data at the point x.
Definition: FitUtil.cxx:1063
#define F1(x, y, z)
Definition: TMD5.cxx:267
unsigned int Size() const
return number of fit points
Definition: FitData.h:302
double Error(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:231
#define MATH_ERROR_MSG(loc, str)
Definition: Error.h:50
Documentation for the abstract class IBaseFunctionMultiDim.
Definition: IFunction.h:62
virtual unsigned int NPar() const =0
Return the number of Parameters.
double Weight(unsigned int ipoint) const
return weight
Definition: UnBinData.h:257
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:473
const double * BinUpEdge(unsigned int ipoint) const
return an array containing the upper edge of the bin for coordinate i In case of empty bin they could...
Definition: BinData.h:485
double EvaluateLogL(const IModelFunction &func, const UnBinData &data, const double *x, int iWeight, bool extended, unsigned int &nPoints)
evaluate the LogL given a model function and the data at the point x.
Definition: FitUtil.cxx:892
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)
Parallel evaluate the Chi2 given a model function and the data at the point x.
Definition: FitUtil.cxx:594
User Class for performing numerical integration of a function in one dimension.
Definition: Integrator.h:94
const DataOptions & Opt() const
Definition: FitData.h:318
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
double EvaluatePoissonLogL(const IModelFunction &func, const BinData &data, const double *x, int iWeight, bool extended, unsigned int &nPoints)
evaluate the Poisson LogL given a model function and the data at the point x.
Definition: FitUtil.cxx:1208
bool HaveCoordErrors() const
flag to control if data provides error on the coordinates
Definition: BinData.h:131
constexpr Double_t E()
Definition: TMath.h:74
ROOT::Math::IParamMultiFunction IModelFunction
Definition: FitUtil.h:35
double EvaluateChi2(const IModelFunction &func, const BinData &data, const double *x, unsigned int &nPoints)
Chi2 Functions.
Definition: FitUtil.cxx:337
double CorrectValue(double rval)
Definition: FitUtil.cxx:282
static const double x1[5]
double f(double x)
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 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:844
Double_t y[n]
Definition: legend1.C:17
double func(double *x, double *p)
Definition: stressTF1.cxx:213
Double_t ey[n]
Definition: legend1.C:17
double Value(unsigned int ipoint) const
return the value for the given fit point
Definition: BinData.h:217
Binding & operator=(OUT(*fun)(void))
double EvalLog(double x)
safe evaluation of log(x) with a protections against negative or zero argument to the log smooth line...
Definition: Util.h:52
User class for performing multidimensional integration.
double f2(const double *x)
const double * GetPointError(unsigned int ipoint, double &errvalue) const
Retrieve the errors on the point (coordinate and value) for the given fit point It must be called onl...
Definition: BinData.h:428
TF1 * f1
Definition: legend1.C:11
unsigned int NDim() const
return coordinate data dimension
Definition: FitData.h:310
const DataRange & Range() const