Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
TFumiliMinimizer.cxx
Go to the documentation of this file.
1// @(#)root/fumili:$Id$
2// Author: L. Moneta Wed Oct 25 16:28:55 2006
3
4/**********************************************************************
5 * *
6 * Copyright (c) 2006 LCG ROOT Math Team, CERN/PH-SFT *
7 * *
8 * *
9 **********************************************************************/
10
11// Implementation file for class TFumiliMinimizer
12
13#include "TFumiliMinimizer.h"
14#include "Math/IFunction.h"
15#include "Math/Util.h"
16#include "TError.h"
17
18#include "TFumili.h"
19
20#include <iostream>
21#include <cassert>
22#include <algorithm>
23#include <functional>
24
25
26// setting USE_FUMILI_FUNCTION will use the Derivatives provided by Fumili
27// instead of what provided in FitUtil::EvalChi2Residual
28// t.d.: use still standard Chi2 but replace model function
29// with a gradient function where gradient is computed by TFumili
30// since TFumili knows the step size can calculate it better
31// Derivative in Fumili are very fast (1 extra call for each parameter)
32// + 1 function evaluation
33//
34//#define USE_FUMILI_FUNCTION
35#ifdef USE_FUMILI_FUNCTION
36bool gUseFumiliFunction = true;
37//#include "FumiliFunction.h"
38// fit method function used in TFumiliMinimizer
39
42#include "Fit/Chi2FCN.h"
43#include "TF1.h"
44#include "TFumili.h"
45
46template<class MethodFunc>
48
50
51public:
53 ROOT::Math::FitMethodFunction(func->NDim(), func->NPoints() ),
54 fFumili(fumili),
55 fObjFunc(0)
56 {
57 fObjFunc = dynamic_cast<const MethodFunc *>(func);
58 assert(fObjFunc != 0);
59
60 // create TF1 class from model function
61 fModFunc = new TF1("modfunc",ROOT::Math::ParamFunctor( &fObjFunc->ModelFunction() ) );
62 fFumili->SetUserFunc(fModFunc);
63 }
64
65 ROOT::Math::FitMethodFunction::Type_t Type() const { return fObjFunc->Type(); }
66
67 FumiliFunction * Clone() const { return new FumiliFunction(fFumili, fObjFunc); }
68
69
70 // recalculate data element using Fumili stuff
71 double DataElement(const double * /*par */, unsigned int i, double * g, double *) const {
72
73 // parameter values are inside TFumili
74
75 // suppose type is bin likelihood
76 unsigned int npar = fObjFunc->NDim();
77 double y = 0;
78 double invError = 0;
79 const double *x = fObjFunc->Data().GetPoint(i,y,invError);
80 double fval = fFumili->EvalTFN(g,const_cast<double *>( x));
81 fFumili->Derivatives(g, const_cast<double *>( x));
82
85 for (unsigned int k = 0; k < npar; ++k) {
86 g[k] *= ( y/fval - 1.) ;
87 }
88
89 return logPdf;
90 }
91 else if (fObjFunc->Type() == ROOT::Math::FitMethodFunction::kLeastSquare ) {
92 double resVal = (y-fval)*invError;
93 for (unsigned int k = 0; k < npar; ++k) {
94 g[k] *= -invError;
95 }
96 return resVal;
97 }
98
99 return 0;
100 }
101
102
103private:
104
105 double DoEval(const double *x ) const {
106 return (*fObjFunc)(x);
107 }
108
109 TFumili * fFumili;
110 const MethodFunc * fObjFunc;
111 TF1 * fModFunc;
112
113};
114#else
116#endif
117//______________________________________________________________________________
118//
119// TFumiliMinimizer class implementing the ROOT::Math::Minimizer interface using
120// TFumili.
121// This class is normally instantiates using the plug-in manager
122// (plug-in with name Fumili or TFumili)
123// In addition the user can choose the minimizer algorithm: Migrad (the default one), Simplex, or Minimize (combined Migrad + Simplex)
124//
125//__________________________________________________________________________________________
126
127// initialize the static instances
128
132
133
134
135
137 fDim(0),
138 fNFree(0),
139 fMinVal(0),
140 fEdm(-1),
141 fFumili(nullptr)
142{
143 // Constructor for TFumiliMinimier class
144
145 // construct with npar = 0 (by default a value of 25 is used in TFumili for allocating the arrays)
146#ifdef USE_STATIC_TMINUIT
147 // allocate here only the first time
148 if (fgFumili == 0) fgFumili = new TFumili(0);
150#else
151 if (fFumili) delete fFumili;
152 fFumili = new TFumili(0);
154#endif
155
156}
157
158
160{
161 // Destructor implementation.
162 if (fFumili) delete fFumili;
163}
164
165
166
168 // Set the objective function to be minimized, by passing a function object implement the
169 // basic multi-dim Function interface. In this case the derivatives will be
170 // calculated by Fumili
171
172 // Here a TFumili instance is created since only at this point we know the number of parameters
173 // needed to create TFumili
174 fDim = func.NDim();
176
177 if(func.HasGradient()) {
178 // In this case the function derivatives are provided
179 // by the user via this interface and there not calculated by Fumili.
180
181 fDim = func.NDim();
183
184 // for Fumili the fit method function interface is required
186 if (!fcnfunc) {
187 Error("SetFunction","Wrong Fit method function type used for Fumili");
188 return;
189 }
190 // assign to the static pointer (NO Thread safety here)
191 fgFunc = nullptr;
194
195 return;
196 }
197
198 // for Fumili the fit method function interface is required
199 const ROOT::Math::FitMethodFunction * fcnfunc = dynamic_cast<const ROOT::Math::FitMethodFunction *>(&func);
200 if (!fcnfunc) {
201 Error("SetFunction","Wrong Fit method function type used for Fumili");
202 return;
203 }
204 // assign to the static pointer (NO Thread safety here)
206 fgGradFunc = nullptr;
208
209#ifdef USE_FUMILI_FUNCTION
210 if (gUseFumiliFunction) {
215 }
216#endif
217
218}
219
220void TFumiliMinimizer::Fcn( int & , double * g , double & f, double * x , int /* iflag */) {
221 // implementation of FCN static function used internally by TFumili.
222 // Adapt IMultiGenFunction interface to TFumili FCN static function
223 f = TFumiliMinimizer::EvaluateFCN(const_cast<double*>(x),g);
224}
225
226// void TFumiliMinimizer::FcnGrad( int &, double * g, double & f, double * x , int iflag ) {
227// // implementation of FCN static function used internally by TFumili.
228// // Adapt IMultiGradFunction interface to TFumili FCN static function in the case of user
229// // provided gradient.
230// ROOT::Math::IMultiGradFunction * gFunc = dynamic_cast<ROOT::Math::IMultiGradFunction *> ( fgFunc);
231
232// assert(gFunc != 0);
233// f = gFunc->operator()(x);
234
235// // calculates also derivatives
236// if (iflag == 2) gFunc->Gradient(x,g);
237// }
238
239double TFumiliMinimizer::EvaluateFCN(const double * x, double * grad) {
240 // function called to evaluate the FCN at the value x
241 // calculates also the matrices of the second derivatives of the objective function needed by FUMILI
242
243
244 //typedef FumiliFCNAdapter::Function Function;
245
246
247
248 // reset
249// assert(grad.size() == npar);
250// grad.assign( npar, 0.0);
251// hess.assign( hess.size(), 0.0);
252
253 double sum = 0;
254 unsigned int ndata = 0;
255 unsigned int npar = 0;
256 if (fgFunc) {
257 ndata = fgFunc->NPoints();
258 npar = fgFunc->NDim();
259 fgFunc->UpdateNCalls();
260 }
261 else if (fgGradFunc) {
262 ndata = fgGradFunc->NPoints();
263 npar = fgGradFunc->NDim();
264 fgGradFunc->UpdateNCalls();
265 }
266
267 // eventually store this matrix as static member to optimize speed
268 std::vector<double> gf(npar);
269 std::vector<double> hess(npar*(npar+1)/2);
270 std::vector<double> h(npar*(npar+1)/2);
271
272 // reset gradients
273 for (unsigned int ipar = 0; ipar < npar; ++ipar)
274 grad[ipar] = 0;
275
276
277 //loop on the data points
278//#define DEBUG
279#ifdef DEBUG
280 std::cout << "=============================================";
281 std::cout << "par = ";
282 for (unsigned int ipar = 0; ipar < npar; ++ipar)
283 std::cout << x[ipar] << "\t";
284 std::cout << std::endl;
285 if (fgFunc) std::cout << "type " << fgFunc->Type() << std::endl;
286#endif
287
288
289 // assume for now least-square
290 // since TFumili does not use errordef we must divide chi2 by 2
293
294 double fval = 0;
295 for (unsigned int i = 0; i < ndata; ++i) {
296 // calculate data element and gradient
297 // DataElement returns (f-y)/s and gf is derivatives of model function multiplied by (-1/sigma)
298 if (gUseFumiliFunction) {
299 fval = fgFunc->DataElement( x, i, &gf[0]);
300 }
301 else {
302 if (fgFunc != nullptr)
303 fval = fgFunc->DataElement(x, i, gf.data(), h.data() );
304 else
305 fval = fgGradFunc->DataElement(x, i, gf.data(), h.data());
306 }
307
308 // t.b.d should protect for bad values of fval
309 sum += 0.5 * fval * fval; // need to divide chi2 by 2
310
311 for (unsigned int j = 0; j < npar; ++j) {
312 grad[j] += fval * gf[j];
313 for (unsigned int k = j; k < npar; ++ k) {
314 int idx = j + k*(k+1)/2;
315 //hess[idx] += gf[j] * gf[k];
316 hess[idx] += 0.5 * h[idx]; // h is gradient of full residual (2 * gf[j] n* gf[k] )
317 }
318 }
319 }
320 }
323
324
325 double fval = 0;
326
327 //std::cout << "\t x " << x[0] << " " << x[1] << " " << x[2] << std::endl;
328
329 for (unsigned int i = 0; i < ndata; ++i) {
330
331 if (gUseFumiliFunction) {
332 fval = fgFunc->DataElement( x, i, &gf[0]);
333 }
334 else {
335 // calculate data element and gradient
336 if (fgFunc != nullptr)
337 fval = fgFunc->DataElement(x, i, &gf[0], h.data());
338 else
339 fval = fgGradFunc->DataElement(x, i, &gf[0], h.data());
340 }
341
342 sum += fval;
343
344 for (unsigned int j = 0; j < npar; ++j) {
345 grad[j] += gf[j]; // maybe a factor of 2 here???
346 for (unsigned int k = j; k < npar; ++ k) {
347 int idx = j + k*(k+1)/2;
348 hess[idx] += h[idx];
349 }
350 }
351 }
354
355 double fval = 0;
356
357 for (unsigned int i = 0; i < ndata; ++i) {
358
359 if (gUseFumiliFunction) {
360 fval = fgFunc->DataElement(x, i, &gf[0]);
361 } else {
362 // calculate data element and gradient
363 if (fgFunc != nullptr)
364 fval = fgFunc->DataElement(x, i, &gf[0]);
365 else
366 fval = fgGradFunc->DataElement(x, i, &gf[0]);
367 }
368 sum -= fval;
369
370 for (unsigned int j = 0; j < npar; ++j) {
371 double gfj = gf[j]; // / fval;
372 grad[j] -= gfj;
373 for (unsigned int k = j; k < npar; ++k) {
374 int idx = j + k * (k + 1) / 2;
375 hess[idx] += gfj * gf[k]; // / (fval );
376 }
377 }
378 }
379 } else {
380 Error("EvaluateFCN", " type of fit method is not supported, it must be chi2 or log-likelihood");
381 }
382
383 // now TFumili excludes fixed parameter in second-derivative matrix
384 // ned to get them using the static instance of TFumili
385 double * zmatrix = fgFumili->GetZ();
386 double * pl0 = fgFumili->GetPL0(); // parameter limits
387 assert(zmatrix != nullptr);
388 assert(pl0 != nullptr);
389 unsigned int k = 0;
390 unsigned int l = 0;
391 for (unsigned int i = 0; i < npar; ++i) {
392 for (unsigned int j = 0; j <= i; ++j) {
393 if (pl0[i] > 0 && pl0[j] > 0) { // only for non-fixed parameters
394 zmatrix[l++] = hess[k];
395 }
396 k++;
397 }
398 }
399
400#ifdef DEBUG
401 std::cout << "FCN value " << sum << " grad ";
402 for (unsigned int ipar = 0; ipar < npar; ++ipar)
403 std::cout << grad[ipar] << "\t";
404 std::cout << std::endl << std::endl;
405#endif
406
407 // evaluate directly function value in case of Poisson likelihoods
408 if (fgFunc && fgFunc->Type() == ROOT::Math::FitMethodFunction::kPoissonLikelihood) return 0.5 * (*fgFunc)(x);
410 return sum;
411
412}
413
414
415
416bool TFumiliMinimizer::SetVariable(unsigned int ivar, const std::string & name, double val, double step) {
417 // set a free variable.
418 if (fFumili == nullptr) {
419 Error("SetVariableValue","invalid TFumili pointer. Set function first ");
420 return false;
421 }
422#ifdef DEBUG
423 std::cout << "set variable " << ivar << " " << name << " value " << val << " step " << step << std::endl;
424#endif
425
426 int ierr = fFumili->SetParameter(ivar , name.c_str(), val, step, 0., 0. );
427 if (ierr) {
428 Error("SetVariable","Error for parameter %d ",ivar);
429 return false;
430 }
431 return true;
432}
433
434bool TFumiliMinimizer::SetLimitedVariable(unsigned int ivar, const std::string & name, double val, double step, double lower, double upper) {
435 // set a limited variable.
436 if (fFumili == nullptr) {
437 Error("SetVariableValue","invalid TFumili pointer. Set function first ");
438 return false;
439 }
440#ifdef DEBUG
441 std::cout << "set limited variable " << ivar << " " << name << " value " << val << " step " << step << std::endl;
442#endif
443 int ierr = fFumili->SetParameter(ivar, name.c_str(), val, step, lower, upper );
444 if (ierr) {
445 Error("SetLimitedVariable","Error for parameter %d ",ivar);
446 return false;
447 }
448 return true;
449}
450#ifdef LATER
451bool Fumili2Minimizer::SetLowerLimitedVariable(unsigned int ivar , const std::string & name , double val , double step , double lower ) {
452 // add a lower bounded variable as a double bound one, using a very large number for the upper limit
453 double s = val-lower;
454 double upper = s*1.0E15;
455 if (s != 0) upper = 1.0E15;
456 return SetLimitedVariable(ivar, name, val, step, lower,upper);
457}
458#endif
459
460
461bool TFumiliMinimizer::SetFixedVariable(unsigned int ivar, const std::string & name, double val) {
462 // set a fixed variable.
463 if (fFumili == nullptr) {
464 Error("SetVariableValue","invalid TFumili pointer. Set function first ");
465 return false;
466 }
467
468
469 int ierr = fFumili->SetParameter(ivar, name.c_str(), val, 0., val, val );
470 fFumili->FixParameter(ivar);
471
472#ifdef DEBUG
473 std::cout << "Fix variable " << ivar << " " << name << " value " << std::endl;
474#endif
475
476 if (ierr) {
477 Error("SetFixedVariable","Error for parameter %d ",ivar);
478 return false;
479 }
480 return true;
481}
482
483bool TFumiliMinimizer::SetVariableValue(unsigned int ivar, double val) {
484 // set the variable value
485 if (fFumili == nullptr) {
486 Error("SetVariableValue","invalid TFumili pointer. Set function first ");
487 return false;
488 }
489 TString name = fFumili->GetParName(ivar);
490 double oldval, verr, vlow, vhigh = 0;
491 int ierr = fFumili->GetParameter( ivar, &name[0], oldval, verr, vlow, vhigh);
492 if (ierr) {
493 Error("SetVariableValue","Error for parameter %d ",ivar);
494 return false;
495 }
496#ifdef DEBUG
497 std::cout << "set variable " << ivar << " " << name << " value "
498 << val << " step " << verr << std::endl;
499#endif
500
501 ierr = fFumili->SetParameter(ivar , name , val, verr, vlow, vhigh );
502 if (ierr) {
503 Error("SetVariableValue","Error for parameter %d ",ivar);
504 return false;
505 }
506 return true;
507}
508
510 // perform the minimization using the algorithm chosen previously by the user
511 // By default Migrad is used.
512 // Return true if the found minimum is valid and update internal cached values of
513 // minimum values, errors and covariance matrix.
514
515 if (fFumili == nullptr) {
516 Error("SetVariableValue","invalid TFumili pointer. Set function first ");
517 return false;
518 }
519
520 // need to set static instance to be used when calling FCN
521 fgFumili = fFumili;
522
523
524 double arglist[10];
525
526 // error cannot be set in TFumili (always the same)
527// arglist[0] = ErrorUp();
528// fFumili->ExecuteCommand("SET Err",arglist,1);
529
530 int printlevel = PrintLevel();
531
532 //arglist[0] = printlevel - 1;
533 //fFumili->ExecuteCommand("SET PRINT",arglist,1,ierr);
534
535 // suppress warning in case Printlevel() == 0
536 if (printlevel <= 0)
537 fFumili->ExecuteCommand("SET NOW",arglist,0);
538 else
539 fFumili->ExecuteCommand("SET WAR",arglist,0);
540
541 if (printlevel < 3)
542 fFumili->ExecuteCommand("SET NOD",arglist,0);
543 else
544 fFumili->ExecuteCommand("SET DEB",arglist,0);
545
546 // minimize: use ExecuteCommand instead of Minimize to set tolerance and maxiter
547
548 arglist[0] = MaxFunctionCalls();
549 arglist[1] = Tolerance();
550
551 if (printlevel > 0)
552 std::cout << "Minimize using TFumili with tolerance = " << Tolerance()
553 << " max calls " << MaxFunctionCalls() << std::endl;
554
555 int iret = fFumili->ExecuteCommand("MIGRAD",arglist,2);
556 fStatus = iret;
557 //int iret = fgFumili->Minimize();
558
559 // Hesse and IMP not implemented
560// // run improved if needed
561// if (ierr == 0 && fType == ROOT::Fumili::kMigradImproved)
562// fFumili->mnexcm("IMPROVE",arglist,1,ierr);
563
564// // check if Hesse needs to be run
565// if (ierr == 0 && IsValidError() ) {
566// fFumili->mnexcm("HESSE",arglist,1,ierr);
567// }
568
569
570 int ntot;
571 int nfree;
572 double errdef = 0; // err def is not used by Fumili
573 fFumili->GetStats(fMinVal,fEdm,errdef,nfree,ntot);
574
575 // recompute function value (because in case of likelihood is not correct)
576
577
578 if (printlevel > 0)
579 fFumili->PrintResults(printlevel,fMinVal);
580
581
582 assert (static_cast<unsigned int>(ntot) == fDim);
583 assert( nfree == fFumili->GetNumberFreeParameters() );
584 fNFree = nfree;
585
586
587 // get parameter values and correlation matrix
588 // fumili stores only lower part of diagonal matrix of the free parameters
589 fParams.resize( fDim);
590 fErrors.resize( fDim);
591 fCovar.resize(fDim*fDim);
592 const double * cv = fFumili->GetCovarianceMatrix();
593 unsigned int l = 0;
594 for (unsigned int i = 0; i < fDim; ++i) {
595 fParams[i] = fFumili->GetParameter( i );
596 fErrors[i] = fFumili->GetParError( i );
597
598 if ( !fFumili->IsFixed(i) ) {
599 for (unsigned int j = 0; j <=i ; ++j) {
600 if ( !fFumili->IsFixed(j) ) {
601 fCovar[i*fDim + j] = cv[l];
602 fCovar[j*fDim + i] = fCovar[i*fDim + j];
603 l++;
604 }
605 }
606 }
607 }
608
609 return (iret==0) ? true : false;
610}
611
612
613// } // end namespace Fit
614
615// } // end namespace ROOT
#define f(i)
Definition RSha256.hxx:104
#define g(i)
Definition RSha256.hxx:105
#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.
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition TError.cxx:208
bool gUseFumiliFunction
char name[80]
Definition TGX11.cxx:110
Chi2FCN class for binned fits using the least square methods.
Definition Chi2FCN.h:46
class evaluating the log likelihood for binned Poisson likelihood fits it is template to distinguish ...
FitMethodFunction class Interface for objective functions (like chi2 and likelihood used in the fit) ...
Type_t
enumeration specifying the possible fit method types
Documentation for the abstract class IBaseFunctionMultiDim.
Definition IFunction.h:63
virtual unsigned int NDim() const =0
Retrieve the dimension of the function.
Param Functor class for Multidimensional functions.
1-Dim function class
Definition TF1.h:182
static ROOT::Math::FitMethodFunction * fgFunc
bool SetFixedVariable(unsigned int, const std::string &, double) override
set fixed variable (override if minimizer supports them )
TFumiliMinimizer(int dummy=0)
Default constructor (an argument is needed by plug-in manager)
bool SetLimitedVariable(unsigned int ivar, const std::string &name, double val, double step, double, double) override
set upper/lower limited variable (override if minimizer supports them )
static double EvaluateFCN(const double *x, double *g)
implementation of FCN for Fumili when user provided gradient is used
~TFumiliMinimizer() override
Destructor (no operations)
bool Minimize() override
method to perform the minimization
static ROOT::Math::FitMethodGradFunction * fgGradFunc
static TFumili * fgFumili
bool SetVariableValue(unsigned int ivar, double val) override
set the value of an existing variable
void SetFunction(const ROOT::Math::IMultiGenFunction &func) override
set the function to minimize
bool SetVariable(unsigned int ivar, const std::string &name, double val, double step) override
set free variable
static void Fcn(int &, double *, double &f, double *, int)
implementation of FCN for Fumili
void SetParNumber(Int_t ParNum)
Definition TFumili.cxx:169
Int_t SetParameter(Int_t ipar, const char *parname, Double_t value, Double_t verr, Double_t vlow, Double_t vhigh) override
Sets for parameter number ipar initial parameter value, name parname, initial error verr and limits v...
Definition TFumili.cxx:1731
Basic string class.
Definition TString.h:138
virtual void SetFCN(void(*fcn)(Int_t &, Double_t *, Double_t &f, Double_t *, Int_t))
To set the address of the minimization objective function called by the native compiler (see function...
Double_t y[n]
Definition legend1.C:17
Double_t x[n]
Definition legend1.C:17
Namespace for new Math classes and functions.
T EvalLog(T x)
safe evaluation of log(x) with a protections against negative or zero argument to the log smooth line...
Definition Util.h:79
BasicFitMethodFunction< ROOT::Math::IMultiGenFunction > FitMethodFunction
Definition Fitter.h:45
TLine l
Definition textangle.C:4
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2339