Logo ROOT  
Reference Guide
WrappedMultiTF1.h
Go to the documentation of this file.
1// @(#)root/mathmore:$Id$
2// Author: L. Moneta Wed Sep 6 09:52:26 2006
3
4/**********************************************************************
5 * *
6 * Copyright (c) 2006 LCG ROOT Math Team, CERN/PH-SFT *
7 * *
8 * *
9 **********************************************************************/
10
11// Header file for class WrappedTFunction
12
13#ifndef ROOT_Math_WrappedMultiTF1
14#define ROOT_Math_WrappedMultiTF1
15
16
17#include "Math/IParamFunction.h"
18
19#include "TF1.h"
20#include <string>
21#include <vector>
22#include <algorithm>
23
24namespace ROOT {
25
26 namespace Math {
27
28 namespace Internal {
29 double DerivPrecision(double eps);
30 TF1 *CopyTF1Ptr(const TF1 *funcToCopy);
31 };
32
33 /**
34 Class to Wrap a ROOT Function class (like TF1) in a IParamMultiFunction interface
35 of multi-dimensions to be used in the ROOT::Math numerical algorithm.
36 This wrapper class does not own the TF1 pointer, so it assumes it exists during the wrapper lifetime.
37 The class copy the TF1 pointer only when it owns it.
38
39 The class from ROOT version 6.03 does not contain anymore a copy of the parameters. The parameters are
40 stored in the TF1 class.
41
42 @ingroup CppFunctions
43 */
44
45 //LM note: are there any issues when cloning the class for the parameters that are not copied anymore ??
46
47 template<class T>
49
50 public:
51
54
55 /**
56 constructor from a function pointer to a TF1
57 If dim = 0 dimension is taken from TF1::GetNdim().
58 In case of multi-dimensional function created using directly TF1 object the dimension
59 returned by TF1::GetNdim is always 1. The user must then pass the correct value of dim
60 */
61 WrappedMultiTF1Templ(TF1 &f, unsigned int dim = 0);
62
63 /**
64 Destructor (no operations). Function pointer is not owned
65 */
67 {
68 if (fOwnFunc && fFunc) delete fFunc;
69 }
70
71 /**
72 Copy constructor
73 */
75
76 /**
77 Assignment operator
78 */
80
81 /** @name interface inherited from IParamFunction */
82
83 /**
84 Clone the wrapper but not the original function
85 */
87 {
88 return new WrappedMultiTF1Templ<T>(*this);
89 }
90
91 /**
92 Retrieve the dimension of the function
93 */
94 unsigned int NDim() const override
95 {
96 return fDim;
97 }
98
99 /// get the parameter values (return values from TF1)
100 const double *Parameters() const override
101 {
102 //return (fParams.size() > 0) ? &fParams.front() : 0;
103 return fFunc->GetParameters();
104 }
105
106 /// set parameter values (only the cached one in this class,leave unchanges those of TF1)
107 void SetParameters(const double *p) override
108 {
109 //std::copy(p,p+fParams.size(),fParams.begin());
111 }
112
113 /// return number of parameters
114 unsigned int NPar() const override
115 {
116 // return fParams.size();
117 return fFunc->GetNpar();
118 }
119
120 /// return parameter name (from TF1)
121 std::string ParameterName(unsigned int i) const override {
122 return std::string(fFunc->GetParName(i));
123 }
124
125 // evaluate the derivative of the function with respect to the parameters
126 void ParameterGradient(const T *x, const double *par, T *grad) const override;
127
128 // evaluate the hessian of the function with respect to the parameters
129 bool ParameterHessian(const T *x, const double *par, T *h) const override;
130
131 // return capability of computing parameter Hessian
132 bool HasParameterHessian() const override;
133
134 // evaluate the 2nd derivatives of the function with respect to the parameters
135 bool ParameterG2(const T *, const double *, T *) const override {
136 return false; // not yet implemented
137 }
138
139 /// precision value used for calculating the derivative step-size
140 /// h = eps * |x|. The default is 0.001, give a smaller in case function changes rapidly
141 static void SetDerivPrecision(double eps);
142
143 /// get precision value used for calculating the derivative step-size
144 static double GetDerivPrecision();
145
146 /// method to retrieve the internal function pointer
147 const TF1 *GetFunction() const
148 {
149 return fFunc;
150 }
151
152 /// method to set a new function pointer and copy it inside.
153 /// By calling this method the class manages now the passed TF1 pointer
154 void SetAndCopyFunction(const TF1 *f = nullptr);
155
156 private:
157 /// evaluate function passing coordinates x and vector of parameters
158 T DoEvalPar(const T *x, const double *p) const override
159 {
160 return fFunc->EvalPar(x, p);
161 }
162
163 /// evaluate function using the cached parameter values (of TF1)
164 /// re-implement for better efficiency
165 T DoEvalVec(const T *x) const
166 {
167 return fFunc->EvalPar(x, 0);
168 }
169
170 /// evaluate function using the cached parameter values (of TF1)
171 /// re-implement for better efficiency
172 T DoEval(const T *x) const override
173 {
174 // no need to call InitArg for interpreted functions (done in ctor)
175
176 //const double * p = (fParams.size() > 0) ? &fParams.front() : 0;
177
178 return fFunc->EvalPar(x, nullptr);
179 }
180
181 /// evaluate the partial derivative with respect to the parameter
182 T DoParameterDerivative(const T *x, const double *p, unsigned int ipar) const override;
183
184 bool fLinear; // flag for linear functions
185 bool fPolynomial; // flag for polynomial functions
186 bool fOwnFunc; // flag to indicate we own the TF1 function pointer
187 TF1 *fFunc; // pointer to ROOT function
188 unsigned int fDim; // cached value of dimension
189 //std::vector<double> fParams; // cached vector with parameter values
190
191 };
192
193 /**
194 * Auxiliar class to bypass the (provisional) lack of vectorization in TFormula::EvalPar.
195 *
196 * WrappedMultiTF1Templ::DoParameterDerivation calls TFormula::EvalPar in the case of a general linear function
197 * built with TFormula using ++; as EvalPar is not vectorized, in order to generalize DoParameterDerivative with
198 * a general type T, we use this auxiliar class to branch the code in compile time with the double
199 * specialization (that can call EvalPar) and the general implementation (that throws an error in the case of
200 * general linear function).
201 */
202 template <class T>
204 static T DoParameterDerivative(const WrappedMultiTF1Templ<T> *, const T *, unsigned int)
205 {
206 Error("DoParameterDerivative", "The vectorized implementation of DoParameterDerivative does not support"
207 "general linear functions built in TFormula with ++");
208
209 return TMath::SignalingNaN();
210 }
211 };
212
213 template <>
215 static double
216 DoParameterDerivative(const WrappedMultiTF1Templ<double> *wrappedFunc, const double *x, unsigned int ipar)
217 {
218 const TFormula *df = dynamic_cast<const TFormula *>(wrappedFunc->GetFunction()->GetLinearPart(ipar));
219 assert(df != nullptr);
220 return (const_cast<TFormula *>(df))->EvalPar(x); // derivatives should not depend on parameters since
221 // function is linear
222 }
223 };
224
225 // implementations for WrappedMultiTF1Templ<T>
226 template<class T>
228 fLinear(false),
229 fPolynomial(false),
230 fOwnFunc(false),
231 fFunc(&f),
232 fDim(dim)
233 //fParams(f.GetParameters(),f.GetParameters()+f.GetNpar())
234 {
235 // constructor of WrappedMultiTF1Templ<T>
236 // pass a dimension if dimension specified in TF1 does not correspond to real dimension
237 // for example in case of multi-dimensional TF1 objects defined as TF1 (i.e. for functions with dims > 3 )
238 if (fDim == 0) fDim = fFunc->GetNdim();
239
240 // check that in case function is linear the linear terms are not zero
241 // function is linear when is a TFormula created with "++"
242 // hyperplane are not yet existing in TFormula
243 if (fFunc->IsLinear()) {
244 int ip = 0;
245 fLinear = true;
246 while (fLinear && ip < fFunc->GetNpar()) {
247 fLinear &= (fFunc->GetLinearPart(ip) != nullptr);
248 ip++;
249 }
250 }
251 // distinguish case of polynomial functions and linear functions
252 if (fDim == 1 && fFunc->GetNumber() >= 300 && fFunc->GetNumber() < 310) {
253 fLinear = true;
254 fPolynomial = true;
255 }
256 }
257
258 template<class T>
260 BaseParamFunc(rhs),
261 fLinear(rhs.fLinear),
262 fPolynomial(rhs.fPolynomial),
263 fOwnFunc(rhs.fOwnFunc),
264 fFunc(rhs.fFunc),
265 fDim(rhs.fDim)
266 //fParams(rhs.fParams)
267 {
268 // copy constructor
270 }
271
272 template<class T>
274 {
275 // Assignment operator
276 if (this == &rhs) return *this; // time saving self-test
277 fLinear = rhs.fLinear;
278 fPolynomial = rhs.fPolynomial;
279 fOwnFunc = rhs.fOwnFunc;
280 fDim = rhs.fDim;
281 //fParams = rhs.fParams;
282 return *this;
283 }
284
285 template <class T>
286 void WrappedMultiTF1Templ<T>::ParameterGradient(const T *x, const double *par, T *grad) const
287 {
288 // evaluate the gradient of the function with respect to the parameters
289 //IMPORTANT NOTE: TF1::GradientPar returns 0 for fixed parameters to avoid computing useless derivatives
290 // BUT the TLinearFitter wants to have the derivatives also for fixed parameters.
291 // so in case of fLinear (or fPolynomial) a non-zero value will be returned for fixed parameters
292
293 if (!fLinear) {
294 // need to set parameter values
295 fFunc->SetParameters(par);
296 // no need to call InitArgs (it is called in TF1::GradientPar)
297 double prec = this->GetDerivPrecision();
298 fFunc->GradientPar(x, grad, prec);
299 } else { // case of linear functions
300 unsigned int np = NPar();
301 for (unsigned int i = 0; i < np; ++i)
302 grad[i] = DoParameterDerivative(x, par, i);
303 }
304 }
305
306 // struct for dealing of generic Hessian computation, since it is available only in TFormula
307 template <class T>
309 static bool Hessian(TF1 *, const T *, const double *, T *)
310 {
311 Error("Hessian", "The vectorized implementation of ParameterHessian is not supported");
312 return false;
313 }
314 static bool IsAvailable(TF1 * ) { return false;}
315 };
316
317 template <>
319 static bool Hessian(TF1 * func, const double *x, const double * par, double * h)
320 {
321 // compute Hessian if TF1 is a formula based
322 unsigned int np = func->GetNpar();
323 auto formula = func->GetFormula();
324 if (!formula) return false;
325 std::vector<double> h2(np*np);
326 func->SetParameters(par);
327 formula->HessianPar(x,h2);
328 for (unsigned int i = 0; i < np; i++) {
329 for (unsigned int j = 0; j <= i; j++) {
330 unsigned int ih = j + i *(i+1)/2; // formula for j <= i
331 unsigned int im = i*np + j;
332 h[ih] = h2[im];
333 }
334 }
335 return true;
336 }
337 static bool IsAvailable(TF1 * func) {
338 auto formula = func->GetFormula();
339 if (!formula) return false;
340 return formula->GenerateHessianPar();
341 }
342 };
343
344
345
346 template <class T>
347 bool WrappedMultiTF1Templ<T>::ParameterHessian(const T *x, const double *par, T *h) const
348 {
349 if (fLinear) {
350 std::fill(h, h + NPar()*(NPar()+1)/2, 0.0);
351 return true;
352 }
353 return GeneralHessianCalc<T>::Hessian(fFunc, x, par, h);
354 }
355
356 template <class T>
358 {
360 }
361
362 template <class T>
363 T WrappedMultiTF1Templ<T>::DoParameterDerivative(const T *x, const double *p, unsigned int ipar) const
364 {
365 // evaluate the derivative of the function with respect to parameter ipar
366 // see note above concerning the fixed parameters
367 if (!fLinear) {
368 fFunc->SetParameters(p);
369 double prec = this->GetDerivPrecision();
370 return fFunc->GradientPar(ipar, x, prec);
371 }
372 if (fPolynomial) {
373 // case of polynomial function (no parameter dependency) (case for dim = 1)
374 assert(fDim == 1);
375 if (ipar == 0) return 1.0;
376#ifdef R__HAS_VECCORE
377 return vecCore::math::Pow(x[0], static_cast<T>(ipar));
378#else
379 return std::pow(x[0], static_cast<int>(ipar));
380#endif
381 } else {
382 // case of general linear function (built in TFormula with ++ )
384 }
385 }
386 template<class T>
388 {
390 }
391
392 template<class T>
394 {
396 }
397
398 template<class T>
400 {
401 const TF1 *funcToCopy = (f) ? f : fFunc;
402 fFunc = ::ROOT::Math::Internal::CopyTF1Ptr(funcToCopy);
403 fOwnFunc = true;
404 }
405
407
408 } // end namespace Math
409
410} // end namespace ROOT
411
412
413#endif /* ROOT_Fit_WrappedMultiTF1 */
#define f(i)
Definition: RSha256.hxx:104
#define h(i)
Definition: RSha256.hxx:106
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition: TError.cxx:188
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t np
Documentation for the abstract class IBaseFunctionMultiDim.
Definition: IFunction.h:62
Interface (abstract class) for parametric gradient multi-dimensional functions providing in addition ...
Class to Wrap a ROOT Function class (like TF1) in a IParamMultiFunction interface of multi-dimensions...
void SetAndCopyFunction(const TF1 *f=nullptr)
method to set a new function pointer and copy it inside.
T DoParameterDerivative(const T *x, const double *p, unsigned int ipar) const override
evaluate the partial derivative with respect to the parameter
bool HasParameterHessian() const override
~WrappedMultiTF1Templ() override
Destructor (no operations).
const double * Parameters() const override
get the parameter values (return values from TF1)
T DoEval(const T *x) const override
evaluate function using the cached parameter values (of TF1) re-implement for better efficiency
bool ParameterG2(const T *, const double *, T *) const override
Evaluate all the second derivatives (diagonal ones) of the function with respect to the parameters at...
static void SetDerivPrecision(double eps)
precision value used for calculating the derivative step-size h = eps * |x|.
static double GetDerivPrecision()
get precision value used for calculating the derivative step-size
ROOT::Math::IParametricGradFunctionMultiDimTempl< T > BaseParamFunc
WrappedMultiTF1Templ & operator=(const WrappedMultiTF1Templ< T > &rhs)
Assignment operator.
std::string ParameterName(unsigned int i) const override
return parameter name (from TF1)
const TF1 * GetFunction() const
method to retrieve the internal function pointer
unsigned int NDim() const override
Retrieve the dimension of the function.
unsigned int NPar() const override
return number of parameters
IMultiGenFunctionTempl< T > * Clone() const override
Clone the wrapper but not the original function.
void SetParameters(const double *p) override
set parameter values (only the cached one in this class,leave unchanges those of TF1)
void ParameterGradient(const T *x, const double *par, T *grad) const override
Evaluate the all the derivatives (gradient vector) of the function with respect to the parameters at ...
ROOT::Math::IParametricFunctionMultiDimTempl< T >::BaseFunc BaseFunc
T DoEvalVec(const T *x) const
evaluate function using the cached parameter values (of TF1) re-implement for better efficiency
T DoEvalPar(const T *x, const double *p) const override
evaluate function passing coordinates x and vector of parameters
WrappedMultiTF1Templ(TF1 &f, unsigned int dim=0)
constructor from a function pointer to a TF1 If dim = 0 dimension is taken from TF1::GetNdim().
bool ParameterHessian(const T *x, const double *par, T *h) const override
Evaluate the all the Hessian (second derivatives matrix) of the function with respect to the paramete...
1-Dim function class
Definition: TF1.h:213
virtual Int_t GetNumber() const
Definition: TF1.h:503
virtual Int_t GetNpar() const
Definition: TF1.h:486
virtual Double_t * GetParameters() const
Definition: TF1.h:525
virtual TFormula * GetFormula()
Definition: TF1.h:458
virtual const char * GetParName(Int_t ipar) const
Definition: TF1.h:534
virtual Double_t EvalPar(const Double_t *x, const Double_t *params=nullptr)
Evaluate function with given coordinates and parameters.
Definition: TF1.cxx:1469
virtual Bool_t IsLinear() const
Definition: TF1.h:607
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:649
virtual Int_t GetNdim() const
Definition: TF1.h:490
virtual const TObject * GetLinearPart(Int_t i) const
Definition: TF1.h:470
The Formula class.
Definition: TFormula.h:87
bool GenerateHessianPar()
Generate hessian computation routine with respect to the parameters.
Definition: TFormula.cxx:3318
RVec< PromoteTypes< T0, T1 > > pow(const T0 &x, const RVec< T1 > &v)
Definition: RVec.hxx:1792
Double_t x[n]
Definition: legend1.C:17
Namespace for new Math classes and functions.
double T(double x)
Definition: ChebyshevPol.h:34
double DerivPrecision(double eps)
Definition: WrappedTF1.cxx:25
TF1 * CopyTF1Ptr(const TF1 *funcToCopy)
Definition: WrappedTF1.cxx:33
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.
Double_t SignalingNaN()
Returns a signaling NaN as defined by IEEE 754](http://en.wikipedia.org/wiki/NaN#Signaling_NaN).
Definition: TMath.h:908
fill
Definition: fit1_py.py:6
static bool Hessian(TF1 *func, const double *x, const double *par, double *h)
static bool Hessian(TF1 *, const T *, const double *, T *)
static double DoParameterDerivative(const WrappedMultiTF1Templ< double > *wrappedFunc, const double *x, unsigned int ipar)
Auxiliar class to bypass the (provisional) lack of vectorization in TFormula::EvalPar.
static T DoParameterDerivative(const WrappedMultiTF1Templ< T > *, const T *, unsigned int)