Logo ROOT  
Reference Guide
 
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
Loading...
Searching...
No Matches
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
14
16
17#include "Minuit2/MnPrint.h"
18
19// #ifndef ROOT_Math_Util
20// #include "Math/Util.h"
21// #endif
22
23#include <cmath>
24#include <cassert>
25#include <vector>
26
27namespace ROOT {
28
29namespace Minuit2 {
30
31/**
32
33
34template wrapped class for adapting to FumiliFCNBase signature
35
36@author Lorenzo Moneta
37
38@ingroup Minuit
39
40*/
41
42template <class Function>
44
45public:
46 // typedef ROOT::Math::FitMethodFunction Function;
47 typedef typename Function::Type_t Type_t;
48
49 FumiliFCNAdapter(const Function &f, unsigned int ndim, double up = 1.) : FumiliFCNBase(ndim), fFunc(f), fUp(up) {}
50
51 double operator()(std::vector<double> const &v) const override { return fFunc.operator()(&v[0]); }
52 double operator()(const double *v) const { return fFunc.operator()(v); }
53 double Up() const override { return fUp; }
54
55 void SetErrorDef(double up) override { fUp = up; }
56
57
58 /**
59 evaluate gradient hessian and function value needed by Fumili
60 */
61 void EvaluateAll(std::vector<double> const &v) override;
62
63private:
65 double fUp;
66};
67
68template <class Function>
69void FumiliFCNAdapter<Function>::EvaluateAll(std::vector<double> const &v)
70{
71 MnPrint print("FumiliFCNAdapter");
72
73 // evaluate all elements
74 unsigned int npar = Dimension();
75 if (npar != v.size())
76 print.Error("npar", npar, "v.size()", v.size());
77 assert(npar == v.size());
78 // must distinguish case of likelihood or LS
79
80 std::vector<double> &grad = Gradient();
81 std::vector<double> &hess = Hessian();
82 // reset
83 assert(grad.size() == npar);
84 grad.assign(npar, 0.0);
85 hess.assign(hess.size(), 0.0);
86
87 unsigned int ndata = fFunc.NPoints();
88
89 std::vector<double> gf(npar);
90 std::vector<double> h(hess.size());
91
92 // loop on the data points
93
94 // if FCN is of type least-square
95 if (fFunc.Type() == Function::kLeastSquare) {
96 print.Debug("Chi2 FCN: Evaluate gradient and Hessian");
97
98 for (unsigned int i = 0; i < ndata; ++i) {
99 // calculate data element and gradient (no need to compute Hessian)
100 double fval = fFunc.DataElement(&v.front(), i, &gf[0]);
101
102 for (unsigned int j = 0; j < npar; ++j) {
103 grad[j] += 2. * fval * gf[j];
104 for (unsigned int k = j; k < npar; ++k) {
105 int idx = j + k * (k + 1) / 2;
106 hess[idx] += 2.0 * gf[j] * gf[k];
107 }
108 }
109 }
110 } else if (fFunc.Type() == Function::kLogLikelihood) {
111 print.Debug("LogLikelihood FCN: Evaluate gradient and Hessian");
112 for (unsigned int i = 0; i < ndata; ++i) {
113
114 // calculate data element and gradient: returns derivative of log(pdf)
115 fFunc.DataElement(&v.front(), i, &gf[0]);
116
117 for (unsigned int j = 0; j < npar; ++j) {
118 double gfj = gf[j];
119 grad[j] -= gfj; // need a minus sign since is a NLL
120 for (unsigned int k = j; k < npar; ++k) {
121 int idx = j + k * (k + 1) / 2;
122 hess[idx] += gfj * gf[k];
123 }
124 }
125 }
126 } else if (fFunc.Type() == Function::kPoissonLikelihood) {
127 print.Debug("Poisson Likelihood FCN: Evaluate gradient and Hessian");
128 // for Poisson need Hessian computed in DataElement since one needs the bin expected value ad bin observed value
129 for (unsigned int i = 0; i < ndata; ++i) {
130 // calculate data element and gradient
131 fFunc.DataElement(&v.front(), i, gf.data(), h.data());
132 for (size_t j = 0; j < npar; ++j) {
133 grad[j] += gf[j];
134 for (unsigned int k = j; k < npar; ++k) {
135 int idx = j + k * (k + 1) / 2;
136 hess[idx] += h[idx];
137 }
138 }
139 }
140 } else {
141 print.Error("Type of fit method is not supported, it must be chi2 or log-likelihood or Poisson Likelihood");
142 }
143}
144
145} // end namespace Minuit2
146
147} // end namespace ROOT
148
149#endif // ROOT_Minuit2_FCNAdapter
#define f(i)
Definition RSha256.hxx:104
#define h(i)
Definition RSha256.hxx:106
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
Double_t(* Function)(Double_t)
Definition Functor.C:4
template wrapped class for adapting to FumiliFCNBase signature
double Up() const override
Error definition of the function.
double operator()(const double *v) const
void SetErrorDef(double up) override
add interface to set dynamically a new error definition Re-implement this function if needed.
void EvaluateAll(std::vector< double > const &v) override
evaluate gradient hessian and function value needed by Fumili
FumiliFCNAdapter(const Function &f, unsigned int ndim, double up=1.)
double operator()(std::vector< double > const &v) const override
The meaning of the vector of parameters is of course defined by the user, who uses the values of thos...
Extension of the FCNBase for the Fumili method.
void Debug(const Ts &... args)
Definition MnPrint.h:135
void Error(const Ts &... args)
Definition MnPrint.h:117
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...