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