ROOT  6.06/09
Reference Guide
FumiliFCNAdapter.h
Go to the documentation of this file.
1 // @(#)root/minuit2:$Id$
2 // Author: L. Moneta 10/2006
3 
4 /**********************************************************************
5  * *
6  * Copyright (c) 2006 ROOT Foundation, CERN/PH-SFT *
7  * *
8  **********************************************************************/
9 
10 #ifndef ROOT_Minuit2_FumiliFCNAdapter
11 #define ROOT_Minuit2_FumiliFCNAdapter
12 
13 #ifndef ROOT_Minuit2_FumiliFCNBase
14 #include "Minuit2/FumiliFCNBase.h"
15 #endif
16 
17 #ifndef ROOT_Math_FitMethodFunction
18 #include "Math/FitMethodFunction.h"
19 #endif
20 
21 #ifndef ROOT_Minuit2_MnPrint
22 #include "Minuit2/MnPrint.h"
23 #endif
24 
25 // #ifndef ROOT_Math_Util
26 // #include "Math/Util.h"
27 // #endif
28 
29 #include <cmath>
30 
31 namespace ROOT {
32 
33  namespace Minuit2 {
34 
35 /**
36 
37 
38 template wrapped class for adapting to FumiliFCNBase signature
39 
40 @author Lorenzo Moneta
41 
42 @ingroup Minuit
43 
44 */
45 
46 template< class Function>
48 
49 public:
50 
51 // typedef ROOT::Math::FitMethodFunction Function;
52  typedef typename Function::Type_t Type_t;
53 
54  FumiliFCNAdapter(const Function & f, unsigned int ndim, double up = 1.) :
55  FumiliFCNBase( ndim ),
56  fFunc(f) ,
57  fUp (up)
58  {}
59 
61 
62 
63  double operator()(const std::vector<double>& v) const {
64  return fFunc.operator()(&v[0]);
65  }
66  double operator()(const double * v) const {
67  return fFunc.operator()(v);
68  }
69  double Up() const {return fUp;}
70 
71  void SetErrorDef(double up) { fUp = up; }
72 
73  //virtual std::vector<double> Gradient(const std::vector<double>&) const;
74 
75  // forward interface
76  //virtual double operator()(int npar, double* params,int iflag = 4) const;
77 
78 
79  /**
80  evaluate gradient hessian and function value needed by fumili
81  */
82  void EvaluateAll( const std::vector<double> & v);
83 
84 private:
85 
86 
87 //data member
88 
89  const Function & fFunc;
90  double fUp;
91 };
92 
93 
94 template<class Function>
95 void FumiliFCNAdapter<Function>::EvaluateAll( const std::vector<double> & v) {
96 
97  //typedef FumiliFCNAdapter::Function Function;
98 
99  //evaluate all elements
100  unsigned int npar = Dimension();
101  if (npar != v.size() ) std::cout << "npar = " << npar << " " << v.size() << std::endl;
102  assert(npar == v.size());
103  //must distinguish case of likelihood or LS
104 
105  std::vector<double> & grad = Gradient();
106  std::vector<double> & hess = Hessian();
107  // reset
108  assert(grad.size() == npar);
109  grad.assign( npar, 0.0);
110  hess.assign( hess.size(), 0.0);
111 
112  double sum = 0;
113  unsigned int ndata = fFunc.NPoints();
114 
115  std::vector<double> gf(npar);
116 
117  //loop on the data points
118 
119 
120  // assume for now least-square
121  if (fFunc.Type() == Function::kLeastSquare) {
122 
123  for (unsigned int i = 0; i < ndata; ++i) {
124  // calculate data element and gradient
125  double fval = fFunc.DataElement(&v.front(), i, &gf[0]);
126 
127  // t.b.d should protect for bad values of fval
128  sum += fval*fval;
129 
130  for (unsigned int j = 0; j < npar; ++j) {
131  grad[j] += 2. * fval * gf[j];
132  for (unsigned int k = j; k < npar; ++ k) {
133  int idx = j + k*(k+1)/2;
134  hess[idx] += 2.0 * gf[j] * gf[k];
135  }
136  }
137  }
138  }
139  else if (fFunc.Type() == Function::kLogLikelihood) {
140 
141 
142  for (unsigned int i = 0; i < ndata; ++i) {
143 
144  // calculate data element and gradient
145  // return value is log of pdf and derivative of the log(Pdf)
146  double fval = fFunc.DataElement(&v.front(), i, &gf[0]);
147 
148  sum -= fval;
149 
150  for (unsigned int j = 0; j < npar; ++j) {
151  double gfj = gf[j] ;
152  grad[j] -= gfj;
153  for (unsigned int k = j; k < npar; ++ k) {
154  int idx = j + k*(k+1)/2;
155  hess[idx] += gfj * gf[k] ;
156  }
157  }
158  }
159  }
160  else {
161  MN_ERROR_MSG("FumiliFCNAdapter: type of fit method is not supported, it must be chi2 or log-likelihood");
162  }
163 }
164 
165 
166  } // end namespace Minuit2
167 
168 } // end namespace ROOT
169 
170 
171 
172 
173 #endif //ROOT_Minuit2_FCNAdapter
const int ndata
#define MN_ERROR_MSG(str)
Definition: MnPrint.h:113
Namespace for new ROOT classes and functions.
Definition: ROOT.py:1
#define assert(cond)
Definition: unittest.h:542
double operator()(const double *v) const
double operator()(const std::vector< double > &v) const
The meaning of the vector of parameters is of course defined by the user, who uses the values of thos...
void SetErrorDef(double up)
add interface to set dynamically a new error definition Re-implement this function if needed...
SVector< double, 2 > v
Definition: Dict.h:5
double Up() const
Error definition of the function.
double f(double x)
template wrapped class for adapting to FumiliFCNBase signature
Extension of the FCNBase for the Fumili method.
Definition: FumiliFCNBase.h:47
FumiliFCNAdapter(const Function &f, unsigned int ndim, double up=1.)
void EvaluateAll(const std::vector< double > &v)
evaluate gradient hessian and function value needed by fumili